Пост 2 для начинающих посвящен описательным статистикам, группированию данных и нормальному распределению. Все эти сведения заложат основу для дальнейшего анализа электоральных данных.

Описательные статистики

Описательные статистические величины, или статистики, это числа, которые используются для обобщения и описания данных. В целях демонстрации того, что мы имеем в виду, посмотрим на столбец с данными об электорате Electorate. Он показывает суммарное число зарегистрированных избирателей в каждом избирательном округе:

def ex_1_6(): '''Число значений в поле "Электорат"''' return load_uk_scrubbed()['Electorate'].count()

650

Мы уже очистили столбец, отфильтровав пустые значения (nan) из набора данных, и поэтому предыдущий пример должен вернуть суммарное число избирательных округов.

Описательные статистики, которые также называются сводными, представляют собой способы измерения свойств последовательностей чисел. Они помогают охарактеризовать последовательность и способны выступать в качестве ориентира для дальнейшего анализа. Начнем с двух самых базовых статистик, которые мы можем вычислить из последовательности чисел ее среднее значение и дисперсию.

• Среднее значение

Наиболее распространенный способ усреднить набор данных взять его среднее значение. Среднее значение на самом деле представляет собой один из нескольких способов измерения центра распределения данных. Среднее значение числового ряда вычисляется на Python следующим образом:

def mean(xs):  '''Среднее значение числового ряда''' return sum(xs) / len(xs) 

Мы можем воспользоваться нашей новой функцией mean для вычисления среднего числа избирателей в Великобритании:

def ex_1_7(): '''Вернуть среднее значение поля "Электорат"''' return mean( load_uk_scrubbed()['Electorate'] )

70149.94

На самом деле, библиотека pandas уже содержит функцию mean, которая гораздо эффективнее вычисляет среднее значение последовательности. В нашем случае ее можно применить следующим образом:

load_uk_scrubbed()['Electorate'].mean()

• Медиана

Медиана это еще одна распространенная описательная статистика для измерения центра распределения последовательности. Если Вы упорядочили все данные от меньшего до наибольшего, то медиана это значение, которое находится ровно по середине. Если в последовательности число точек данных четное, то медиана определяется, как полусумма двух срединных значений.

def median(xs): '''Медиана числового ряда''' n = len(xs) mid = n // 2 if n % 2 == 1: return sorted(xs)[mid] else: return mean( sorted(xs)[mid-1:][:2] )

Медианное значение электората Великобритании составляет:

def ex_1_8(): '''Вернуть медиану поля "Электорат"''' return median( load_uk_scrubbed()['Electorate'] )

70813.5

Библиотека pandas тоже располагает встроенной функцией для вычисления медианного значения, которая так и называется median.

• Дисперсия

Среднее арифметическое и медиана являются двумя альтернативными способами описания среднего значения последовательности, но сами по себе они мало что говорят о содержащихся в ней значениях. Например, если известно, что среднее последовательности из девяноста девяти значений равно 50, то мы почти ничего не скажем о том, какого рода значения последовательность содержит.

Она может содержать целые числа от 1 до 99 или два ряда чисел, состоящих из 49 нулей и 50 девяносто-девяток, а может быть и так, что она содержит ряд из 98 чисел, равных -1 и одно единственное значение 5048, или же вообще все значения могут быть равны 50.

Дисперсия последовательности чисел показывает "разброс" данных вокруг среднего значения. К примеру, данные, приведенные выше, имели бы разную дисперсию. На языке математики дисперсия обозначается следующим образом:

где s² это математический символ, который часто используют для обозначения дисперсии.

Выражение

можно преобразовать в функцию квадрата отклонения, square_deviation, которую можно применить к каждому элементу последовательности xs. При этом мы также воспользуемся созданной ранее функцией вычисления среднего значения mean, которая суммирует значения последовательности и делит полученную сумму на их количество.

def variance(xs): '''Дисперсия числового ряда, несмещенная дисперсия при n <= 30''' mu = mean(xs) n = len(xs) n = n-1 if n in range(1, 30) else n  square_deviation = lambda x : (x - mu) ** 2  return sum( map(square_deviation, xs) ) / n

Для вычисления квадрата выражения используется оператор языка Python возведения в степень **.

Поскольку мы взяли средний квадрат отклонения, т.е. получили квадрат отклонения и затем его среднее, то единицы измерения дисперсии тоже будут в квадрате, т.е. дисперсия электората Великобритании будет измеряться "людьми в квадрате". Несколько неестественно рассуждать об избирателях в таком виде. Единицу измерения можно привести к более естественному виду, снова обозначающему "людей", путем извлечения квадратного корня из дисперсии. В результате получим так называемое стандартное отклонение, или среднеквадратичное отклонение:

def standard_deviation(xs): '''Стандартное отклонение числового ряда''' return sp.sqrt( variance(xs) ) def ex_1_9(): '''Стандартное отклонение поля "Электорат"''' return standard_deviation( load_uk_scrubbed()['Electorate'] )

7672.77

В библиотеке pandas функции для вычисления дисперсии и стандартного отклонения имплементированы соответственно, как var и std. При этом последняя по умолчанию вычисляет несмещенное значение, поэтому, чтобы получить тот же самый результат, нужно применить именованный аргумент ddof=0, который сообщает, что требуется вычислить смещенное значение стандартного отклонения:

load_uk_scrubbed()['Electorate'].std( ddof=0 )

• Квантили

Медиана это один из способов вычислить из списка срединное значение, т.е. находящееся ровно по середине, дисперсия же предоставляет способ измерить разброс данных вокруг среднего значения. Если весь разброс данных представить на шкале от 0 до 1, то значение 0.5 будет медианным.

Для примера рассмотрим следующую ниже последовательность чисел:

[10 11 15 21 22.5 28 30]

Отсортированная последовательность состоит из семи чисел, поэтому медианой является число 21 четвертое в ряду. Его также называют 0.5-квантилем. Мы можем получить более полную картину последовательности чисел, взглянув на 0.0 (нулевой), 0.25, 0.5, 0.75 и 1.0 квантили. Все вместе эти цифры не только показывают медиану, но также обобщают диапазон данных и сообщат о характере распределения чисел внутри него. Они иногда упоминаются в связи с пятичисловой сводкой.

Один из способов составления пятичисловой сводки для данных об электорате Великобритании показан ниже. Квантили можно вычислить непосредственно в pandas при помощи функции quantile. Последовательность требующихся квантилей передается в виде списка.

def ex_1_10(): '''Вычислить квантили: возвращает значение в последовательности xs,  соответствующее p-ому проценту''' q = [0, 1/4, 1/2, 3/4, 1] return load_uk_scrubbed()['Electorate'].quantile(q=q)

0.00 21780.000.25 65929.250.50 70813.500.75 74948.501.00 109922.00Name: Electorate, dtype: float64

Когда квантили делят диапазон на четыре равных диапазона, как показано выше, то они называются квартилями. Разница между нижним (0.25) и верхним (0.75) квартилями называется межквартильным размахом, или иногда сокращенно МКР. Аналогично дисперсии вокруг среднего значения, межквартильный размах измеряет разброс данных вокруг медианы.

Группирование данных в корзины

В целях развития интуитивного понимания в отношении того, что именно все эти расчеты разброса значений измеряют, мы можем применить метод под названием группировка в частотные корзины (binning). Когда данные имеют непрерывный характер, использование специального словаря для подсчета частот Counter (подобно тому, как он использовался при подсчете количества пустых значений в наборе данных об электорате) становится нецелесообразным, поскольку никакие два значения не могут быть одинаковыми. Между тем, общее представление о структуре данных можно все-равно получить, сгруппировав для этого данные в частотные корзины (bins).

Процедура образования корзин заключается в разбиении диапазона значений на ряд последовательных, равноразмерных и меньших интервалов. Каждое значение в исходном ряду попадает строго в одну корзину. Подсчитав количества точек, попадающих в каждую корзину, мы можем получить представление о разбросе данных:

На приведенном выше рисунке показано 15 значений x, разбитых на 5 равноразмерных корзин. Подсчитав количество точек, попадающих в каждую корзину, мы можем четко увидеть, что большинство точек попадают в корзину по середине, а меньшинство в корзины по краям. Следующая ниже функция Python nbin позволяет добиться того же самого результата:

def nbin(n, xs):  '''Разбивка данных на частотные корзины''' min_x, max_x = min(xs), max(xs) range_x = max_x - min_x fn = lambda x: min( int((abs(x) - min_x) / range_x * n), n-1 ) return map(fn, xs)

Например, мы можем разбить диапазон 0-14 на 5 корзин следующим образом:

list( nbin(5, range(15)) )

[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4]

После того, как мы разбили значения на корзины, мы можем в очередной раз воспользоваться словарем Counter, чтобы подсчитать количество точек в каждой корзине. В следующем ниже примере мы воспользуемся этим словарем для разбиения данных об электорате Великобритании на пять корзин:

def ex_1_11(): '''Разбиmь электорат Великобритании на 5 корзин''' series = load_uk_scrubbed()['Electorate'] return Counter( nbin(5, series) )

Counter({2: 450, 3: 171, 1: 26, 0: 2, 4: 1})

Количество точек в крайних корзинах (0 и 4) значительно ниже, чем в корзинах в середине количества, судя по всему, растут по направлению к медиане, а затем снова снижаются. В следующем разделе мы займемся визуализацией формы этих количеств.

Гистограммы

Гистограмма это один из способов визуализации распределения одной последовательности значений. Гистограммы попросту берут непрерывное распределение, разбивают его на корзины, и изображают частоты точек, попадающих в каждую корзину, в виде столбцов. Высота каждого столбца гистограммы показывает количество точек данных, которые содержатся в этой корзине.

Мы уже увидели, каким образом можно выполнить разбиение данных на корзины самостоятельно, однако в библиотеке pandas уже содержится функция hist, которая разбивает данные и визуализирует их в виде гистограммы.

def ex_1_12(): '''Построить гистограмму частотных корзин        электората Великобритании''' load_uk_scrubbed()['Electorate'].hist() plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

Приведенный выше пример сгенерирует следующий ниже график:

Число корзин, на которые данные разбиваются, можно сконфигурировать, передав в функцию при построении гистограммы именованный аргумент bins:

def ex_1_13(): '''Построить гистограмму частотных корзин  электората Великобритании с 200 корзинами''' load_uk_scrubbed()['Electorate'].hist(bins=200) plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

Приведенный выше график показывает единственный высокий пик, однако он выражает форму данных довольно грубо. Следующий ниже график показывает мелкие детали, но величина столбцов делает неясной форму распределения, в особенности в хвостах:

При выборе количества корзин для представления данных следует найти точку равновесия с малым количеством корзин форма данных будет представлена лишь приблизительно, а слишком большое их число приведет к тому, что шумовые признаки могут заслонить лежащую в основании структуру.

def ex_1_14(): '''Построить гистограмму частотных корзин  электората Великобритании с 20 корзинами''' load_uk_scrubbed()['Electorate'].hist(bins=20) plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

Ниже показана гистограмма теперь уже из 20 корзин:

Окончательный график, состоящий из 20 корзин, судя по всему, пока лучше всего представляет эти данные.

Наряду со средним значением и медианой, есть еще один способ измерить среднюю величину последовательности. Это мода. Мода это значение, встречающееся в последовательности наиболее часто. Она определена исключительно только для последовательностей, имеющих по меньшей мере одно дублирующееся значение; во многих статистических распределениях это не так, и поэтому для них мода не определена. Тем не менее, пик гистограммы часто называют модой, поскольку он соответствует наиболее распространенной корзине.

Из графика ясно видно, что распределение вполне симметрично относительно моды, и его значения резко падают по обе стороны от нее вдоль тонких хвостов. Эти данные приближенно подчиняются нормальному распределению.

Нормальное распределение

Гистограмма дает приблизительное представление о том, каким образом данные распределены по всему диапазону, и является визуальным средством, которое позволяет квалифицировать данные как относящиеся к одному из немногих популярных распределений. В анализе данных многие распределения встречаются часто, но ни одно не встречается также часто, как нормальное распределение, именуемое также гауссовым распределением.

^{Распределение названо нормальным распределением из-за того, что оно очень часто встречается в природе. Галилей заметил, что ошибки в его астрономических измерениях подчинялись распределению, где малые отклонения от среднего значения встречались чаще, чем большие. Вклад великого математика Гаусса в описание математической формы этих ошибок привел к тому, что это распределение стали называть в его честь распределением Гаусса.}

Любое распределение похоже на алгоритм сжатия: оно позволяет очень эффективно резюмировать потенциально большой объем данных. Нормальное распределение требует только два параметра, исходя из которых можно аппроксимировать остальные данные. Это среднее значение и стандартное отклонение.

Центральная предельная теорема

Высокая встречаемость нормального распределения отчасти объясняется центральной предельной теоремой. Дело в том, что значения, полученные из разнообразных статистических распределений, при определенных обстоятельствах имеют тенденцию сходиться к нормальному распределению, и мы это покажем далее.

В программировании типичным распределением является равномерное распределение. Оно представлено распределением чисел, генерируемых функцией библиотеки SciPy stats.uniform.rvs: в справедливом генераторе случайных чисел все числа имеют равные шансы быть сгенерированными. Мы можем увидеть это на гистограмме, многократно генерируя серию случайных чисел между 0 и 1 и затем построив график с результатами.

def ex_1_15(): '''Показать гистограмму равномерного распределения  синтетического набора данных''' xs = stats.uniform.rvs(0, 1, 10000) pd.Series(xs).hist(bins=20) plt.xlabel('Равномерное распределение') plt.ylabel('Частота') plt.show()

Обратите внимание, что в этом примере мы впервые использовали тип Series библиотеки pandas для числового ряда данных.

Приведенный выше пример создаст следующую гистограмму:

Каждый столбец гистограммы имеет примерно одинаковую высоту, что соответствует равновероятности генерирования числа, которое попадает в каждую корзину. Столбцы имеют не совсем одинаковую высоту, потому что равномерное распределение описывает теоретический результат, который наша случайная выборка не может отразить в точности. В следующих нескольких главах мы будем изучать способы точной количественной оценки расхождения между теорией и практикой, чтобы определить, являются ли расхождения достаточно большими, чтобы обратить на это внимание. В данном случае они таковыми не являются.

Если напротив сгенерировать гистограмму средних значений последовательностей чисел, то в результате получится распределение, которое выглядит совсем непохоже.

def bootstrap(xs, n, replace=True):  '''Вернуть список массивов меньших размеров  по n элементов каждый''' return np.random.choice(xs, (len(xs), n), replace=replace) def ex_1_16(): '''Построить гистограмму средних значений''' xs = stats.uniform.rvs(loc=0, scale=1, size=10000) pd.Series( map(sp.mean, bootstrap(xs, 10)) ).hist(bins=20) plt.xlabel('Распределение средних значений')  plt.ylabel('Частота') plt.show()

Приведенный выше пример сгенерирует результат, аналогичный следующей ниже гистограмме:

Хотя величина среднего значения близкая к 0 или 1 не является невозможной, она является чрезвычайно невероятной и становится менее вероятной по мере роста числа усредненных чисел и числа выборочных средних. Фактически, на выходе получается результат очень близкий к нормальному распределению.

Этот результат, когда средний эффект множества мелких случайных колебаний в итоге приводит к нормальному распределению, называется центральной предельной теоремой, иногда сокращенно ЦПТ, и играет важную роль для объяснения, почему нормальное распределение встречается так часто в природных явлениях.

До 20-ого века самого термина еще не существовало, хотя этот эффект был зафиксирован еще в 1733 г. французским математиком Абрахамом де Mуавром (Abraham de Moivre), который использовал нормальное распределение, чтобы аппроксимировать число орлов в результате бросания уравновешенной монеты. Исход бросков монеты лучше всего моделировать при помощи биномиального распределения, с которым мы познакомимся в главе 4, Классификация. В отличие от центральной предельной теоремы, которая позволяет получать выборки из приближенно нормального распределения, библиотека ScyPy содержит функции для эффективного генерирования выборок из самых разнообразных статистических распределений, включая нормальное:

def ex_1_17(): '''Показать гистограмму нормального распределения  синтетического набора данных''' xs = stats.norm.rvs(loc=0, scale=1, size=10000) pd.Series(xs).hist(bins=20) plt.xlabel('Нормальное распределение') plt.ylabel('Частота') plt.show()

Отметим, что в функции sp.random.normal параметр loc это среднее значение, scale дисперсия и size размер выборки. Приведенный выше пример сгенерирует следующую гистограмму нормального распределения:

По умолчанию среднее значение и стандартное отклонение для получения нормального распределения равны соответственно 0 и 1.

Следующая часть, часть 3, серии постов Python, наука о данных и выборы посвящена генерированию распределений, их свойствам, а также графикам для их сопоставительного анализа

Заключительный пост 5 для начинающих посвящен сопоставительной визуализации электоральных данных.

Сопоставительная визуализация электоральных данных

Теперь рассмотрим набор данных других всеобщих выборов, на этот раз Российских, проходивших в 2011 г. Россия гораздо более крупная страна, и поэтому данные о проголосовавших на выборах там гораздо объемнее. Для этого мы загрузим в оперативную память один большой TSV-файл с разделением полей данных символом табуляции.

def load_ru(): '''Загрузить данные по России''' return pd.read_csv('data/ch01/Russia2011.tsv', '\t')

Посмотрим, какие имена столбцов имеются в российских данных:

def ex_1_29(): '''Показать список полей электоральных  данных по России''' return load_ru().columns

Будет выведен следующий список столбцов:

Index(['Код ОИК', 'ОИК ', 'Имя участка','Число избирателей, внесенных в список избирателей',...'Политическая партия СПРАВЕДЛИВАЯ РОССИЯ','Политическая партия ЛДПР - Либерально-демократическая партия России','Политическая партия "ПАТРИОТ РОССИИ"','Политическая партия КОММУНИСТИЧЕСКАЯ ПАРТИЯ КОММУНИСТ РОССИИ','Политическая партия "Российская объединенная демократическая партия "ЯБЛОКО"','Политическая партия "ЕДИНАЯ РОССИЯ"','Всероссийская политическая партия "ПАРТИЯ РОСТА"'],dtype='object')

Имена столбцов в российском наборе данных очень описательны, но, пожалуй, длиннее, чем нужно. Также, было бы удобно, если столбцы, представляющие те же самые атрибуты, которые мы уже видели в данных по выборам в Великобритании (к примеру, доля победителя и явка на выборы) были промаркированы одинаковым образом в обоих наборах данных. Переименуем их надлежащим образом.

Наряду с набором данных функция библиотеки pandas rename ожидает словарь, в котором ключам с текущими именами столбцов поставлены в соответствие значения с новыми именами. Если объединить ее с данными, которые мы уже рассматривали, то мы получим следующее:

def load_ru_victors(): '''Загрузить данные по России,  выбрать, переименовать и вычислить поля''' new_cols_dict = { 'Число избирателей, внесенных в список избирателей':'Электорат', 'Число действительных избирательных бюллетеней': 'Действительные бюллетени', 'Политическая партия "ЕДИНАЯ РОССИЯ"':'Победитель'  } newcols = list(new_cols_dict.values())  df = load_ru().rename( columns=new_cols_dict )[newcols]  df['Доля победителя'] = df['Победитель'] / df['Действительные бюллетени']  df['Явка'] = df['Действительные бюллетени'] / df['Электорат']  return df

Библиотека pandas располагает функцией безопасного деления divide, которая идентична операции /, но защищает от деления на ноль. Она вместо пропущенного значения (nan) в одном из полей подставляет значение, передаваемое в именованном аргументе fill_value. Если же оба значения поля равны nan, то результат будет отсутствовать. Поэтому операцию деления можно было бы переписать следующим образом:

 df[ 'Доля победителя' ] = \ df[ 'Победитель' ].divide( df[ 'Действительные бюллетени' ], \ fill_value=1 )

Визуализация электоральных данных РФ

Мы ранее видели, что гистограмма явки на выборы в Великобритании была приближенно нормально распределенной (хотя и с легкими хвостами). Теперь, когда мы загрузили и преобразовали данные о выборах в России, посмотрим, насколько они сопоставимы:

def ex_1_30(): '''Показать гистограмму  электоральных данных по России''' load_ru_victors()['Явка'].hist(bins=20) plt.xlabel('Явка в России')  plt.ylabel('Частота') plt.show()

Приведенный выше пример сгенерирует следующую гистограмму:

Эта гистограмма совсем не похожа на классические колоколообразные кривые, которые мы видели до сих пор. Имеется явно выраженная положительная асимметрия, и явка избирателей в действительности увеличивается с 80% в сторону 100% совсем не то, что мы ожидали бы от нормально распределенных данных.

Учитывая ожидания, заданные данными из Британии и центральной предельной теоремой (ЦПТ), такой результат любопытен. Для начала покажем данные на квантильном графике:

def ex_1_31(): '''Показать квантильный график  победителя на выборах в РФ''' qqplot( load_ru_victors()['Доля победителя'].dropna() ) plt.show()

Этот пример вернет следующий график:

На квантильном графике показана линия, которая не является ни прямой, ни одной из S-образных кривых. По существу, квантильный график говорит о наличии легкого хвоста в верхнем конце распределения и тяжелого хвоста в нижнем. Это почти противоположно тому, что мы видим на гистограмме, которая четко указывает на крайне тяжелый правый хвост.

На самом деле, этот квантильный график дезориентирует, и происходит этот именно потому, что хвост очень тяжелый: плотность точек между 0.5 и 1.0 на гистограмме говорит о том, что пик должен составлять порядка 0.7 с последующим правым хвостом за пределами 1.0. Наличие значения, превышающего 100% явно выходит за рамки логики, но квантильный график не объясняет это (он не учитывает, что речь идет о процентах), так что внезапное отсутствие данных за пределами 1.0 интерпретируется как подрезанный правый хвост.

С учетом центральной предельной теоремы и того, что мы наблюдали в данных выборов в Великобритании, тенденция к 100% явке избирателей на выборы выглядит очень любопытно. Давайте выполним параллельный сопоставительный анализ наборов данных по Великобритании и России.

Сравнительная визуализация

Предположим, мы хотели бы сравнить распределение электоральных данных между Великобританией и Россией. Мы уже видели в этой главе, как использовать ИФР и коробчатые диаграммы, поэтому теперь займемся исследованием альтернативного варианта, который аналогичен гистограмме.

Мы могли бы попытаться изобразить оба набора данных на гистограмме, но это будет безуспешной затеей, поскольку результаты не поддаются интерпретации по двум следующим причинам:

• Размеры избирательных округов, и, следовательно, средних значений распределений сильно отличаются

• Абсолютные количества избирательных округов настолько отличаются, что столбцы гистограмм будут иметь разную высоту

Вместо гистограммы альтернативным вариантом, позволяющим решить обе эти задачи, является функция массы вероятности.

Функции массы вероятности

Функция массы вероятности (ФМВ), от англ. Probability Mass Function (PMF), чаще именуемая функцией вероятности дискретной случайной величины, имеет много общего с гистограммой. Однако, вместо того, чтобы показывать количества значений, попадающих в группы, она показывает вероятность, что взятое из распределения число будет в точности равно заданному значению. Поскольку функция закрепляет вероятность за каждым значением, которое может быть возвращено распределением, и поскольку вероятности измеряются по шкале от 0 до 1, (где 1 соответствует полной определенности), то площадь под функцией массы вероятности равна 1.

Таким образом функция массы вероятности гарантирует, что площадь под нашими графиками будет между наборами данных сопоставима. Однако у нас все еще имеется одно затруднение, которое заключается в том, что размеры избирательных округов и поэтому средних значений распределений несопоставимы. Это затруднение решается отдельно при помощи нормализации.

Существует неисчислимое количество способов нормализации данных, однако один из самых основных обеспечивает, чтобы каждый числовой ряд находился в диапазоне от 0 до 1. Ни одно наше значение не находится в отрицательном диапазоне, поэтому мы можем выполнить нормализацию, попросту разделив каждое индивидуальное значение на самое большое:

def plot_as_pmf(dt, label, ax): '''График функции вероятности дискретной случайной величины (или функции массы вероятности)''' s = pd.cut(dt, bins=40, labels=False) # разбить на 40 корзин pmf = s.value_counts().sort_index() / len(s) # подсчитать кво в корзинах newax = pmf.plot(label=label, grid=True, ax=ax)  return newax

Имея в распоряжении приведенную выше функцию, мы теперь можем нормализовать данные по Великобритании и России и изобразить их рядом на тех же осях:

def ex_1_32(): '''Сопоставление данных явки по Великобритании и РФ, данные нормализованы на основе функции массы вероятностей''' ax = plot_as_pmf(load_uk_victors()['Явка'], 'Великобритания', None) plot_as_pmf(load_ru_victors()['Явка'], 'Россия', ax) plt.xlabel('Интервальные группы явки') # Частотные корзины plt.ylabel('Вероятность') plt.legend(loc='best') plt.show()

Приведенный выше пример сгенерирует следующий график:

После нормализации эти два распределения вполне готовы для проведения сопоставительного анализа. Теперь становится совершенно очевидным, каким образом несмотря на более низкую среднюю явку, чем в Великобритании (0.6366 против 0.6523) на российских выборах произошел массивный подъем явки близкий к 100%. Поскольку результаты голосования представляют собой объединенный эффект многих независимых волеизъявлений, они ожидаемо будут соответствовать центральной предельной теореме и будут приближенно нормально распределенными. В сущности, за редким исключением, как в Канаде, например, где население имеет гетерогенный характер (там французскоговорящая и англоговорящая группы населения в результате дают бимодальную кривую), результаты выборов по всему миру такому ожиданию обычно соответствуют.

Данные российских выборов показывают чрезвычайно аномальный результат, хотя и не настолько высокий, как модальный пик в центре распределения, который приблизительно соответствует 50% явке. Исследователь Питер Климек (Peter Klimek) и его коллеги в Венском медицинском университете пошли дальше и предположили, что этот результат является явным признаком подтасовки результатов голосования.

Диаграммы рассеяния

Мы обнаружили любопытные результаты, связанные с явкой на российских выборах и установили, что они имеют сигнатуру, отличающуюся от британских выборов. Теперь посмотрим, каким образом доля голосов за побеждающего кандидата связана с явкой. В конце концов, если неожиданно высокая явка действительно является признаком нечестной игры в пользу действующего президента и правительства, то можно ожидать, что они будут голосовать за себя, а не за кого-либо еще. Таким образом, большинство, если не все, из этих дополнительных голосов ожидаемо будут отданы за итоговых победителей выборов.

Раздел статистики посвященный корреляции довольно подробно рассматривает величины, лежащие в основе взаимосвязи двух переменных, однако на данный момент было бы интересно попросту визуализировать связь между явкой и долей голосов за побеждающую партию.

Заключительный технический прием визуализации, с которым мы познакомим в этой главе, представлен диаграммой рассеяния. Диаграммы рассеяния очень хорошо подходят для визуализации взаимосвязей между двумя переменными: там, где существует линейная взаимосвязь, на графике она будет видна, как диагональная направленность. Библиотека pandas содержит для этого вида графиков функцию scatter с такими же аргументами, что и для функции двумерных графиков plot.

def ex_1_33(): '''Показать диаграмму рассеяния  выборов в Великобритании''' df = load_uk_victors()[ ['Явка', 'Доля победителей'] ] df.plot.scatter(0, 1, s=3) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.show()

Приведенный выше пример сгенерирует следующую ниже диаграмму:

Хотя точки широко разбросаны в виде нечеткого эллипса, четко видна общая диагональная направленность к верхнему правому углу графика рассеяния. Она указывает на интересный результат явка на выборы коррелирует с долей голосов, отданных за окончательных победителей на выборах, в то время, как ожидалось получить обратный результат: наличия так называемого самодовольства избирателей, которое приводит к более низкой явке в случае, когда победитель гонки очевиден.

^{Как отмечалось ранее, британские выборы 2010 г. были далеко необычными: они привели к "подвисшему" парламенту и коалиционному правительству. Фактически, "победители" в данном случае представлены обеими сторонами, которые были противниками, вплоть до дня выборов. И поэтому голосование за любую из партий считается как голосование за победителя.}

Затем, мы создадим такую же диаграмму рассеяния для выборов в России:

def ex_1_34(): '''Показать диаграмму рассеяния выборов в РФ''' df = load_ru_victors()[ ['Явка', 'Доля победителя'] ] df.plot.scatter(0, 1, s=3) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.show()

Этот пример сгенерирует следующую диаграмму:

Хотя из расположения точек в данных по России четко прослеживается диагональная направленность, сплошной объем данных затеняет внутреннюю структуру. В последнем разделе этой главы мы покажем простой метод, который позволяет с помощью настройки прозрачности графического изображения вычленять структуру из графиков, подобных этому.

Настройка прозрачности рассеяния

В ситуациях, аналогичных приведенной выше, когда диаграмма рассеяния переполнена точками, можно улучшить визуализацию внутренней структуры данных путем настройки уровня прозрачности. Поскольку наложение полупрозрачных точек будет менее прозрачным, а области с меньшим количеством точек будут более прозрачными, то диаграмма рассеяния с полупрозрачными точками может изобразить плотность данных намного лучше, чем сплошные точки.

Выполнить настройку альфа-канала, регулирующего прозрачность изображаемых на графике pandas точек можно при помощи именованного аргумента alpha в функции scatter в виде числа между 0 и 1, где 1 означает полную непрозрачность, 0 полную прозрачность.

def ex_1_35(): '''Показать диаграмму рассеяния (с прозрачностью) выборов в РФ''' df = load_ru_victors()[ ['Явка', 'Доля победителя'] ] rows = sp.random.choice(df.index.values, 10000) df.loc[rows].plot.scatter(0, 1, s=3, alpha=0.1) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.axis([0, 1.05, 0, 1.05]) plt.show()

Приведенный выше пример сгенерирует следующую диаграмму:

Приведенная выше диаграмма рассеяния показывает общую направленность совместного изменения доли победителя и явки на выборы. Мы видим корреляцию между двумя значениями и "горячую точку" в правом верхнем углу графика, которая соответствует явке близкой к 100% и 100%-ому голосованию в пользу побеждающей стороны. Как раз эта особенность в частности является признаком того, что исследователи из Венского медицинского университета обозначили как сигнатура фальсификации выборов. Этот факт также подтверждается результатами других спорных выборов по всему миру, например, таких как президентские выборы 2011 г. в Уганде.

^{Результаты многих других выборов по всему миру на уровне округов доступны здесь. На указанном веб-сайте можно получить ссылку на научно-исследовательскую работу и скачать другие наборы данных, на которых можно попрактиковать то, что Вы узнали в этой главе об исправлении и преобразовании реальных данных.}

Примеры исходного кода для этого поста находится в моем репо на Github.

Выводы

Эта серия постов была посвящена сводным статистическим величинам и значениям статистических распределений. Мы увидели, каким образом даже простой анализ может предъявить свидетельства о потенциальной фальсификационной активности.

В частности, мы познакомились с центральной предельной теоремой и причиной, почему она играет такую важную роль в объяснении широкого применения нормального распределения в науке о данных. Подходящее статистическое распределение способно всего в нескольких статистиках выразить сущность большой последовательности чисел, некоторые из которых были имплементированы тут на основе встроенных функций языка Python и функций библиотеки scipy. Кроме того, были показаны возможности библиотеки pandas, которая была задействована для загрузки, преобразования и визуального сопоставления нескольких наборов данных. В ходе анализа было обнаружено любопытное расхождение между двумя статистическими распределениями.

Внизу поста можно проголосовать за или против размещения следующей серии постов. Их тема - инференциальная статистика, которая позволяет давать количественную оценку измеренному расхождению между двумя или несколькими статистическими распределениями и принимать решение о статистической значимости этого расхождения. Указанная серия постов также посвящена методам проверки статистических гипотез, которые лежат в основе проведения робастных статистических экспериментов, позволяющих делать выводы на основе имеющихся данных.

Пост 4 для начинающих посвящен техническим приемам визуализации данных.

Важность визуализации

Простые приемы визуализации, подобные тем, которые были показаны ранее, позволяют лаконично передать большое количество информации. Они дополняют сводные статистики, которые мы рассчитали ранее в этой главе, и поэтому очень важно уметь ими пользоваться. Такие статистики, как среднее значение и стандартное отклонение, неизбежно скрывают много информации по той причине, что сворачивают последовательность в одно единственное число.

Английский математик Фрэнсис Энскомб составил коллекцию из четырех точечных графиков, ныне известную как квартет Энскомба, которые обладают практически идентичными статистическими свойствами (включая среднее, дисперсию и стандартное отклонение). Несмотря на это, они четко показывают, что распределение значений последовательностей и сильно расходится:

При этом наборы данных не должны быть подстроенными, потому как ценные аналитические выводы будут непременно выявлены при построении графиков. Возьмем для примера гистограмму оценок, полученных выпускниками по национальному экзамену на аттестат зрелости в Польше в 2013 г.:

Способности учащихся вполне ожидаемо должны быть нормально распределены, и действительно за исключением крутого всплеска вокруг 30% так и есть. То, что мы ясно наблюдаем это результат по-человечески понятного натягивания экзаменаторами оценок учащегося на проходной балл.

На самом деле статистические распределения для последовательностей, взятых из крупных выборок, могут быть настолько надежными, что даже незначительное отклонение от них может являться свидетельством противоправной деятельности. Закон Бенфорда, или закон первой цифры, показывает любопытную особенность случайных чисел, которые генерируются в широком диапазоне. Единица появляется в качестве ведущей цифры примерно в 30% случаев, в то время как цифры крупнее появляется все реже и реже. Например, девятка появляется в виде первой цифры менее чем в 5% случаев.

^{Закон Бенфорда назван в честь физика Фрэнка Бенфорда (Frank Benford), который сформулировал его в 1938 г., показав его состоятельность на различных источниках данных. Проявление этого закона было ранее отмечено американским астрономом Саймоном Ньюкомом (Simon Newcomb), который еще более 50 лет назад до него обратил внимание на страницы своих логарифмических справочников: страницы с номерами, начинавшихся с цифры 1, имели более потрепанный вид.}

Бенфорд показал, что его закон применяется к самым разнообразным данным, таким как счета за электричество, адреса домов, цены на акции, численность населения, уровень смертности и длины рек. Для наборов данных, которые охватывают большие диапазоны значений, закон настолько состоятелен, что случаи отклонения от него стали принимать в качестве подтверждающих данных в судебной практике по финансовым махинациям.

Визуализация данных об электорате

Вернемся к данным выборов и сравним электоральную последовательность, которую мы создали ранее, относительно теоретической нормальной ИФР. Для создания нормальной ИФР из последовательности значений можно воспользоваться функцией sp.random.normal библиотеки SciPy, как уже было показано выше. Среднее значение и стандартное отклонение по умолчанию равны соответственно 0 и 1, поэтому нам нужно предоставить измеренные среднее значение и стандартное отклонение, взятые из электоральных данных. Эти значения для наших электоральных данных составляют соответственно 70150 и 7679.

Ранее в этой главе мы уже генерировали эмпирическую ИФР. Следующий ниже пример просто сгенерирует обе ИФР и выведет их на одном двумерном графике:

def ex_1_24(): '''Показать эмпирическую и подогнанную ИФР  электората Великобритании''' emp = load_uk_scrubbed()['Electorate'] fitted = stats.norm.rvs(emp.mean(), emp.std(ddof=0), len(emp)) df = empirical_cdf(emp) df2 = empirical_cdf(fitted) ax = df.plot(0, 1, label='эмпирическая')  df2.plot(0, 1, label='подогнанная', grid=True, ax=ax)  plt.xlabel('Электорат') plt.ylabel('Вероятность') plt.legend(loc='best') plt.show()

Приведенный выше пример создаст следующий график:

Несмотря на наличие незначительной асимметрии, близкая расположенность двух кривых друг к другу говорит о нормальности исходных данных. Асимметрия выражена в противоположном направлении относительно построенной ранее кривой ИФР нечестного булочника, то есть наши данные об электорате слегка смещены влево.

Поскольку мы сравниваем наше распределение с теоретическим нормальным распределением, то можно воспользоваться квантильным графиком, который делает это по умолчанию:

def ex_1_25(): '''Показать квантильный график  электората Великобритании''' qqplot( load_uk_scrubbed()['Electorate'] ) plt.show()

Следующий ниже квантильный график еще лучше показывает левую асимметрию, с очевидностью присутствующую в данных:

Как мы и ожидали, кривая изгибается в противоположном направлении по отношению к построенному ранее в этой главе квантильному графику нечестного булочника. Это свидетельствует о том, что число более мелких избирательных округов на самом деле больше, чем можно было бы ожидать, если бы данные были более тесно нормально распределены.

Добавление производных столбцов

В целях выяснения процента электората, который проголосовал за одну из двух партий, требуется вычислить сумму голосов, отданных за каждую из них. Для этого нам понадобится создать новое поле данных Victors (Победители) из данных, которые соответствуют Консервативной (Con) и Либерально-демократической (LD) партиям и заодно проверим, имеются ли пропущенные значения.

def ex_1_26(): '''Вычислить производное поле данных "Победители" и  число имеющихся в нем пропущенных значений''' df = load_uk_scrubbed() df['Победители'] = df['Con'] + df['LD'] freq = Counter(df['Con'].apply( lambda x: x > 0 )) print('Поле "Победители": %d, в т.ч. пропущено %d'  % (freq[True], freq[False]))

Поле "Победители": 631, в т.ч. пропущено 19

Результат показывает, что в 19 случаях данные отсутствуют. Очевидно, что в каком-то из столбцов: столбце Con либо столбце LD (либо обоих), данные отсутствуют, но в каком именно? Снова воспользуемся словарем Counter, чтобы увидеть масштаб проблемы:

'''Проверить пропущенные значения в полях "Консервативная партия" (Con) и   "Либерально-демократическая партия" (LD)'''df = load_uk_scrubbed()Counter(df['Con'].apply(lambda x: x > 0)),  Counter(df['LD'].apply(lambda x: x > 0))

(Counter({False: 19, True: 631}), Counter({False: 19, True: 631}))

В обоих случаях будет выведено одинаковое число строк, в которых значения есть, и в которых они отсутствуют. Воспользуемся предикативной функцией isnull, которую мы уже встречали ранее в этой главе, чтобы узнать, какие строки не содержат значений:

def ex_1_27(): '''Выборка полей данных по условию, что поля "Консервативная партия" (Con) и  "Либерально-демократическая" (LD) не пустые''' df = load_uk_scrubbed() rule = df['Con'].isnull() & df['LD'].isnull() return df[rule][['Region', 'Electorate', 'Con', 'LD']]

Небольшое обследование данных должно определить причину, почему эти поля оказались пустыми. Как оказалось, кандидаты в соответствующих избирательных округах не выдвигались. Следует ли эти строки отфильтровать или же оставить как есть, равными нулю? Это интересный вопрос. Давайте их отфильтруем, поскольку вообще-то невозможно, чтобы в этих округах избиратели выбрали какого-либо кандидата, неважно от либеральных демократов или консерваторов. Если же мы напротив допустили, что они равны нулю, то среднее количество людей, которое при заданных вариантах выбора проголосовало за одну из этих партий, было бы искусственно занижено.

Зная, как фильтровать проблемные строки, теперь добавим производные столбцы, которые будут представлять победителя, долю голосов за победителя и явку на выборы. Отфильтруем строки так, чтобы показать только те, где были выдвинуты кандидаты от обеих партий:

def load_uk_victors(): '''Загрузить данные по Великобритании,  выбрать поля и отфильтровать''' df = load_uk_scrubbed() rule = df['Con'].notnull() df = df[rule][['Con', 'LD', 'Votes', 'Electorate']]  df['Победители'] = df['Con'] + df['LD']  df['Доля победителей'] = df['Победители'] / df['Votes']  df['Явка'] = df['Votes'] / df['Electorate'] return df

В результате в нашем наборе данных теперь имеется три дополнительных столбца: Victors, Victors Share и Turnout, т.е. победители, доля победителей и явка на выборы. Покажем на квантильном графике долю голосов за победителя, чтобы увидеть, как она соотносится с теоретическим нормальным распределением:

def ex_1_28(): '''Показать квантильный график победителей  на выборах в Великобритании''' qqplot( load_uk_victors()['Доля победителей'] ) plt.show()

Приведенный выше пример создаст следующий ниже график:

Основываясь на сводном графике разных форм кривых квантильных графиков, показанном ранее в этой главе, можно заключить, что доля голосов, отданных за победителя, имеет "легкие хвосты" по сравнению с нормальным распределением. Это означает, что ближе к среднему значению расположено больше данных, чем можно было бы ожидать, исходя из действительно нормально распределенных данных.

Примеры исходного кода для этого поста находится в моем репо на Github.

Следующая заключительная часть, часть 5, серии постов Python, наука о данных и выборы посвящена сопоставительной визуализации электоральных данных.

Предыдущий пост см. здесь.

Выборки и популяции

В статистической науке термины выборка и популяция имеют особое значение. Популяция, или генеральная совокупность, это все множество объектов, которые исследователь хочет понять или в отношении которых сделать выводы. Например, во второй половине 19-го века основоположник генетики Грегор Йохан Мендель) записывал наблюдения о растениях гороха. Несмотря на то, что он изучал в лабораторных условиях вполне конкретные сорта растения, его задача состояла в том, чтобы понять базовые механизмы, лежащие в основе наследственности абсолютно всех возможных сортов гороха.

^{В статистической науке о группе объектов, из которых извлекается выборка, говорят, как о популяции, независимо от того, являются изучаемые объекты живыми существами или нет.}

Поскольку популяция может быть крупной или бесконечной, как в случае растений гороха Менделя мы должны изучать репрезентативные выборки, и на их основе делать выводы о всей популяции в целом. В целях проведения четкого различия между поддающимися измерению атрибутами выборок и недоступными атрибутами популяции, мы используем термин статистики, имея при этом в виду атрибуты выборки, и говорим о параметрах, имея в виду атрибуты популяции.

^{Статистики это атрибуты, которые мы можем измерить на основе выборок. Параметры это атрибуты популяции, которые мы пытаемся вывести статистически.}

В действительности, статистики и параметры различаются в силу применения разных символов в математических формулах:

Если вы вернетесь к уравнению стандартной ошибки, то заметите, что она вычисляется не из выборочного стандартного отклонения S_x, а из популяционного стандартного отклонения _x. Это создает парадоксальную ситуацию мы не можем вычислить выборочную статистику, используя для этого популяционные параметры, которые мы пытаемся вывести. На практике, однако, предполагается, что выборочное и популяционное стандартные отклонения одинаковы при размере выборки порядка n 30.

Теперь вычислим стандартную ошибку средних значений за определенный день. Например, возьмем конкретный день, скажем, 1 мая:

def ex_2_8():  '''Вычислить стандартную ошибку  средних значений за определенный день''' may_1 = '2015-05-01' df = with_parsed_date( load_data('dwell-times.tsv') )  filtered = df.set_index( ['date'] )[may_1] se = standard_error( filtered['dwell-time'] ) print('Стандартная ошибка:', se)

Стандартная ошибка: 3.627340273094217

Хотя мы взяли выборку всего из одного дня, вычисляемая нами стандартная ошибка очень близка к стандартному отклонению всех выборочных средних 3.6 сек. против 3.7 сек. Это, как если бы, подобно клетке, содержащей ДНК, в каждой выборке была закодирована информация обо всей находящейся внутри нее популяции.

Интервалы уверенности

Поскольку стандартная ошибка выборки измеряет степень близости, с которой, по нашим ожиданиям, выборочное среднее соответствует среднему популяционному, то мы можем также рассмотреть обратное стандартная ошибка измеряет степень близосто, с которой, по нашим ожиданиям, популяционное среднее соответствует измеренному среднему выборочному. Другими словами, на основе стандартной ошибки мы можем вывести, что популяционное среднее находится в пределах некого ожидаемого диапазона выборочного среднего с некоторой степенью уверенности.

Понятия степень уверенности и ожидаемый диапазон, взятые вместе, дают определение термину интервал уверенности.

Примечание. В большинстве языков под термином confidence в контексте инференциальной статистики понимается именно уверенность в отличие от отечественной статистики, где принято говорить о доверительном интервале. На самом деле речь идет не о доверии (trust), а об уверенности исследователя в полученных результатах. Это яркий пример мягкой подмены понятия. Подобного рода ошибки вполне объяснимы - первые переводы появились еще в доколумбову доинтернетовскую эпоху, когда источников было мало, и приходилось домысливать в силу своего понимания. Сегодня же, когда существует масса профильных глоссариев, словарей и источников, такого рода искажения не оправданы. ИМХО.

При установлении интервалов уверенности обычной практикой является задание интервала размером 95% мы на 95% уверены, что популяционный параметр находится внутри интервала. Разумеется, еще остается 5%-я возможность, что он там не находится.

Какой бы ни была стандартная ошибка, 95% популяционного среднего значения будет находиться между -1.96 и 1.96 стандартных отклонений от выборочного среднего. И, следовательно, число 1.96 является критическим значением для 95%-ого интервала уверенности. Это критическое значение носит название z-значения.

^{Название z-значение вызвано тем, что нормальное распределение называется z-распределением. Впрочем, иногда z-значение так и называют гауссовым значением.}

Число 1.96 используется так широко, что его стоит запомнить. Впрочем, критическое значение мы можем вычислить сами, воспользовавшись функцией scipy stats.norm.ppf. Приведенная ниже функция confidence_interval ожидает значение для p между 0 и 1. Для нашего 95%-ого интервала уверенности оно будет равно 0.95. В целях вычисления положения каждого из двух хвостов нам нужно вычесть это число из единицы и разделить на 2 (2.5% для интервала уверенности шириной 95%):

def confidence_interval(p, xs): '''Интервал уверенности''' mu = xs.mean() se = standard_error(xs) z_crit = stats.norm.ppf(1 - (1-p) / 2)  return [mu - z_crit * se, mu + z_crit * se]def ex_2_9(): '''Вычислить интервал уверенности для данных за определенный день''' may_1 = '2015-05-01' df = with_parsed_date( load_data('dwell-times.tsv') )  filtered = df.set_index( ['date'] )[may_1] ci = confidence_interval(0.95, filtered['dwell-time']) print('Интервал уверенности: ', ci)

Интервал уверенности: [83.53415272762004, 97.753065317492741]

Полученный результат говорит о том, что можно на 95% быть уверенным в том, что популяционное среднее находится между 83.53 и 97.75 сек. И действительно, популяционное среднее, которое мы вычислили ранее, вполне укладывается в этот диапазон.

Сравнение выборок

После вирусной маркетинговой кампании веб-команда в AcmeContent извлекает для нас однодневную выборку времени пребывания посетителей на веб-сайте для проведения анализа. Они хотели бы узнать, не привлекла ли их недавняя кампания более активных посетителей веб-сайта. Интервалы уверенности предоставляют нам интуитивно понятный подход к сравнению двух выборок.

Точно так же, как мы делали ранее, мы загружаем значения времени пребывания, полученные в результате маркетинговой кампании, и их резюмируем:

def ex_2_10():    '''Сводные статистики данных, полученных        в результате вирусной кампании'''    ts = load_data('campaign-sample.tsv')['dwell-time']         print('n:                      ', ts.count())    print('Среднее:                ', ts.mean())    print('Медиана:                ', ts.median())    print('Стандартное отклонение: ', ts.std())    print('Стандартная ошибка:     ', standard_error(ts))    ex_2_10()

n:                       300Среднее:                 130.22Медиана:                 84.0Стандартное отклонение:  136.13370714388034Стандартная ошибка:      7.846572839994115

Среднее значение выглядит намного больше, чем то, которое мы видели ранее 130 сек. по сравнению с 90 сек. Вполне возможно, здесь имеется некое значимое расхождение, хотя стандартная ошибка более чем в 2 раза больше той, которая была в предыдущей однодневной выборке, в силу меньшего размера выборки и большего стандартного отклонения. Основываясь на этих данных, можно вычислить 95%-й интервал уверенности для популяционного среднего, воспользовавшись для этого той же самой функцией confidence_interval, что и прежде:

def ex_2_11(): '''Интервал уверенности для данных, полученных в результате вирусной кампании''' ts = load_data('campaign-sample.tsv')['dwell-time']  print('Интервал уверенности:', confidence_interval(0.95, ts))

Интервал уверенности: [114.84099983154137, 145.59900016845864]

95%-ый интервал уверенности для популяционного среднего лежит между 114.8 и 145.6 сек. Он вообще не пересекается с вычисленным нами ранее популяционным средним в размере 90 сек. Похоже, имеется какое-то крупное расхождение с опорной популяцией, которое едва бы произошло по причине одной лишь ошибки выборочного обследования. Наша задача теперь состоит в том, чтобы выяснить почему это происходит.

^{Ошибка выборочного обследования, также систематическая ошибка при взятии выборки, возникает, когда статистические характеристики популяции оцениваются исходя из подмножества этой популяции.}

Искаженность

Выборка должна быть репрезентативной, то есть представлять популяцию, из которой она взята. Другими словами, при взятии выборки необходимо избегать искажения, которое происходит в результате того, что отдельно взятые члены популяции систематически из нее исключаются (либо в нее включаются) по сравнению с другими.

Широко известным примером искажения при взятии выборки является опрос населения, проведенный в США еженедельным журналом Литературный Дайджест (Literary Digest) по поводу президентских выборов 1936 г. Это был один из самых больших и самых дорогостоящих когда-либо проводившихся опросов: тогда по почте было опрошено 2.4 млн. человек. Результаты были однозначными губернатор-республиканец от шт. Канзас Альфред Лэндон должен был победить Франклина Д. Рузвельта с 57% голосов. Как известно, в конечном счете на выборах победил Рузвельт с 62% голосов.

Первопричина допущенной журналом огромной ошибки выборочного обследования состояла в искажении при отборе. В своем стремлении собрать как можно больше адресов избирателей журнал Литературный Дайджест буквально выскреб все телефонные справочники, подписные перечни журнала и списки членов клубов. В эру, когда телефоны все еще во многом оставались предметом роскоши, такая процедура гарантировано имела избыточный вес в пользу избирателей, принадлежавших к верхнему и среднему классам, и не был представительным для электората в целом. Вторичной причиной искажения стала искаженность в результате неответов в опросе фактически согласились участвовать всего менее четверти тех, к кому обратились. В этом виде искаженности при отборе предпочтение отдается только тем респондентам, которые действительно желают принять участие в голосовании.

Распространенный способ избежать искаженности при отборе состоит в том, чтобы каким-либо образом рандомизировать выборку. Введение в процедуру случайности делает вряд ли возможным, что экспериментальные факторы окажут неправомерное влияние на качество выборки. Опрос населения еженедельником Литературный Дайджест был сосредоточен на получении максимально возможной выборки, однако неискаженная малая выборка была бы намного полезнее, чем плохо отобранная большая выборка.

Если мы откроем файл campaign_sample.tsv, то обнаружим, что наша выборка приходится исключительно на 6 июня 2015 года. Это был выходной день, и этот факт мы можем легко подтвердить при помощи функции pandas:

'''Проверка даты''' d = pd.to_datetime('2015 6 6') d.weekday() in [5,6]

True

Все наши сводные статистики до сих пор основывались на данных, которые мы отфильтровывали для получения только рабочих дней. Искаженность в нашей выборке вызвана именно этим фактом, и, если окажется, что поведение посетителей в выходные отличается от поведения в будние дни вполне возможный сценарий тогда мы скажем, что выборки представляют две разные популяции.

Визуализация разных популяций

Теперь снимем фильтр для рабочих дней и построим график среднесуточного времени пребывания для всех дней недели рабочих и выходных:

def ex_2_12():  '''Построить график времени ожидания  по всем дням, без фильтра''' df = load_data('dwell-times.tsv') means = mean_dwell_times_by_date(df)['dwell-time'] means.hist(bins=20) plt.xlabel('Ежедневное время ожидания неотфильтрованное, сек.') plt.ylabel('Частота') plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

Распределение уже не является нормальным. Оно фактически является бимодальным, поскольку у него два пика. Второй меньший пик соответствует вновь добавленным выходным дням, и он ниже потому что количество выходных дней гораздо меньше количества рабочих дней, а также потому что стандартная ошибка распределения больше.

^{Распределения более чем с одним пиком обычно называются мультимодальными. Они могут указывать на совмещение двух или более нормальных распределений, и, следовательно, возможно, на совмещение двух или более популяций. Классическим примером бимодальности является распределение показателей роста людей, поскольку модальный рост мужчин выше, чем у женщин.}

Данные выходных дней имеют другие характеристики в отличие от данных будних дней. Мы должны удостовериться, что мы сравниваем подобное с подобным. Отфильтруем наш первоначальный набор данных только по выходным дням:

def ex_2_13():    '''Сводные статистики данных,       отфильтрованных только по выходным дням'''    df = with_parsed_date( load_data('dwell-times.tsv') )    df.index = df['date']    df = df[df['date'].index.dayofweek > 4]   # суббота-воскресенье    weekend_times = df['dwell-time']      print('n:                      ', weekend_times.count())    print('Среднее:                ', weekend_times.mean())    print('Медиана:                ', weekend_times.median())    print('Стандартное отклонение: ', weekend_times.std())    print('Стандартная ошибка:     ', standard_error(weekend_times))        

n:                       5860Среднее:                 117.78686006825939Медиана:                 81.0Стандартное отклонение:  120.65234077179436Стандартная ошибка:      1.5759770362547678

Итоговое среднее значение в выходные дни (на основе 6-ти месячных данных) составляет 117.8 сек. и попадает в пределы 95%-ого интервала уверенности для маркетинговой выборки. Другими словами, хотя среднее значение времени пребывания в размере 130 сек. является высоким даже для выходных, расхождение не настолько большое, что его нельзя было бы приписать простой случайной изменчивости в выборке.

Мы только что применили подход к установлению подлинного расхождения в популяциях (между посетителями веб-сайта в выходные по сравнению с посетителями в будние дни), который при проведении проверки обычно не используется. Более традиционный подход начинается с выдвижения гипотетического предположения, после чего это предположение сверяется с данными. Для этих целей статистический метод анализа определяет строгий подход, который называется проверкой статистических гипотез.

Это и будет темой следующего поста, поста 3.

Примеры исходного кода для этого поста находятся в моемрепона Github. Все исходные данные взяты врепозиторииавтора книги.

Этот заключительный пост посвящен анализу дисперсии. Предыдущий пост см. здесь.

Анализ дисперсии

Анализ дисперсии (варианса), который в специальной литературе также обозначается как ANOVA от англ. ANalysis Of VAriance, это ряд статистических методов, используемых для измерения статистической значимости расхождений между группами. Он был разработан чрезвычайно одаренным статистиком Рональдом Фишером, который также популяризировал процедуру проверки статистической значимости в своих исследовательских работах по биологическому тестированию.

Примечание. В предыдущей и этой серии постов для термина variance использовался принятый у нас термин дисперсия и в скобках местами указывался термин варианс. Это не случайно. За рубежом существуют парные термины variance и covariance, и они по идее должны переводиться с одним корнем, например, как варианс и коварианс, однако на деле у нас парная связь разорвана, и они переводятся как совершенно разные дисперсия и ковариация. Но это еще не все. Dispersion (статистическая дисперсия) за рубежом является отдельным родовым понятием разбросанности, т.е. степени, с которой распределение растягивается или сжимается, а мерами статистической дисперсии являются варианс, стандартное отклонение и межквартильный размах. Dispersion, как родовое понятие разбросанности, и variance, как одна из ее мер, измеряющая расстояние от среднего значения - это два разных понятия. Далее в тексте для variance везде будет использоваться общепринятый термин дисперсия. Однако данное расхождение в терминологии следует учитывать.

Наши тесты на основе z-статистики и t-статистики были сосредоточены на выборочных средних значениях как первостепенном механизме проведения разграничения между двумя выборками. В каждом случае мы искали расхождение в средних значениях, деленных на уровень расхождения, который мы могли обоснованно ожидать, и количественно выражали стандартной ошибкой.

Среднее значение не является единственным выборочным индикатором, который может указывать на расхождение между выборками. На самом деле, в качестве индикатора статистического расхождения можно также использовать выборочную дисперсию.

В целях иллюстрации того, как это могло бы работать, рассмотрим приведенную выше диаграмму. Каждая из трех групп слева может представлять выборки времени пребывания на конкретном веб-сайте с его собственным средним значением и стандартным отклонением. Если время пребывания для всех трех групп объединить в одну, то дисперсия будет больше средней дисперсии для групп, взятых отдельно.

Статистическая значимость на основе анализа дисперсии вытекает из соотношения двух дисперсий дисперсии между исследуемыми группами и дисперсии внутри исследуемых групп. Если существует значимое межгрупповое расхождение, которое не отражено внутри групп, то эти группы помогают объяснить часть дисперсии между группами. И напротив, если внутригрупповая дисперсия идентична межгрупповой дисперсии, то группы статистически от друг друга неразличимы.

F-распределение

F-распределение параметризуется двумя степенями свободы степенями свободы размера выборки и числа групп.

Первая степень свободы это количество групп минус 1, и вторая степень свободы размер выборки минус число групп. Если k представляет число групп, и n объем выборки, то получаем:

Мы можем визуализировать разные F-распределения на графике при помощи функции библиотеки pandas plot:

def ex_2_Fisher(): '''Визуализация разных F-распределений на графике''' mu = 0 d1_values, d2_values = [4, 9, 49], [95, 90, 50] linestyles = ['-', '--', ':', '-.'] x = sp.linspace(0, 5, 101)[1:]  ax = None for (d1, d2, ls) in zip(d1_values, d2_values, linestyles): dist = stats.f(d1, d2, mu) df = pd.DataFrame( {0:x, 1:dist.pdf(x)} )  ax = df.plot(0, 1, ls=ls,  label=r'$d_1=%i,\ d_2=%i$' % (d1,d2), ax=ax) plt.xlabel('$x$\nF-статистика') plt.ylabel('Плотность вероятности \n$p(x|d_1, d_2)$') plt.show()

Кривые приведенного выше графика показывают разные F-распределения для выборки, состоящей из 100 точек, разбитых на 5, 10 и 50 групп.

F-статистика

Тестовая выборочная величина, которая представляет соотношение дисперсии внутри и между группами, называется F-статистикой. Чем ближе значение F-статистики к единице, тем более похожи обе дисперсии. F-статистика вычисляется очень просто следующим образом:

где S²_b это межгрупповая дисперсия, и S²_w внутригрупповая дисперсия.

По мере увеличения соотношения F межгрупповая дисперсия увеличивается по сравнению с внутригрупповой. Это означает, что это разбиение на группы хорошо справляется с задачей объяснения дисперсии, наблюдавшейся по всей выборке в целом. Там, где это соотношение превышает критический порог, мы можем сказать, что расхождение является статистически значимым.

^{F-тест всегда является односторонним, потому что любая дисперсия среди групп демонстрирует тенденцию увеличивать F. При этом F не может уменьшаться ниже нуля.}

Внутригрупповая дисперсия для F-теста вычисляется как среднеквадратичное отклонение от среднего значения. Мы вычисляем ее как сумму квадратов отклонений от среднего значения, деленную на первую степень свободы. Например, если имеется kгрупп, каждая со средним значением x_k, то мы можем вычислить внутригрупповую дисперсию следующим образом:

где SSW это внутригрупповая сумма квадратов, и x_jk это значение j-ого элемента в группе .

Приведенная выше формула для вычисления SSW имеет грозный вид, но на деле довольно легко имплементируется на Python, как сумма квадратичных отклонений от среднего значения ssdev, делающая вычисление внутригрупповой суммы квадратов тривиальным:

def ssdev( xs ): '''Сумма квадратов отклонений между  каждым элементом и средним по выборке''' mu = xs.mean()  square_deviation = lambda x : (x - mu) ** 2  return sum( map(square_deviation, xs) )

Межгрупповая дисперсия для F-теста имеет похожую формулу:

где SST это полная сумма квадратов отклонений, и SSW значение, которое мы только что вычислили. Полная сумма квадратов является суммой квадратичных расхождений от итогового среднего значения, которую можно вычислить следующим образом:

Отсюда, SST это попросту полная сумма квадратов без какого-либо разбиения на группы. На языке Python значения SST и SSW вычисляются элементарно, как будет показано ниже.

ssw = sum( groups.apply( lambda g: ssdev(g) ) ) # внутригрупповая сумма # квадратов отклонений sst = ssdev( df['dwell-time'] ) # полная сумма квадратов по всему наборуssb = sst  ssw # межгрупповая сумма квадратов отклонений

F-статистика вычисляется как отношение межгрупповой дисперсии к внутригрупповой. Объединив определенные ранее функции ssb и ssw и две степени свободы, мы можем вычислить F-статистика.

На языке Python F-статистика из групп и двух степеней свободы вычисляется следующим образом:

msb = ssb / df1 # усредненная межгрупповаяmsw = ssw / df2 # усредненная внутригрупповаяf_stat = msb / msw

Имея возможность вычислить F-статистику из групп, мы теперь готовы использовать его в соответствующем F-тесте.

F-тест

Как и в случае со всеми проверками статистических гипотез, которые мы рассмотрели в этой серии постов, как только имеется выборочная величина (статистика) и распределение, нам попросту нужно подобрать значение уровня значимости и посмотреть, не превысили ли наши данные критическое для теста значение.

Библиотека scipy предлагает функцию stats.f.sf, но она измеряет дисперсию между и внутри всего двух групп. В целях выполнения F-теста на наших 20 разных группах, нам придется имплементировать для нее нашу собственную функцию. К счастью, мы уже проделали всю тяжелую работу в предыдущих разделах, вычислив надлежащую F-статистику. Мы можем выполнить F-тест, отыскав F-статистику в F-распределении, параметризованном правильными степенями свободы. В следующем ниже примере мы напишем функцию f_test, которая все это использует для выполнения теста на произвольном числе групп:

def f_test(groups): m, n = len(groups), sum(groups.count()) df1, df2 = m - 1, n - m  ssw = sum( groups.apply(lambda g: ssdev(g)) )  sst = ssdev( df['dwell-time'] )  ssb = sst - ssw  msb = ssb / df1  msw = ssw / df2  f_stat = msb / msw return stats.f.sf(f_stat, df1, df2)    def ex_2_24(): '''Проверка вариантов дизайна веб-сайта на основе F-теста''' df = load_data('multiple-sites.tsv') groups = df.groupby('site')['dwell-time'] return f_test(groups)

0.014031745203658217

В последней строке приведенной выше функции мы преобразуем значение F-статистики в p-значение, пользуясь функцией scipy stats.f.sf, параметризованной правильными степенями свободы. P-значение является мерой всей модели, т.е. насколько хорошо разные веб-сайты объясняют дисперсию времени пребывания в целом. Нам остается только выбрать уровень значимости и выполнить проверку. Будем придерживаться 5%-ого уровня значимости.

Проверка возвращает p-значение, равное 0.014, т.е. значимый результат. Разные варианты веб-сайта действительно имеют разные дисперсии, которые нельзя просто объяснить одной лишь случайной ошибкой в выборке.

Для визуализации распределений всех вариантов дизайна веб-сайта на одном графике мы можем воспользоваться коробчатой диаграммой, разместив распределения для сопоставления рядом друг с другом:

def ex_2_25(): '''Визуализация распределений всех вариантов  дизайна веб-сайта на одной коробчатой диаграмме''' df = load_data('multiple-sites.tsv') df.boxplot(by='site', showmeans=True) plt.xlabel('Номер дизайна веб-сайта') plt.ylabel('Время пребывания, сек.') plt.title('') plt.suptitle('') plt.show()

В приведенном выше примере показана работа функции boxplot, которая вычисляет свертку на группах и сортирует группы по номеру варианта дизайна веб-сайта. Соответственно наш изначальный веб-сайт с номером 0, на графике расположен крайним слева.

Может создастся впечатление, что вариант дизайна веб-сайта с номером 10 имеет самое длительное время пребывания, поскольку его межквартильный размах простирается вверх выше других. Однако, если вы присмотритесь повнимательнее, то увидите, что его среднее значение меньше, чем у варианта дизайна с номером 6, имеющего среднее время пребывания более 144 сек.:

def ex_2_26(): '''T-проверка вариантов 0 и 10 дизайна веб-сайта''' df = load_data('multiple-sites.tsv') groups = df.groupby('site')['dwell-time'] site_0 = groups.get_group(0)  site_10 = groups.get_group(10) _, p_val = stats.ttest_ind(site_0, site_10, equal_var=False) return p_val

0.0068811940138903786

Подтвердив статистически значимый эффект при помощи F-теста, теперь мы вправе утверждать, что вариант дизайна веб-сайта с номером 6 статистически отличается от изначального значения:

def ex_2_27(): '''t-тест вариантов 0 и 6 дизайна веб-сайта''' df = load_data('multiple-sites.tsv') groups = df.groupby('site')['dwell-time'] site_0 = groups.get_group(0)  site_6 = groups.get_group(6) _, p_val = stats.ttest_ind(site_0, site_6, equal_var=False) return p_val

0.005534181712508717

Наконец, у нас есть подтверждающие данные, из которых вытекает, что веб-сайт с номером 6 является подлинным улучшением существующего веб-сайта. В результате нашего анализа исполнительный директор AcmeContent санкционирует запуск обновленного дизайна веб-сайта. Веб-команда - в восторге!

Размер эффекта

В этой серии постов мы сосредоточили наше внимание на широко используемых в статистической науке методах проверки статистической значимости, которые обеспечивают обнаружение статистического расхождения, которое не может быть легко объяснено случайной изменчивостью. Мы должны всегда помнить, что выявление значимого эффекта это не одно и то же, что и выявление большого эффекта. В случае очень больших выборок значимым будет считаться даже крошечное расхождение в выборочных средних. В целях более глубокого понимания того, является ли наше открытие значимым и важным, мы так же должны констатировать величину эффекта.

Интервальный индекс d Коэна

Индекс dКоэна это поправка, которую применяется для того, чтобы увидеть, не является ли наблюдавшееся расхождение не просто статистически значимым, но и действительно большим. Как и поправка Бонферрони, она проста:

Здесь S_ab это объединенное стандартное отклонение (не объединенная стандартная ошибка) выборок. Она вычисляется аналогично вычислению объединенной стандартной ошибки:

def pooled_standard_deviation(a, b): '''Объединенное стандартное отклонение  (не объединенная стандартная ошибка)''' return sp.sqrt( standard_deviation(a) ** 2 + standard_deviation(b) ** 2)

Так, для варианта под номером 6 дизайна нашего веб-сайта мы можем вычислить индекс dКоэна следующим образом:

def ex_2_28(): '''Вычисление интервального индекса d Коэна  для варианта дизайна веб-сайта под номером 6''' df = load_data('multiple-sites.tsv') groups = df.groupby('site')['dwell-time'] a = groups.get_group(0) b = groups.get_group(6) return (b.mean() - a.mean()) / pooled_standard_deviation(a, b)

0.38913648705499848

В отличие от p-значений, абсолютный порог для индекса d Коэна отсутствует. Считать ли эффект большим или нет частично зависит от контекста, однако этот индекс действительно предоставляет полезную, нормализованную меру величины эффекта. Значения выше 0.5, как правило, считаются большими, поэтому значение 0.38 это умеренный эффект. Он определенно говорит о значительном увеличении времени пребывания на нашем веб-сайте и что усилия, потраченные на обновление веб-сайта, определенно не были бесполезными.

Примеры исходного кода для этого поста находятся в моемрепона Github. Все исходные данные взяты врепозиторииавтора книги.

Резюме

В этой серии постов мы узнали о разнице между описательной и инференциальной статистикой, занимающейся методами статистического вывода. Мы еще раз убедились в важности нормального распределения и центральной предельной теоремы, и научились квантифицировать расхождение с популяциями, используя проверку статистических гипотез на основе z-теста, t-теста и F-теста.

Мы узнали о том, каким именно образом методология инференциальной статистики анализирует непосредственно сами выборки, чтобы выдвигать утверждения о всей популяции в целом, из которой они были отобраны. Мы познакомились с целым рядом методов интервалами уверенности, размножением выборок путем бутстрапирования и проверкой статистической значимости которые помогают разглядеть глубинные параметры популяции. Симулируя многократные испытания, мы также получили представление о трудности проверки статистической значимости в условиях многократных сравнений и увидели, как F-тест помогает решать эту задачу и устанавливать равновесие между ошибками 1-го и 2-го рода.

Мы также коснулись терминологических болевых точек и выяснили некоторые нюансы смыслового дрейфа в отечественной статистике.

В следующей серии постов, если читатели пожелают, мы применим полученные знания о дисперсии и F-тесте к одиночным выборкам. Мы представим метод регрессионного анализа и воспользуемся им для обнаружения корреляции между переменными в выборке из спортсменов-олимпийцев.

Предыдущий пост см. здесь.

Проверка статистических гипотез

Для статистиков и исследователей данных проверка статистической гипотезы представляет собой формальную процедуру. Стандартный подход к проверке статистической гипотезы подразумевает определение области исследования, принятие решения в отношении того, какие переменные необходимы для измерения предмета изучения, и затем выдвижение двух конкурирующих гипотез. Во избежание рассмотрения только тех данных, которые подтверждают наши субъективные оценки, исследователи четко констатируют свою гипотезу заранее. Затем, основываясь на данных, они применяют выборочные статистики с целью подтвердить либо отклонить эту гипотезу.

Проверка статистической гипотезы подразумевает использование тестовой статистики, т.е. выборочной величины, как функции от результатов наблюдений. Тестовая статистика (test statistic) - это вычисленная из выборочных данных величина, которая используется для оценивания прочности данных, подтверждающих нулевую статистическую гипотезу и служит для выявления меры расхождения между эмпирическими и гипотетическими значениями. Конкретные методы проверки называются тестами, например, z-тест, t-тест (соответственно z-тест Фишера, t-тест Студента) и т.д. в зависимости от применяемых в них тестовых статистик.

Примечание. В отечественной статистической науке используется туманный термин статистика критерия. Туманный потому здесь мы снова наблюдаем мягкую подмену: вместо теста возникает критерий. Если уж на то пошло, то критерий - это принцип или правило. Например, выполняя z-тест, t-тест и т.д., мы соответственно используем z-статистику, t-статистику и т.д. в правиле отклонения гипотезы. Это хорошо резюмируется следующей ниже таблицей:

Для того, чтобы помочь сохранить поток посетителей веб-сайта, дизайнеры приступают к работе над вариантом веб-сайта с использованием всех новейших методов по поддержанию внимания аудитории. Мы хотели бы удостовериться, что наши усилия не напрасны, и поэтому стараемся увеличить время пребывания посетителей на обновленном веб-сайте.

Отсюда главный вопрос нашего исследования состоит в том, приводит ли обновленный вид веб-сайта к увеличению времени пребывания на нем посетителей? Мы принимаем решение проверить его относительно среднего значения времени пребывания. Теперь, мы должны изложить две наши гипотезы. По традиции считается, что изучаемые данные не содержат того, что исследователь ищет. Таким образом, консервативное мнение заключается в том, что данные не покажут ничего необычного. Все это называется нулевой гипотезой и обычно обозначается как H₀.

^{При тестировании статистической гипотезы исходят из того, что нулевая гипотеза является истинной до тех пор, пока вес представленных данных, подтверждающих обратное, не сделает ее неправдоподобной. Этот подход к поиску доказательств в обратную сторону частично вытекает из простого психологического факта, что, когда люди пускаются на поиски чего-либо, они, как правило, это находят.}

Затем исследователь формулирует альтернативную гипотезу, обозначаемую как H₁. Она может попросту заключаться в том, что популяционное среднее отличается от базового уровня. Или же, что популяционное среднее больше или меньше базового уровня, либо больше или меньше на некоторую указанную величину. Мы хотели бы проверить, не увеличивает ли обновленный дизайн веб-сайта время пребывания, и поэтому нашей нулевой и альтернативной гипотезами будут следующие:

• H₀: Время пребывания для обновленного веб-сайта не отличается от времени пребывания для существующего веб-сайта

• H₁: Время пребывания для обновленного веб-сайта больше по сравнению с временем пребывания для существующего веб-сайта

Наше консервативное допущение состоит в том, что обновленный веб-сайт никак не влияет на время пребывания посетителей на веб-сайте. Нулевая гипотеза не обязательно должна быть нулевой (т.е. эффект отсутствует), но в данном случае, у нас нет никакого разумного оправдания, чтобы считать иначе. Если выборочные данные не поддержат нулевую гипотезу (т.е. если данные расходятся с ее допущением на слишком большую величину, чтобы носить случайный характер), то мы отклоним нулевую гипотезу и предложим альтернативную в качестве наилучшего альтернативного объяснения.

Указав нулевую и альтернативную гипотезы, мы должны установить уровень значимости, на котором мы ищем эффект.

Статистическая значимость

Проверка статистической значимости изначально разрабатывалась независимо от проверки статистических гипотез, однако сегодня оба подхода очень часто используются во взаимодействии друг с другом. Задача проверки статистической значимости состоит в том, чтобы установить порог, за пределами которого мы решаем, что наблюдаемые данные больше не поддерживают нулевую гипотезу.

Следовательно, существует два риска:

• Мы можем принять расхождение как значимое, когда на самом деле оно возникло случайным образом

• Мы можем приписать расхождение случайности, когда на самом деле оно показывает истинное расхождение с популяцией

Эти две возможности обозначаются соответственно, как ошибки 1-го и 2-го рода:

Чем больше мы уменьшаем риск совершения ошибок 1-го рода, тем больше мы увеличиваем риск совершения ошибок 2-го рода. Другими словами, с чем большей уверенностью мы хотим не заявлять о наличии расхождения, когда его нет, тем большее расхождение между выборками нам потребуется, чтобы заявить о статистической значимости. Эта ситуация увеличивает вероятность того, что мы проигнорируем подлинное расхождение, когда мы с ним столкнемся.

В статистической науке обычно используются два порога значимости. Это уровни в 5% и 1%. Расхождение в 5% обычно называют значимым, а расхождение в 1% крайне значимым. В формулах этот порог часто обозначается греческой буквой (альфа) и называется уровнем значимости. Поскольку, отсутствие эффекта по результатам эксперимента может рассматриваться как неуспех (эксперимента либо обновленного веб-сайта, как в нашем случае), то может возникнуть желание корректировать уровень значимости до тех пор, пока эффект не будет найден. По этой причине классический подход к проверке статистической значимости требует, чтобы мы устанавливали уровень значимости до того, как обратимся к нашим данным. Часто выбирается уровень в 5%, и поэтому мы на нем и остановимся.

Проверка обновленного дизайна веб-сайта

Веб-команда в AcmeContent была поглощена работой, конструируя обновленный веб-сайт, который будет стимулировать посетителей оставаться на нем в течение более длительного времени. Она употребила все новейшие методы и, в результате мы вполне уверены, что веб-сайт покажет заметное улучшение показателя времени пребывания.

Вместо того, чтобы запустить его для всех пользователей сразу, в AcmeContent хотели бы сначала проверить веб-сайт на небольшой выборке посетителей. Мы познакомили веб-команду с понятием искаженности выборки, и в результате там решили в течение одного дня перенаправлять случайные 5% трафика на обновленный веб-сайт. Результат с дневным трафиком был нам предоставлен одним текстовым файлом. Каждая строка показывает время пребывания посетителей. При этом, если посетитель пользовался исходным дизайном, ему присваивалось значение "0", и если он пользовался обновленным (и надеемся, улучшенным) дизайном, то ему присваивалось значение "1".

Выполнение z-теста

Ранее при тестировании с интервалами уверенности мы располагали лишь одним популяционным средним, с которым и выполнялось сравнение.

При тестировании нулевой гипотезы с помощью z-теста мы имеем возможность сравнивать две выборки. Посетители, которые видели обновленный веб-сайт, отбирались случайно, и данные для обеих групп были собраны в тот же день, чтобы исключить другие факторы с временной зависимостью.

Поскольку в нашем распоряжении имеется две выборки, то и стандартных ошибок у нас тоже две. Z-тест выполняется относительно объединенной стандартной ошибки, т.е. квадратного корня суммы дисперсий (вариансов), деленных на размеры выборок. Она будет такой же, что и результат, который мы получим, если взять стандартную ошибку обеих выборок вместе:

Здесь ²_a это дисперсия выборки a, ²_b дисперсия выборки bи соответственно n_a и n_b размеры выборок a и b. На Python объединенная стандартная ошибка вычисляется следующим образом:

def pooled_standard_error(a, b, unbias=False): '''Объединенная стандартная ошибка''' std1 = a.std(ddof=0) if unbias==False else a.std()  std2 = b.std(ddof=0) if unbias==False else b.std() x = std1 ** 2 / a.count() y = std2 ** 2 / b.count() return sp.sqrt(x + y)

С целью выявления того, является ли видимое нами расхождение неожиданно большим, можно взять наблюдавшиеся расхождения между средними значениями на объединенной стандартной ошибке. Эту статистическую величину принято обозначать переменной z:

Используя функции pooled_standard_error, которая вычисляет объединенную стандартную ошибку, z-статистику можно получить следующим образом:

def z_stat(a, b, unbias=False): return (a.mean() - b.mean()) / pooled_standard_error(a, b, unbias)

Соотношение z объясняет, насколько средние значения отличаются относительно величины, которую мы ожидаем при заданной стандартной ошибке. Следовательно, z-статистика сообщает нам о том, на какое количество стандартных ошибок расходятся средние значения. Поскольку стандартная ошибка имеет нормальное распределение вероятностей, мы можем связать это расхождение с вероятностью, отыскав z-статистику в нормальной ИФР:

def z_test(a, b):  return stats.norm.cdf([ z_stat(a, b) ])

В следующем ниже примере z-тест используется для сравнения результативность двух веб-сайтов. Это делается путем группировки строк по номеру веб-сайта, в результате чего возвращается коллекция, в которой конкретному веб-сайту соответствует набор строк. Мы вызываем groupby('site')['dwell-time'] для конвертирования набора строк в набор значений времени пребывания. Затем вызываем функцию get_group с номером группы, соответствующей номеру веб-сайта:

def ex_2_14():    '''Сравнение результативности двух вариантов       дизайна веб-сайта на основе z-теста'''    groups = load_data('new-site.tsv').groupby('site')['dwell-time']    a = groups.get_group(0)    b = groups.get_group(1)         print('a n:         ', a.count())    print('b n:         ', b.count())    print('z-статистика:', z_stat(a, b))    print('p-значение:  ', z_test(a, b))

a n:          284b n:          16z-статистика: -1.6467438180091214p-значение:   [0.04980536]

Установление уровня значимости в размере 5% во многом аналогично установлению интервала уверенности шириной 95%. В сущности, мы надеемся убедиться, что наблюдавшееся расхождение попадает за пределы 95%-го интервала уверенности. Если это так, то мы можем утверждать, что нашли результат с 5%-ым уровнем значимости.

^{P-значение это вероятность совершения ошибки 1-го рода в результате неправильного отклонения нулевой гипотезы, которая в действительности является истинной. Чем меньше p-значение, тем больше определенность в том, что нулевая гипотеза является ложной, и что мы нашли подлинный эффект.}

Этот пример возвращает значение 0.0498, или 4.98%. Поскольку оно немногим меньше нашего 5% порога значимости, мы можем утверждать, что нашли нечто значимое.

Приведем еще раз нулевую и альтернативную гипотезы:

• H₀: Время пребывания на обновленном веб-сайте не отличается от времени пребывания на существующем веб-сайте

• H₁: Время пребывания на обновленном веб-сайте превышает время пребывания на существующем веб-сайте.

Наша альтернативная гипотеза состоит в том, что время пребывания на обновленном веб-сайте больше.

Мы готовы заявить о статистической значимости, и что время пребывания на обновленном веб-сайте больше по сравнению с существующим веб-сайтом, но тут есть одна трудность чем меньше размер выборки, тем больше неопределенность в том, что выборочное стандартное отклонение совпадет с популяционным. Как показано в результатах предыдущего примера, наша выборка из обновленного веб-сайта содержит всего 16 посетителей. Столь малые выборки делают невалидным допущение о том, что стандартная ошибка нормально распределена.

К счастью, существует тест и связанное с ним распределение, которое моделирует увеличенную неопределенность стандартных ошибок для выборок меньших размеров.

t-распределение Студента

Популяризатором t-распределения был химик, работавший на пивоварню Гиннес в Ирландии, Уилльям Госсетт, который включил его в свой анализ темного пива Стаут.

^{В 1908 Уильям Госсет опубликовал статью об этой проверке в журнале Биометрика, но при этом по распоряжению своего работодателя, который рассматривал использованную Госсеттом статистику как коммерческую тайну, был вынужден использовать псевдоним. Госсет выбрал псевдоним Студент.}

В то время как нормальное распределение полностью описывается двумя параметрами средним значением и стандартным отклонением, t-распределение описывается лишь одним параметром, так называемыми степенями свободы. Чем больше степеней свободы, тем больше t-распределение похоже на нормальное распределение с нулевым средним и стандартным отклонением, равным 1. По мере уменьшения степеней свободы, это распределение становится более широким с более толстыми чем у нормального распределения, хвостами.

Приведенный выше рисунок показывает, как t-распределение изменяется относительно нормального распределения при наличии разных степеней свободы. Более толстые хвосты для выборок меньших размеров соответствуют увеличенной возможности наблюдать более крупные отклонения от среднего значения.

Степени свободы

Степени свободы, часто обозначаемые сокращенно df от англ. degrees of freedom, тесно связаны с размером выборки. Это полезная статистика и интуитивно понятное свойство числового ряда, которое можно легко продемонстрировать на примере.

Если бы вам сказали, что среднее, состоящее из двух значений, равно 10 и что одно из значений равно 8, то Вам бы не потребовалась никакая дополнительная информация для того, чтобы суметь заключить, что другое значение равно 12. Другими словами, для размера выборки, равного двум, и заданного среднего значения одно из значений ограничивается, если другое известно.

Если напротив вам говорят, что среднее, состоящее из трех значений, равно 10, и первое значение тоже равно 10, то Вы были бы не в состоянии вывести оставшиеся два значения. Поскольку число множеств из трех чисел, начинающихся с 10, и чье среднее равно 10, является бесконечным, то прежде чем вы сможете вывести значение третьего, второе тоже должно быть указано.

Для любого множества из трех чисел ограничение простое: вы можете свободно выбрать первые два числа, но заключительное число ограничено. Степени свободы могут таким образом быть обобщены следующим образом: количество степеней свободы любой отдельной выборки на единицу меньше размера выборки.

При сопоставлении двух выборок степени свободы на две единицы меньше суммы размеров этих выборок, что равно сумме их индивидуальных степеней свободы.

t-статистика

При использовании t-распределения мы обращаемся к t-статистике. Как и z-статистика, эта величина количественно выражает степень маловероятности отдельно взятого наблюдавшегося отклонения. Для двухвыборочного t-теста соответствующая t-статистика вычисляется следующим образом:

Здесь S_ab это объединенная стандартная ошибка. Объединенная стандартная ошибка вычисляется таким же образом, как и раньше:

Однако это уравнение допускает наличие информации о популяционных параметрах _aи _b, которые можно аппроксимировать только на основе крупных выборок. t-тест предназначен для малых выборок и не требует от нас принимать допущения о поплуляционной дисперсии (вариансе).

Как следствие, объединенная стандартная ошибка для t-теста записывается как квадратный корень суммы стандартных ошибок:

На практике оба приведенных выше уравнения для объединенной стандартной ошибки дают идентичные результаты при заданных одинаковых входных последовательностях. Разница в математической записи всего лишь служит для иллюстрации того, что в условиях t-теста мы на входе зависим только от выборочных статистик. Объединенная стандартная ошибка может быть вычислена следующим образом:

def pooled_standard_error_t(a, b):  '''Объединенная стандартная ошибка для t-теста''' return sp.sqrt(standard_error(a) ** 2 +  standard_error(b) ** 2)

Хотя в математическом плане t-статистика и z-статистика представлены по-разному, на практике процедура вычисления обоих идентичная:

t_stat = z_statdef ex_2_15():    '''Вычисление t-статистики        двух вариантов дизайна веб-сайта'''    groups = load_data('new-site.tsv').groupby('site')['dwell-time']    a = groups.get_group(0)    b = groups.get_group(1)        return t_stat(a, b)

-1.6467438180091214

Различие между двумя выборочными показателями является не алгоритмическим, а концептуальным z-статистика применима только тогда, когда выборки подчинены нормальному распределению.

t-тест

Разница в характере работы t-теста вытекает из распределения вероятностей, из которого вычисляется наше p-значение. Вычислив t-статистику, мы должны отыскать ее значение в t-распределении, параметризованном степенями свободы наших данных:

def t_test(a, b): df = len(a) + len(b) - 2 return stats.t.sf([ abs(t_stat(a, b)) ], df)

Значение степени свободы обеих выборок на две единицы меньше их размеров, и для наших выборок составляет 298.

Напомним, что мы выполняем проверку статистической гипотезы. Поэтому выдвинем нашу нулевую и альтернативную гипотезы:

• H₀: Эта выборка взята из популяции с предоставленным средним значением

• H₁: Эта выборка взята из популяции со средним значением большего размера

Выполним следующий ниже пример:

def ex_2_16(): '''Сравнение результативности двух вариантов  дизайна веб-сайта на основе t-теста''' groups = load_data('new-site.tsv').groupby('site')['dwell-time'] a = groups.get_group(0) b = groups.get_group(1)  return t_test(a, b)

array([ 0.05033241])

Этот пример вернет p-значение, составляющее более 0.05. Поскольку оно больше , равного 5%, который мы установили для проверки нулевой гипотезы, то мы не можем ее отклонить. Наша проверка с использованием t-теста значимого расхождения между средними значениями не обнаружила. Следовательно, наш едва значимый результат z-теста отчасти объясняется наличием слишком малой выборки.

Двухсторонние тесты

В нашей альтернативной гипотезе было принято неявное допущение, что обновленный веб-сайт будет работать лучше существующего. В процедуре проверки нулевой статистической гипотезы предпринимаются особые усилия для обеспечения того, чтобы при поиске статистической значимости мы не делали никаких скрытых допущений.

Проверки, при выполнении которых мы ищем только значимое количественное увеличение или уменьшение, называются односторонними и обычно не приветствуются, кроме случая, когда изменение в противоположном направлении было бы невозможным. Название термина односторонний обусловлено тем, что односторонняя проверка размещает всю в одном хвосте распределения. Не делая проверок в другом направлении, проверка имеет больше мощности отклонить нулевую гипотезу в отдельно взятом направлении и, в сущности, понижает порог, по которому мы судим о результате как значимом.

^{Статистическая мощность это вероятность правильного принятия альтернативной гипотезы. Она может рассматриваться как способность проверки обнаруживать эффект там, где имеется искомый эффект.}

Хотя более высокая статистическая мощность выглядит желательной, она получается за счет наличия большей вероятности совершить ошибку 1-го рода. Правильнее было бы допустить возможность того, что обновленный веб-сайт может в действительности оказаться хуже существующего. Этот подход распределяет нашу одинаково по обоим хвостам распределения и обеспечивает значимый результат, не искаженный под воздействием априорного допущения об улучшении работы обновленного веб-сайта.

В действительности в модуле stats библиотеки scipy уже предусмотрены функции для выполнения двухвыборочных t-проверок. Это функция stats.ttest_ind. В качестве первого аргумента мы предоставляем выборку данных и в качестве второго - выборку для сопоставления. Если именованный аргумент equal_var равен True, то выполняется стандартная независимая проверка двух выборок, которая предполагает равные популяционные дисперсии, в противном случае выполняется проверка Уэлша (обратите внимание на служебную функцию t_test_verbose, (которую можно найти среди примеров исходного кода в репо):

def ex_2_17(): '''Двухсторонний t-тест''' groups = load_data('new-site.tsv').groupby('site')['dwell-time'] a = groups.get_group(0) b = groups.get_group(1)  return t_test_verbose(a, sample2=b, fn=stats.ttest_ind) #t-тест Уэлша

{'p-значение': 0.12756432502462475, 'степени свободы     ': 17.761382349686098, 'интервал уверенности': (76.00263198799597, 99.89877646270826), 'n1          ': 284, 'n2          ': 16, 'среднее x   ': 87.95070422535211, 'среднее y   ': 122.0, 'дисперсия x ': 10463.941024237296, 'дисперсия y ': 6669.866666666667, 't-статистика': -1.5985205593851322}

По результатам t-теста служебная функция t_test_verbose возвращает много информации и в том числе p-значение. P-значение примерно в 2 раза больше того, которое мы вычислили для односторонней проверки. На деле, единственная причина, почему оно не совсем в два раза больше, состоит в том, что в модуле stats имплементирован легкий вариант t-теста, именуемый t-тестом Уэлша, который немного более робастен, когда две выборки имеют разные стандартные отклонения. Поскольку мы знаем, что для экспоненциальных распределений среднее значение и дисперсия тесно связаны, то этот тест немного более строг в применении и даже возвращает более низкую значимость.

Одновыборочный t-тест

Независимые выборки в рамках t-тестов являются наиболее распространенным видом статистического анализа, который обеспечивает очень гибкий и обобщенный способ установления, что две выборки представляют одинаковую либо разную популяцию. Однако в случаях, когда популяционное среднее уже известно, существует еще более простая проверка, представленная функцией библиотеки sciзy stats.ttest_1samp.

Мы передаем выборку и популяционное среднее относительно которого выполняется проверка. Так, если мы просто хотим узнать, не отличается ли обновленный веб-сайт значимо от существующего популяционного среднего времени пребывания, равного 90 сек., то подобную проверку можно выполнить следующим образом:

def ex_2_18(): groups = load_data('new-site.tsv').groupby('site')['dwell-time'] b = groups.get_group(1)  return t_test_verbose(b, mean=90, fn=stats.ttest_1samp) 

{'p-значение          ': 0.13789520958229415, 'степени свободы df  ': 15.0, 'интервал уверенности': (78.4815276659039, 165.5184723340961), 'n1                  ': 16, 'среднее x           ': 122.0, 'дисперсия x         ': 6669.866666666667, 't-статистика        ': 1.5672973291495713}

Служебная функция t_test_verbose не только возвращает p-значение для выполненной проверки, но и интервал уверенности для популяционного среднего. Интервал имеет широкий диапазон между 78.5 и 165.5 сек., и, разумеется, перекрывается 90 сек. нашего теста. Как раз он и объясняет, почему мы не смогли отклонить нулевую гипотезу.

Многократные выборки

В целях развития интуитивного понимания относительно того, каким образом t-тест способен подтвердить и вычислить эти статистики из столь малых данных, мы можем применить подход, который связан с многократными выборками, от англ. resampling. Извлечение многократных выборок основывается на допущении о том, что каждая выборка является лишь одной из бесконечного числа возможных выборок из популяции. Мы можем лучше понять природу того, какими могли бы быть эти другие выборки, и, следовательно, добиться лучшего понимания опорной популяции, путем извлечения большого числа новых выборок из нашей существующей выборки.

На самом деле существует несколько методов взятия многократных выборок, и мы обсудим один из самых простых бутстрапирование. При бустрапировании мы генерируем новую выборку, неоднократно извлекая из исходной выборки случайное значение с возвратом до тех пор, пока не сгенерируем выборку, имеющую тот же размер, что и оригинал. Поскольку выбранные значения возвращаются назад после каждого случайного отбора, то в новой выборке то же самое исходное значение может появляться многократно. Это как если бы мы неоднократно вынимали случайную карту из колоды игральных карт и каждый раз возвращали вынутую карту назад в колоду. В результате время от времени мы будем иметь карту, которую мы уже вынимали.

^{Бутстраповская выборка, или бутстрап, синтетический набор данных, полученный в результате генерирования повторных выборок (с возвратом) из исследуемой выборки, используемой в качестве суррогатной популяции, в целях аппроксимации выборочного распределения статистики (такой как, среднее, медиана и др.).}

В библиотеке pandas при помощи функции sample можно легко извлекать бутстраповские выборки и генерировать большое число многократных выборок. Эта функция принимает ряд опциональных аргументов, в т.ч. n (число элементов, которые нужно вернуть из числового ряда), axis (ось, из которой извлекать выборку) и replace (выборка с возвратом или без), по умолчанию равный False. После этой функции можно задать метод агрегирования, вычисляющий сводную статистику в отношении бутстраповских выборок:

def ex_2_19(): '''Построение графика синтетических времен пребывания  путем извлечения бутстраповских выборок''' groups = load_data('new-site.tsv').groupby('site')['dwell-time'] b = groups.get_group(1)  xs = [b.sample(len(b), replace=True).mean() for _ in range(1000)]  pd.Series(xs).hist(bins=20) plt.xlabel('Бутстрапированные средние значения времени пребывания, сек.') plt.ylabel('Частота')  plt.show()

Приведенный выше пример наглядно показывает результаты на гистограмме:

Гистограмма демонстрирует то, как средние значения изменялись вместе с многократными выборками, взятыми из времени пребывания на обновленном веб-сайте. Хотя на входе имелась лишь одна выборка, состоящая из 16 посетителей, бутстрапированные выборки очень четко просимулировали стандартную ошибку изначальной выборки и позволили визуализировать интервал уверенности (между 78 и 165 сек.), вычисленный ранее в результате одновыборочного t-теста.

Благодаря бутстрапированию мы просимулировали взятие многократных выборок, при том, что у нас на входе имелась всего одна выборка. Этот метод обычно применяется для оценивания параметров, которые мы не способны или не знаем, как вычислить аналитически.

Проверка многочисленных вариантов дизайна

Было разочарованием обнаружить отсутствие статистической значимости на фоне увеличенного времени пребывания пользователей на обновленном веб-сайте. Хотя хорошо, что мы обнаружили это на малой выборке пользователей, прежде чем выкладывать его на всеобщее обозрение.

Не позволяя себя обескуражить, веб-команда AcmeContent берется за сверхурочную работу и создает комплект альтернативных вариантов дизайна веб-сайта. Беря лучшие элементы из других проектов, они разрабатывают 19 вариантов для проверки. Вместе с нашим изначальным веб-сайтом, который будет действовать в качестве контрольного, всего имеется 20 разных вариантов дизайна веб-сайта, куда посетители будут перенаправляться.

Вычисление выборочных средних

Веб-команда разворачивает 19 вариантов дизайна обновленного веб-сайта наряду с изначальным. Как отмечалось ранее, каждый вариант дизайна получает случайные 5% посетителей, и при этом наше испытание проводится в течение 24 часов.

На следующий день мы получаем файл, показывающий значения времени пребывания посетителей на каждом варианте веб-сайта. Все они были промаркированы числами, при этом число 0 соответствовало веб-сайту с исходным дизайном, а числа от 1 до 19 представляли другие варианты дизайна:

def ex_2_20(): df = load_data('multiple-sites.tsv') return df.groupby('site').aggregate(sp.mean)

Этот пример сгенерирует следующую ниже таблицу:

Мы хотели бы проверить каждый вариант дизайна веб-сайта, чтобы увидеть, не генерирует ли какой-либо из них статистически значимый результат. Для этого можно сравнить варианты дизайна веб-сайта друг с другом следующим образом, причем нам потребуется вспомогательный модуль Python itertools, который содержит набор функций, создающих итераторы для эффективной циклической обработки:

import itertoolsdef ex_2_21(): '''Проверка вариантов дизайна веб-сайта на основе t-теста по принципу "каждый с каждым"''' groups = load_data('multiple-sites.tsv').groupby('site') alpha = 0.05 pairs = [list(x) # найти сочетания из n по k for x in itertools.combinations(range(len(groups)), 2)]  for pair in pairs: gr, gr2 = groups.get_group( pair[0] ), groups.get_group( pair[1] ) site_a, site_b = pair[0], pair[1] a, b = gr['dwell-time'], gr2['dwell-time']  p_val = stats.ttest_ind(a, b, equal_var = False).pvalue  if p_val < alpha:  print('Варианты веб-сайта %i и %i значимо различаются: %f'  % (site_a, site_b, p_val))

Однако это было бы неправильно. Мы скорее всего увидим статистическое расхождение между вариантами дизайна, показавшими себя в особенности хорошо по сравнению с вариантами, показавшими себя в особенности плохо, даже если эти расхождения носили случайный характер. Если вы выполните приведенный выше пример, то увидите, что многие варианты дизайна веб-сайта статистически друг от друга отличаются.

С другой стороны, мы можем сравнить каждый вариант дизайна веб-сайта с нашим текущим изначальным значением средним значением времени пребывания, равным 90 сек., измеренным на данный момент для существующего веб-сайта:

def ex_2_22(): groups = load_data('multiple-sites.tsv').groupby('site') alpha = 0.05  baseline = groups.get_group(0)['dwell-time'] for site_a in range(1, len(groups)): a = groups.get_group( site_a )['dwell-time'] p_val = stats.ttest_ind(a, baseline, equal_var = False).pvalue  if p_val < alpha:  print('Вариант %i веб-сайта значимо отличается: %f'  % (site_a, p_val))

В результате этой проверки будут идентифицированы два варианта дизайна веб-сайта, которые существенно отличаются:

Вариант 6 веб-сайта значимо отличается: 0.005534Вариант 10 веб-сайта 10 значимо отличается: 0.006881

Малые p-значения (меньше 1%) указывают на то, что существует статистически очень значимые расхождения. Этот результат представляется весьма многообещающим, однако тут есть одна проблема. Мы выполнили t-тест по 20 выборкам данных с уровнем значимости , равным 0.05. Уровень значимости определяется, как вероятность неправильного отказа от нулевой гипотезы. На самом деле после 20-кратного выполнения t-теста становится вероятным, что мы неправильно отклоним нулевую гипотезу по крайней мере для одного варианта веб-сайта из 20.

Сравнивая таким одновременным образом многочисленные страницы, мы делаем результаты t-теста невалидными. Существует целый ряд альтернативных технических приемов решения проблемы выполнения многократных сравнений в статистических тестах. Эти методы будут рассмотрены в следующем разделе.

Поправка Бонферрони

Для проведения многократных проверок используется подход, который объясняет увеличенную вероятность обнаружить значимый эффект в силу многократных испытаний. Поправка Бонферрони это очень простая корректировка, которая обеспечивает, чтобы мы вряд ли совершили ошибки 1-го рода. Она выполняется путем настройки значения уровня значимости для тестов.

Настройка очень простая поправка Бонферрони попросту делит требуемое значение на число тестов. Например, если для теста имелось kвариантов дизайна веб-сайта, и эксперимента равно 0.05, то поправка Бонферрони выражается следующим образом:

Она представляет собой безопасный способ смягчить увеличение вероятности совершения ошибки 1-го рода при многократной проверке. Следующий пример идентичен примеру ex-2-22, за исключением того, что значение разделено на число групп:

def ex_2_23(): '''Проверка вариантов дизайна веб-сайта на основе t-теста против исходного (0) с поправкой Бонферрони''' groups = load_data('multiple-sites.tsv').groupby('site') alpha = 0.05 / len(groups) baseline = groups.get_group(0)['dwell-time'] for site_a in range(1, len(groups)): a = groups.get_group(site_a)['dwell-time'] p_val = stats.ttest_ind(a, baseline, equal_var = False).pvalue  if p_val < alpha:  print('Вариант %i веб-сайта значимо отличается от исходного: %f'  % (site_a, p_val))

Если вы выполните приведенный выше пример, то увидите, что при использовании поправки Бонферрони ни один из веб-сайтов больше не считается статистически значимым.

Метод проверки статистической значимости связан с поддержанием равновесия чем меньше шансы совершения ошибки 1-го рода, тем больше риск совершения ошибки 2-го рода. Поправка Бонферрони очень консервативна, и весьма возможно, что из-за излишней осторожности мы пропускаем подлинное расхождение.

Примеры исходного кода для этого поста находятся в моемрепона Github. Все исходные данные взяты врепозиторииавтора книги.

В заключительном посте, посте 4, этой серии постов мы проведем исследование альтернативного подхода к проверке статистической значимости, который позволяет устанавливать равновесие между совершением ошибок 1-го и 2-го рода, давая нам возможность проверить все 20 вариантов веб-сайта одновременно.

Чем больше я узнаю людей, тем больше мне нравится моя собака.

Марк Твен

В предыдущих сериях постов из ремикса книги Генри Гарнера Clojure для исследования данных (Clojure for Data Science) на языке Python мы рассмотрели методы описания выборок с точки зрения сводных статистик и методов статистического вывода из них параметров популяции. Такой анализ сообщает нам нечто о популяции в целом и о выборке в частности, но он не позволяет нам делать очень точные утверждения об их отдельных элементах. Это связано с тем, что в результате сведения данных всего к двум статистикам - среднему значению и стандартному отклонению - теряется огромный объем информации.

Нам часто требуется пойти дальше и установить связь между двумя или несколькими переменными либо предсказать одну переменную при наличии другой. И это подводит нас к теме данной серии из 5 постов - исследованию корреляции и регрессии. Корреляция имеет дело с силой и направленностью связи между двумя или более переменными. Регрессия определяет природу этой связи и позволяет делать предсказания на ее основе.

В этой серии постов будет рассмотрена линейная регрессия. При наличии выборки данных наша модель усвоит линейное уравнение, позволяющее ей делать предсказания о новых, не встречавшихся ранее данных. Для этого мы снова обратимся к библиотеке pandas и изучим связь между ростом и весом спортсменов-олимпийцев. Мы введем понятие матриц и покажем способы управления ими с использованием библиотеки pandas.

О данных

В этой серии постов используются данные, любезно предоставленные компанией Guardian News and Media Ltd., о спортсменах, принимавших участие в Олимпийских Играх 2012 г. в Лондоне. Эти данные изначально были взяты из блога газеты Гардиан.

Обследование данных

Когда вы сталкиваетесь с новым набором данных, первая задача состоит в том, чтобы его обследовать с целью понять, что именно он содержит.

Файл all-london-2012-athletes.tsv достаточно небольшой. Мы можем обследовать данные при помощи pandas, как мы делали в первой серии постов Python, исследование данных и выборы, воспользовавшись функцией read_csv:

def load_data(): return pd.read_csv('data/ch03/all-london-2012-athletes-ru.tsv', '\t') def ex_3_1(): '''Загрузка данных об участниках  олимпийских игр в Лондоне 2012 г.''' return load_data()

Если выполнить этот пример в консоли интерпретатора Python либо в блокноте Jupyter, то вы должны увидеть следующий ниже результат:

Столбцы данных (нам повезло, что они ясно озаглавлены) содержат следующую информацию:

• ФИО атлета

• страна, за которую он выступает

• возраст, лет

• рост, см.

• вес, кг.

• пол "М" или "Ж"

• дата рождения в виде строки

• место рождения в виде строки (со страной)

• число выигранных золотых медалей

• число выигранных серебряных медалей

• число выигранных бронзовых медалей

• всего выигранных золотых, серебряных и бронзовых медалей

• вид спорта, в котором он соревновался

• состязание в виде списка, разделенного запятыми

Даже с учетом того, что данные четко озаглавлены, очевидно присутствие пустых мест в столбцах с ростом, весом и местом рождения. При наличии таких данных следует проявлять осторожность, чтобы они не сбили с толку.

Визуализация данных

В первую очередь мы рассмотрим разброс роста спортсменов на Олимпийских играх 2012 г. в Лондоне. Изобразим эти значения роста в виде гистограммы, чтобы увидеть характер распределения данных, не забыв сначала отфильтровать пропущенные значения:

def ex_3_2(): '''Визуализация разброса значений  роста спортсменов на гистограмме''' df = load_data() df['Рост, см'].hist(bins=20) plt.xlabel('Рост, см.') plt.ylabel('Частота') plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

Как мы и ожидали, данные приближенно нормально распределены. Средний рост спортсменов составляет примерно 177 см. Теперь посмотрим на распределение веса олимпийских спортсменов:

def ex_3_3(): '''Визуализация разброса значений веса спортсменов''' df = load_data() df['Вес'].hist(bins=20) plt.xlabel('Вес') plt.ylabel('Частота') plt.show()

Приведенный выше пример сгенерирует следующую ниже гистограмму:

Данные показывают четко выраженную асимметрию. Хвост с правой стороны намного длиннее, чем с левой, и поэтому мы говорим, что асимметрия - положительная. Мы можем оценить асимметрию данных количественно при помощи функции библиотеки pandas skew:

def ex_3_4(): '''Вычисление асимметрии веса спортсменов''' df = load_data() swimmers = df[ df['Вид спорта'] == 'Swimming'] return swimmers['Вес'].skew()

0.23441459903001483

К счастью, эта асимметрия может быть эффективным образом смягчена путем взятия логарифма веса при помощи функции библиотеки numpy np.log:

def ex_3_5(): '''Визуализация разброса значений веса спортсменов на полулогарифмической гистограмме с целью удаления  асимметрии''' df = load_data() df['Вес'].apply(np.log).hist(bins=20) plt.xlabel('Логарифмический вес') plt.ylabel('Частота') plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

Теперь данные намного ближе к нормальному распределению. Из этого следует, что вес распределяется согласно логнормальному распределению.

Логнормальное распределение

Логнормальное распределение это распределение набора значений, чей логарифм нормально распределен. Основание логарифма может быть любым положительным числом за исключением единицы. Как и нормальное распределение, логнормальное распределение играет важную роль для описания многих естественных явлений.

Логарифм показывает степень, в которую должно быть возведено фиксированное число (основание) для получения данного числа. Изобразив логарифмы на графике в виде гистограммы, мы показали, что эти степени приближенно нормально распределены. Логарифмы обычно берутся по основанию 10 или основанию e, трансцендентному числу, приближенно равному 2.718. В функции библиотеки numpy np.log и ее инверсии np.exp используется основание e. Выражение log_e также называется натуральным логарифмом, или ln, из-за свойств, делающих его особенно удобным в исчислении.

Логнормальное распределение обычно имеет место в процессах роста, где темп роста не зависит от размера. Этот феномен известен как закон Джибрэта, который был cформулирован в 1931 г. Робертом Джибрэтом, заметившим, что он применим к росту фирм. Поскольку темп роста пропорционален размеру, более крупные фирмы демонстрируют тенденцию расти быстрее, чем фирмы меньшего размера.

^{Нормальное распределение случается в ситуациях, где много мелких колебаний, или вариаций, носит суммирующий эффект, тогда как логнормальное распределение происходит там, где много мелких вариаций имеет мультипликативный эффект.}

С тех пор выяснилось, что закон Джибрэта применим к большому числу ситуаций, включая размеры городов и, согласно обширному математическому ресурсу Wolfram MathWorld, к количеству слов в предложениях шотландского писателя Джорджа Бернарда Шоу.

В остальной части этой серии постов мы будем использовать натуральный логарифм веса спортсменов, чтобы наши данные были приближенно нормально распределены. Мы выберем популяцию спортсменов примерно с одинаковыми типами телосложения, к примеру, олимпийских пловцов.

Визуализация корреляции

Один из самых быстрых и самых простых способов определить наличие корреляции между двумя переменными состоит в том, чтобы рассмотреть их на графике рассеяния. Мы отфильтруем данные, выбрав только пловцов, и затем построим график роста относительно веса спортсменов:

def swimmer_data(): '''Загрузка данных роста и веса только олимпийских пловцов''' df = load_data() return df[df['Вид спорта'] == 'Swimming'].dropna()def ex_3_6(): '''Визуализация корреляции между ростом и весом''' df = swimmer_data() xs = df['Рост, см'] ys = df['Вес'].apply( np.log ) pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True) plt.xlabel('Рост, см.') plt.ylabel('Логарифмический вес') plt.show()

Этот пример сгенерирует следующий ниже график:

Результат ясно показывает, что между этими двумя переменными имеется связь. График имеет характерно смещенную эллиптическую форму двух коррелируемых, нормально распределенных переменных с центром вокруг среднего значения. Следующая ниже диаграмма сравнивает график рассеяния с распределениями вероятностей роста и логарифма веса:

Точки, близко расположенные к хвосту одного распределения, также демонстрируют тенденцию близко располагаться к тому же хвосту другого распределения, и наоборот. Таким образом, между двумя распределениями существует связь, которую в ближайших нескольких разделах мы покажем, как определять количественно. Впрочем, если мы внимательно посмотрим на предыдущий график рассеяния, то увидим, что из-за округления измерений точки уложены в столбцы и строки (в см. и кг. соответственно для роста и веса). Там, где это происходит, иногда желательно внести в данные искажения, которые также называются сдвигом или джиттером с тем, чтобы яснее показать силу связи. Без генерирования джиттера (в виде случайных отклонений) может оказаться, что, то, что по внешнему виду составляет одну точку, фактически представляет много точек, которые обозначены одинаковой парой значений. Внесение нескольких случайных помех делает эту ситуацию вряд ли возможной.

Генерирование джиттера

Поскольку каждое значение округлено до ближайшего сантиметра или килограмма, то значение, записанное как 180 см, на самом деле может быть каким угодно между 179.5 и 180.5 см, тогда как значение 80 кг на самом деле может быть каким угодно между 79.5 и 80.5 кг. Для создания случайных искажений, мы можем добавить случайные помехи в каждую точку данных роста в диапазоне между -0.5 и 0.5 и в том же самом диапазоне проделать с точками данных веса (разумеется, это нужно cделать до того, как мы возьмем логарифм значений веса):

def jitter(limit): '''Генератор джиттера (произвольного сдвига точек данных)''' return lambda x: random.uniform(-limit, limit) + xdef ex_3_7(): '''Визуализация корреляции между ростом и весом с джиттером''' df = swimmer_data() xs = df['Рост, см'].apply(jitter(0.5)) ys = df['Вес'].apply(jitter(0.5)).apply(np.log) pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True) plt.xlabel('Рост, см.') plt.ylabel('Логарифмический вес') plt.show()

График с джиттером выглядит следующим образом:

Как и в случае с внесением прозрачности в график рассеяния в первой серии постов об описательной статистике, генерирование джиттера это механизм, который обеспечивает исключение несущественных факторов, таких как объем данных или артефакты округления, которые могут заслонить от нас возможность увидеть закономерности в данных.

Ковариация

Одним из способов количественного определения силы связи между двумя переменными является их ковариация. Она измеряет тенденцию двух переменных изменяться вместе.

Если у нас имеется два ряда чисел, X и Y, то их отклонения от среднего значения составляют:

Здесь x_i это значение X с индексом i, y_i значение Y с индексом i, x среднее значение X, и y среднее значение Y. Если X и Y проявляют тенденцию изменяться вместе, то их отклонения от среднего будет иметь одинаковый знак: отрицательный, если они меньше среднего, положительный, если они больше среднего. Если мы их перемножим, то произведение будет положительным, когда у них одинаковый знак, и отрицательным, когда у них разные знаки. Сложение произведений дает меру тенденции этих двух переменных отклоняться от среднего значения в одинаковом направлении для каждой заданной выборки.

Ковариация определяется как среднее этих произведений:

На чистом Python ковариация вычисляется следующим образом:

def covariance(xs, ys): '''Вычисление ковариации (несмещенная, т.е. n-1)''' dx = xs - xs.mean()  dy = ys - ys.mean() return (dx * dy).sum() / (dx.count() - 1)

В качестве альтернативы, мы можем воспользоваться функцией pandas cov:

df['Рост, см'].cov(df['Вес'])

1.3559273321696459

Ковариация роста и логарифма веса для наших олимпийских пловцов равна 1.356, однако это число сложно интерпретировать. Единицы измерения здесь представлены произведением единиц на входе.

По этой причине о ковариации редко сообщают как об отдельной сводной статистике. Сделать число более понятным можно, разделив отклонения на произведение стандартных отклонений. Это позволяет трансформировать единицы измерения в стандартные оценки и ограничить выход числом в диапазоне между -1 и +1. Этот результат называется корреляцией Пирсона.

^{Стандартная оценка, англ. standard score, также z-оценка это относительное число стандартных отклонений, на которые значение переменной отстоит от среднего значения. Положительная оценка показывает, что переменная находится выше среднего, отрицательная ниже среднего. Это безразмерная величина, получаемая при вычитании популяционного среднего из индивидуальных значений и деления разности на популяционное стандартное отклонение.}

Корреляция Пирсона

Корреляция Пирсона часто обозначается переменной rи вычисляется следующим образом, где отклонения от среднего dx_iи dy_iвычисляются как и прежде:

Поскольку для переменных X и Y стандартные отклонения являются константными, уравнение может быть упрощено до следующего, где _xи _y это стандартные отклонения соответственно X и Y:

В таком виде формула иногда упоминается как коэффициент корреляции смешанных моментов Пирсона или попросту коэффициент корреляции и, как правило, обозначается буквой r.

Ранее мы уже написали функции для вычисления стандартного отклонения. В сочетании с нашей функцией с вычислением ковариации получится следующая реализация корреляции Пирсона:

def variance(xs): '''Вычисление корреляции, несмещенная дисперсия при n <= 30''' x_hat = xs.mean() n = xs.count() n = n - 1 if n in range( 1, 30 ) else n  return sum((xs - x_hat) ** 2) / ndef standard_deviation(xs): '''Вычисление стандартного отклонения''' return np.sqrt(variance(xs))def correlation(xs, ys):  '''Вычисление корреляции''' return covariance(xs, ys) / (standard_deviation(xs) *  standard_deviation(ys))

В качестве альтернативы мы можем воспользоваться функцией pandas corr:

df['Рост, см'].corr(df['Вес'])

Поскольку стандартные оценки безразмерны, то и коэффициент корреляции rтоже безразмерен. Если rравен -1.0 либо 1.0, то переменные идеально антикоррелируют либо идеально коррелируют.

Правда, если r = 0, то с необходимостью вовсе не следует, что переменные не коррелируют. Корреляция Пирсона измеряет лишь линейные связи. Как продемонстрировано на следующих графиках, между переменными может существовать еще некая нелинейная связь, которую rне объясняет:

Отметим, что корреляция центрального примера не определена, потому что стандартное отклонение y = 0. Поскольку наше уравнение для rсодержало бы деление ковариации на 0, то результат получается бессмысленным. В этом случае между переменными не может быть никакой корреляции; yвсегда будет иметь среднее значение. Простое обследование стандартных отклонений это подтвердит.

Мы можем вычислить коэффициент корреляции для данных роста и логарифма веса наших пловцов следующим образом:

def ex_3_8(): '''Вычисление корреляции средствами pandas на примере данных роста и веса''' df = swimmer_data() return df['Рост, см'].corr( df['Вес'].apply(np.log))

0.86748249283924894

В результате получим ответ 0.867, который количественно выражает сильную, положительную корреляцию, уже наблюдавшуюся нами на точечном графике.

Выборочный rи популяционный

Аналогично среднему значению и стандартному отклонению, коэффициент корреляции является сводной статистикой. Он описывает выборку; в данном случае, выборку спаренных значений: роста и веса. Коэффициент корреляции известной выборки обозначается буквой r, тогда как коэффициент корреляции неизвестной популяции обозначается греческой буквой (рхо).

Как мы убедились в предыдущей серии постов о тестировании гипотез, мы не должны исходить из того, что результаты, полученные в ходе измерения нашей выборки, применимы к популяции в целом. К примеру, наша популяция может состоять из всех пловцов всех недавних Олимпийских игр. И будет совершенно недопустимо обобщать, например, на другие олимпийские виды спорта, такие как тяжелая атлетика или фитнес-плавание.

Даже в допустимой популяции такой как пловцы, выступавшие на недавних Олимпийских играх, наша выборка коэффициента корреляции является всего лишь одной из многих потенциально возможных. То, насколько мы можем доверять нашему r, как оценке параметра , зависит от двух факторов:

• Размера выборки

• Величины r

Безусловно, чем больше выборка, тем больше мы ей доверяем в том, что она представляет всю совокупность в целом. Возможно, не совсем интуитивно очевидно, но величина тоже оказывает влияние на степень нашей уверенности в том, что выборка представляет параметр . Это вызвано тем, что большие коэффициенты вряд ли возникли случайным образом или вследствие случайной ошибки при отборе.

Проверка статистических гипотез

В предыдущей серии постов мы познакомились с проверкой статистических гипотез, как средством количественной оценки вероятности, что конкретная гипотеза (как, например, что две выборки взяты из одной и той же популяции) истинная. Чтобы количественно оценить вероятность, что корреляция существует в более широкой популяции, мы воспользуемся той же самой процедурой.

В первую очередь, мы должны сформулировать две гипотезы, нулевую гипотезу и альтернативную:

H₀ - это гипотеза, что корреляция в популяции нулевая. Другими словами, наше консервативное представление состоит в том, что измеренная корреляция целиком вызвана случайной ошибкой при отборе.

H₁ - это альтернативная возможность, что корреляция в популяции не нулевая. Отметим, что мы не определяем направление корреляции, а только что она существует. Это означает, что мы выполняем двустороннюю проверку.

Стандартная ошибка коэффициента корреляции rпо выборке задается следующей формулой:

Эта формула точна, только когда rнаходится близко к нулю (напомним, что величина влияет на нашу уверенность), но к счастью, это именно то, что мы допускаем согласно нашей нулевой гипотезы.

Мы можем снова воспользоваться t-распределением и вычислить t-статистику:

В приведенной формуле df это степень свободы наших данных. Для проверки корреляции степень свободы равна n - 2, где n это размер выборки. Подставив это значение в формулу, получим:

В итоге получим t-значение 102.21. В целях его преобразования в p-значение мы должны обратиться к t-распределению. Библиотека scipy предоставляет интегральную функцию распределения (ИФР) для t-распределения в виде функции stats.t.cdf, и комплементарной ей (1-cdf) функции выживания stats.t.sf. Значение функции выживания соответствует p-значению для односторонней проверки. Мы умножаем его на 2, потому что выполняем двустороннюю проверку:

def t_statistic(xs, ys): '''Вычисление t-статистики''' r = xs.corr(ys) # как вариант, correlation(xs, ys) df = xs.count() - 2 return r * np.sqrt(df / 1 - r ** 2)def ex_3_9(): '''Выполнение двухстороннего t-теста''' df = swimmer_data() xs = df['Рост, см'] ys = df['Вес'].apply(np.log) t_value = t_statistic(xs, ys) df = xs.count() - 2  p = 2 * stats.t.sf(t_value, df) # функция выживания  return {'t-значение':t_value, 'p-значение':p}

{'p-значение': 1.8980236317815443e-106, 't-значение': 25.384018200627057}

P-значение настолько мало, что в сущности равно 0, означая, что шанс, что нулевая гипотеза является истинной, фактически не существует. Мы вынуждены принять альтернативную гипотезу о существовании корреляции.

Интервалы уверенности

Установив, что в более широкой популяции, безусловно, существует корреляция, мы, возможно, захотим количественно выразить диапазон значений, внутри которого, как мы ожидаем, будет лежать параметр , вычислив для этого интервал уверенности. Как и в случае со средним значением в предыдущей серии постов, интервал уверенности для rвыражает вероятность (выраженную в %), что параметр популяции находится между двумя конкретными значениями.

Однако при попытке вычислить стандартную ошибку коэффициента корреляции возникает сложность, которой не было в случае со средним значением. Поскольку абсолютное значение коэффициента корреляции rне может превышать 1, распределение возможных выборок коэффициентов корреляции rсмещается по мере приближения rк пределу своего диапазона.

Приведенный выше график показывает отрицательно скошенное распределение r-выборок для параметра , равного 0.6.

К счастью, трансформация под названием z-преобразование Фишера стабилизирует дисперсию r по своему диапазону. Она аналогична тому, как наши данные о весе спортсменов стали нормально распределенными, когда мы взяли их логарифм.

Уравнение для z-преобразования следующее:

Стандартная ошибка z равна:

Таким образом, процедура вычисления интервалов уверенности состоит в преобразовании rв z с использованием z-преобразования, вычислении интервала уверенности в терминах стандартной ошибки SE_zи затем преобразовании интервала уверенности в r.

В целях вычисления интервала уверенности в терминах SE_z, мы можем взять число стандартных отклонений от среднего, которое дает нам требуемый уровень доверия. Обычно используют число 1.96, так как оно является числом стандартных отклонений от среднего, которое содержит 95% площади под кривой. Другими словами, 1.96 стандартных ошибок от среднего значения выборочного rсодержит истинную популяционную корреляцию с 95%-ой определенностью.

Мы можем убедиться в этом, воспользовавшись функцией scipy stats.norm.ppf. Она вернет стандартную оценку, связанную с заданной интегральной вероятностью в условиях односторонней проверки.

Однако, как показано на приведенном выше графике, мы хотели бы вычесть ту же самую величину, т.е. 2.5%, из каждого хвоста с тем, чтобы 95%-й интервал уверенности был центрирован на нуле. Для этого при выполнении двусторонней проверки нужно просто уменьшить разность наполовину и вычесть результат из 100%. Так что, требуемый уровень доверия в 95% означает, что мы обращаемся к критическому значению 97.5%:

def critical_value(confidence, ntails): # ДИ и число хвостов '''Расчет критического значения путем вычисления квантиля и получения  для него нормального значения''' lookup = 1 - ((1 - confidence) / ntails)  return stats.norm.ppf(lookup, 0, 1) # mu=0, sigma=1critical_value(0.95, 2)

1.959963984540054

Поэтому наш 95%-й интервал уверенности в z-пространстве для задается следующей формулой:

Подставив в нашу формулу z_rи SE_z, получим:

Для r=0.867и n=859она даст нижнюю и верхнюю границу соответственно 1.137 и 1.722. В целях их преобразования из z-оценок в r-значения, мы используем следующее обратное уравнение z-преобразования:

Преобразования и интервал уверенности можно вычислить при помощи следующего исходного кода:

def z_to_r(z): '''Преобразование z-оценки обратно в r-значение''' return (np.exp(z*2) - 1) / (np.exp(z*2) + 1)def r_confidence_interval(crit, xs, ys):  '''Расчет интервала уверенности для критического значения и данных''' r = xs.corr(ys) n = xs.count() zr = 0.5 * np.log((1 + r) / (1 - r))  sez = 1 / np.sqrt(n - 3) return (z_to_r(zr - (crit * sez))), (z_to_r(zr + (crit * sez)))def ex_3_10(): '''Расчет интервала уверенности на примере данных роста и веса''' df = swimmer_data() X = df['Рост, см'] y = df['Вес'].apply(np.log) interval = r_confidence_interval(1.96, X, y)  print('Интервал уверенности (95%):', interval)

Интервал уверенности (95%): (0.8499088588880347, 0.8831284878884087)

В результате получаем 95%-й интервал уверенности для , расположенный между 0.850 и 0.883. Мы можем быть абсолютно уверены в том, что в более широкой популяции олимпийских пловцов существует сильная положительная корреляция между ростом и весом.

В следующем посте, посте 2, будет рассмотрена сама тема серии постов - регрессия и приемы оценивания ее качества.

Предыдущий пост см. здесь.

Регрессия

Хотя, возможно, и полезно знать, что две переменные коррелируют, мы не можем использовать лишь одну эту информацию для предсказания веса олимпийских пловцов при наличии данных об их росте или наоборот. При установлении корреляции мы измерили силу и знак связи, но не наклон, т.е. угловой коэффициент. Для генерирования предсказания необходимо знать ожидаемый темп изменения одной переменной при заданном единичном изменении в другой.

Мы хотели бы вывести уравнение, связывающее конкретную величину одной переменной, так называемой независимой переменной, с ожидаемым значением другой, зависимой переменной. Например, если наше линейное уравнение предсказывает вес при заданном росте, то рост является нашей независимой переменной, а вес зависимой.

^{Описываемые этими уравнениями линии называются линиями регрессии. Этот Термин был введен британским эрудитом 19-ого века сэром Фрэнсисом Гэлтоном. Он и его студент Карл Пирсон, который вывел коэффициент корреляции, в 19-ом веке разработали большое количество методов, применяемых для изучения линейных связей, которые коллективно стали известны как методы регрессионного анализа.}

Вспомним, что из корреляции не следует причинная обусловленность, причем термины зависимый и независимый не означают никакой неявной причинной обусловленности. Они представляют собой всего лишь имена для входных и выходных математических значений. Классическим примером является крайне положительная корреляция между числом отправленных на тушение пожара пожарных машин и нанесенным пожаром ущербом. Безусловно, отправка пожарных машин на тушение пожара сама по себе не наносит ущерб. Никто не будет советовать сократить число машин, отправляемых на тушение пожара, как способ уменьшения ущерба. В подобных ситуациях мы должны искать дополнительную переменную, которая была бы связана с другими переменными причинной связью и объясняла корреляцию между ними. В данном примере это может быть размер пожара. Такие скрытые причины называются спутывающими переменными, потому что они искажают нашу возможность определять связь между зависимыми переменными.

Линейные уравнения

Две переменные, которые мы можем обозначить как xи y, могут быть связаны друг с другом строго или нестрого. Самая простая связь между независимой переменной xи зависимой переменной yявляется прямолинейной, которая выражается следующей формулой:

Здесь значения параметров aи bопределяют соответственно точную высоту и крутизну прямой. Параметр aназывается пересечением с вертикальной осью или константой, а b градиентом, наклоном линии или угловым коэффициентом. Например, в соотнесенности между температурными шкалами по Цельсию и по Фаренгейту a = 32и b = 1.8. Подставив в наше уравнение значения aи b, получим:

Для вычисления 10С по Фаренгейту мы вместо xподставляем 10:

Таким образом, наше уравнение сообщает, что 10С равно 50F, и это действительно так. Используя Python и возможности визуализации pandas, мы можем легко написать функцию, которая переводит градусы из Цельсия в градусы Фаренгейта и выводит результат на график:

'''Функция перевода из градусов Цельсия в градусы Фаренгейта'''celsius_to_fahrenheit = lambda x: 32 + (x * 1.8)def ex_3_11(): '''График линейной зависимости температурных шкал''' df = pd.DataFrame({'C':s, 'F':s.map(celsius_to_fahrenheit)}) df.plot('C', 'F', legend=False, grid=True) plt.xlabel('Градусы Цельсия') plt.ylabel('Градусы Фаренгейта') plt.show()

Этот пример сгенерирует следующий ниже линейный график:

Обратите внимание, как синяя линия пересекает 0 на шкале Цельсия при величине 32 на шкале Фаренгейта. Пересечение a это значение y, при котором значение xравно 0.

Наклон линии с неким угловым коэффициентом определяется параметром b; в этом уравнении его значение близко к 2. Как видно, диапазон шкалы Фаренгейта почти вдвое шире диапазона шкалы Цельсия. Другими словами, прямая устремляется вверх по вертикали почти вдвое быстрее, чем по горизонтали.

Остатки

К сожалению, немногие связи столь же чистые, как перевод между градусами Цельсия и Фаренгейта. Прямолинейное уравнение редко позволяет нам определять yстрого в терминах x. Как правило, будет иметься ошибка, и, таким образом, уравнение примет следующий вид:

Здесь, это ошибка или остаточный член, обозначающий расхождение между значением, вычисленным параметрами aи bдля данного значения xи фактическим значением y. Если предсказанное значение y это y, то ошибка это разность между обоими:

Такая ошибка называется остатком. Остаток может возникать из-за случайных факторов, таких как погрешность измерения, либо неслучайных факторов, которые неизвестны. Например, если мы пытаемся предсказать вес как функцию роста, то неизвестные факторы могут состоять из диеты, уровня физической подготовки и типа телосложения (либо просто эффекта округления до самого близкого килограмма).

Если для aи bмы выберем неидеальные параметры, то остаток для каждого xбудет больше, чем нужно. Из этого следует, что параметры, которые мы бы хотели найти, должны минимизировать остатки во всех значениях xи y.

Обычные наименьшие квадраты

Для того, чтобы оптимизировать параметры линейной модели, мы бы хотели создать функцию стоимости, так называемую функцией потери, которая количественно выражает то, насколько близко наши предсказания укладывается в данные. Мы не можем просто взять и просуммировать положительные и отрицательные остатки, потому что даже самые большие остатки обнулят друг друга, если их знаки противоположны.

Прежде, чем вычислить сумму, мы можем возвести значения в квадрат, чтобы положительные и отрицательные остатки учитывались в стоимости. Возведение в квадрат также создает эффект наложения большего штрафа на большие ошибки, чем на меньшие ошибки, но не настолько много, чтобы самый большой остаток всегда доминировал.

Выражаясь в терминах задачи оптимизации, мы стремимся выявить коэффициенты, которые минимизируют сумму квадратов остатков. Этот метод называется обычными наименьшими квадратами, от англ. Ordinary Least Squares (OLS), и формула для вычисления наклона линии регрессии по указанному методу выглядит так:

Хотя она выглядит сложнее предыдущих уравнений, на самом деле, эта формула представляет собой всего лишь сумму квадратов остатков, деленную на сумму квадратов отклонений от среднего значения. В данном уравнении используется несколько членов из других уравнений, которые уже рассматривались, и мы можем его упростить, приведя к следующему виду:

Пересечение (a) это член, позволяющий прямой с заданным наклоном проходить через среднее значение X и Y:

Значения aи b это коэффициенты, получаемые в результате оценки методом обычных наименьших квадратов.

Наклон и пересечение

Мы уже рассматривали функции covariance, variance и mean, которые нужны для вычисления наклона прямой и точки пересечения для данных роста и веса пловцов. Поэтому вычисление наклона и пересечения имеют тривиальный вид:

def slope(xs, ys): '''Вычисление наклона линии (углового коэффициента)''' return xs.cov(ys) / xs.var()def intercept(xs, ys):  '''Вычисление точки пересечения (с осью Y)''' return ys.mean() - (xs.mean() * slope(xs, ys))def ex_3_12(): '''Вычисление пересечения и наклона (углового коэффициента)  на примере данных роста и веса''' df = swimmer_data() X = df['Рост, см'] y = df['Вес'].apply(np.log) a = intercept(X, y) b = slope(X, y)  print('Пересечение: %f, наклон: %f' % (a,b))

Пересечение: 1.691033, наклон: 0.014296

В результате будет получен наклон приблизительно 0.0143 и пересечение приблизительно 1.6910.

Интерпретация

Величина пересечения это значение зависимой переменной (логарифмический вес), когда независимая переменная (рост) равна нулю. Для получения этого значения в килограммах мы можем воспользоваться функцией np.exp, обратной для функции np.log. Наша модель дает основания предполагать, что вероятнее всего вес олимпийского пловца с нулевым ростом будет 5.42 кг. Разумеется, такое предположение лишено всякого смысла, к тому же экстраполяция за пределы границ тренировочных данных является не самым разумным решением.

Величина наклона показывает, насколько yизменяется для каждой единицы изменения в x. Модель исходит из того, что каждый дополнительный сантиметр роста прибавляет в среднем 1.014 кг. веса олимпийских пловцов. Поскольку наша модель основывается на данных о всех олимпийских пловцах, она представляет собой усредненный эффект от увеличения в росте на единицу без учета любого другого фактора, такого как возраст, пол или тип телосложения.

Визуализация

Результат линейного уравнения можно визуализировать при помощи имплементированной ранее функции regression_line и простой функции от x, которая вычисляет yна основе коэффициентов aи b.

'''Функция линии регрессии'''regression_line = lambda a, b: lambda x: a + (b * x) # вызовы fn(a,b)(x)def ex_3_13(): '''Визуализация линейного уравнения на примере данных роста и веса''' df = swimmer_data() X = df['Рост, см'].apply( jitter(0.5) ) y = df['Вес'].apply(np.log) a, b = intercept(X, y), slope(X, y)  ax = pd.DataFrame(np.array([X, y]).T).plot.scatter(0, 1, s=7) s = pd.Series(range(150,210)) df = pd.DataFrame( {0:s, 1:s.map(regression_line(a, b))} )  df.plot(0, 1, legend=False, grid=True, ax=ax) plt.xlabel('Рост, см.') plt.ylabel('Логарифмический вес') plt.show()

Функция regression_line возвращает функцию от x, которая вычисляет a + bx.

Указанная функция может также использоваться для вычисления каждого остатка, показывая степень, с которой наша оценка yотклоняется от каждого измеренного значения y.

def residuals(a, b, xs, ys): '''Вычисление остатков''' estimate = regression_line(a, b) # частичное применение return pd.Series( map(lambda x, y: y - estimate(x), xs, ys) )constantly = lambda x: 0def ex_3_14(): '''Построение графика остатков на примере данных роста и веса''' df = swimmer_data() X = df['Рост, см'].apply( jitter(0.5) ) y = df['Вес'].apply(np.log) a, b = intercept(X, y), slope(X, y)  y = residuals(a, b, X, y) ax = pd.DataFrame(np.array([X, y]).T).plot.scatter(0, 1, s=12) s = pd.Series(range(150,210)) df = pd.DataFrame( {0:s, 1:s.map(constantly)} )  df.plot(0, 1, legend=False, grid=True, ax=ax) plt.xlabel('Рост, см.') plt.ylabel('Остатки') plt.show()

График остатков это график, который показывает остатки на оси Y и независимую переменную на оси X. Если точки на графике остатков разбросаны произвольно по обе стороны от горизонтальной оси, то линейная модель хорошо подогнана к нашим данным:

За исключением нескольких выбросов на левой стороне графика, график остатков, по-видимому, показывает, что линейная модель хорошо подогнана к данным. Построение графика остатков имеет важное значение для получения подтверждения, что линейная модель применима. В линейной модели используются некоторые допущения относительно данных, которые при их нарушении делают не валидными модели, которые вы строите.

Допущения

Первостепенное допущение линейной регрессии состоит в том, что, безусловно, существует линейная зависимость между зависимой и независимой переменной. Кроме того, остатки не должны коррелировать друг с другом либо с независимой переменной. Другими словами, мы ожидаем, что ошибки будут иметь нулевое среднее и постоянную дисперсию по отношению к зависимой и независимой переменной. График остатков позволяет быстро устанавливать, является ли это действительно так.

Левая сторона нашего графика имеет более крупные значения остатков, чем правая сторона. Это соответствует большей дисперсии веса среди более низкорослых спортсменов. Когда дисперсия одной переменной изменяется относительно другой, говорят, что переменные гетероскедастичны, т.е. их дисперсия неоднородна. Этот факт представляет в регрессионном анализе проблему, потому что делает не валидным допущение в том, что модельные ошибки не коррелируют и нормально распределены, и что их дисперсии не варьируются вместе с моделируемыми эффектами.

Гетероскедастичность остатков здесь довольно мала и особо не должна повлиять на качество нашей модели. Если дисперсия на левой стороне графика была бы более выраженной, то она привела бы к неправильной оценке дисперсии методом наименьших квадратов, что в свою очередь повлияло бы на выводы, которые мы делаем, основываясь на стандартной ошибке.

Качество подгонки и R-квадрат

Хотя из графика остатков видно, что линейная модель хорошо вписывается в данные, т.е. хорошо к ним подогнана, было бы желательно количественно измерить качество этой подгонки. R², или R-квадрат, варьируется в интервале между 0 и 1 и обозначает объяснительную мощность линейной регрессионной модели. Он вычисляет объясненную долю изменчивости в зависимой переменной.

Обычно, чем ближе R²к 1, тем лучше линия регрессии подогнана к точкам данных и больше изменчивости в Y объясняется независимой переменной X. R² можно вычислить с помощью следующей ниже формулы:

Здесь var() это дисперсия остатков и var(Y) дисперсия в Y. В целях понимания смысла этой формулы допустим, что вы пытаетесь угадать чей-то вес. Если вам больше ничего неизвестно об испытуемых, то наилучшей стратегией будет угадывать среднее значение весовых данных внутри популяции в целом. Таким путем средневзвешенная квадратичная ошибка вашей догадки в сравнении с истинным весом будет var(Y), т.е. дисперсией данных веса в популяции.

Но если бы я сообщил вам их рост, то в соответствии с регрессионной моделью вы бы предположили, что a + bx. В этом случае вашей средневзвешенной квадратичной ошибкой было бы или дисперсия остатков модели.

Компонент формулы var()/var(Y) это соотношение средневзвешенной квадратичной ошибки с объяснительной переменной и без нее, т. е. доля изменчивости, оставленная моделью без объяснения. Дополнение R²до единицы это доля изменчивости, объясненная моделью.

^{Как и в случае с r, низкий R2не означает, что две переменные не коррелированы. Просто может оказаться, что их связь не является линейной.}

Значение R² описывает качество подгонки линии регрессии к данным. Оптимально подогнанная линия это линия, которая минимизирует значение R². По мере удаления либо приближения от своих оптимальных значений R² всегда будет расти.

Левый график показывает дисперсию модели, которая всегда угадывает среднее значение для , правый же показывает меньшие по размеру квадраты, связанные с остатками, которые остались необъясненными моделью f. С чисто геометрической точки зрения можно увидеть, как модель объяснила большинство дисперсии в y. Приведенный ниже пример вычисляет R²путем деления дисперсии остатков на дисперсию значений y:

def r_squared(a, b, xs, ys): '''Рассчитать коэффициент детерминации (R-квадрат)''' r_var = residuals(a, b, xs, ys).var()  y_var = ys.var() return 1 - (r_var / y_var)def ex_3_15(): '''Рассчитать коэффициент R-квадрат  на примере данных роста и веса''' df = swimmer_data() X = df['Рост, см'].apply( jitter(0.5) ) y = df['Вес'].apply(np.log) a, b = intercept(X, y), slope(X, y) return r_squared(a, b, X, y)

0.75268223613272323

В результате получим значение 0.753. Другими словами, более 75% дисперсии веса пловцов, выступавших на Олимпийских играх 2012 г., можно объяснить ростом.

В случае простой регрессионной модели (с одной независимой переменной), связь между коэффициентом детерминации R² и коэффициентом корреляции rявляется прямолинейной:

Коэффициент корреляции rможет означать, что половина изменчивости в переменной Y объясняется переменной X, но фактически R² составит 0.5², т.е. 0.25.

Множественная линейная регрессия

Пока что в этой серии постов мы видели, как строится линия регрессии с одной независимой переменной. Однако, нередко желательно построить модель с несколькими независимыми переменными. Такая модель называется множественной линейной регрессией.

Каждой независимой переменной потребуется свой собственный коэффициент. Вместо того, чтобы для каждой из них пытаться подобрать букву в алфавите, зададим новую переменную (бета), которая будет содержать все наши коэффициенты:

Такая модель эквивалентна двухфакторной линейно-регрессионной модели, где ₁= a и ₂= b при условии, что x₁всегда гарантированно равен 1, вследствие чего ₁ это всегда константная составляющая, которая представляет наше пересечение, при этом x₁называется постоянным смещением уравнения регрессии, или членом смещения.

Обобщив линейное уравнение в терминах , его легко расширить на столько коэффициентов, насколько нам нужно:

Каждое значение от x₁до x_nсоответствует независимой переменной, которая могла бы объяснить значение y. Каждое значение от ₁до _nсоответствует коэффициенту, который устанавливает относительный вклад независимой переменной.

Простая линейная регрессия преследовала цель объяснить вес исключительно с точки зрения роста, однако объяснить вес людей помогает много других факторов: их возраст, пол, питание, тип телосложения. Мы располагаем сведениями о возрасте олимпийских пловцов, поэтому мы смогли бы построить модель, которая учитывает и эти дополнительные данные.

До настоящего момента мы предоставляли независимую переменную в виде одной последовательности значений, однако при наличии двух и более параметров нам нужно предоставлять несколько значений для каждого x. Мы можем воспользоваться функциональностью библиотеки pandas, чтобы выбрать два и более столбцов и управлять каждым как списком, но есть способ получше: матрицы.

Темой следующего поста, поста 3, будут матричные операции, нормальное уравнение и коллинеарность.

Предыдущий пост см. здесь.

Предсказание

В заключении, мы подходим к одному из наиболее важных применений линейной регрессии: предсказанию. Мы натренировали модель способную предсказывать вес олимпийских пловцов при наличии данных об их росте, половой принадлежности и годе рождения.

9-кратный олимпийский чемпион по плаванию Марк Шпитц завоевал 7 золотых медалей на Олимпийских играх 1972 г. Он родился в 1950 г. и, согласно веб-страницы Википедии, имеет рост 183 см. и вес 73 кг. Посмотрим, что наша модель предсказывает в отношении его веса.

Наша множественная регрессионная модель требует предоставить эти значения в матричной форме. Каждый параметр нужно передать в том порядке, в котором модель усвоила признаки, чтобы применить правильный коэффициент. После смещения вектор признаков должен содержать рост, пол и год рождения в тех же единицах измерения, в которых модель была натренирована:

Матрица содержит коэффициенты для каждого из этих признаков:

Предсказанием модели будет сумма произведений коэффициентов и признаков xв каждой строке:

Поскольку матричное умножение производит каждый элемент в результате соответствующего сложения произведений строк и столбцов каждой матрицы, получение результата попросту сводится к перемножению транспонированной матрицы с вектором xspitz.

Напомним, что размерность итоговой матрицы будет числом строк из первой матрицы и числом столбцов из второй матрицы:

Здесь ^Tx это произведение матрицы размера 1 nи матрицы размера n 1. Результатом является матрица размера 1 1:

Исходный код вычислений очень прост:

def predict(coefs, x):     '''функция предсказания'''    return np.matmul(coefs, x.values) 

def ex_3_29():    '''Вычисление ожидаемого веса спортсмена'''    df = swimmer_data()    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int)     df['Год рождения'] = df['Дата рождения'].map(str_to_year)    X = df[['Рост, см', 'бин_Пол', 'Год рождения']]     X.insert(0, 'константа', 1.0)    y = df['Вес'].apply(np.log)     beta = linear_model(X, y)    xspitz = pd.Series([1.0, 183, 1, 1950]) # параметры Марка Шпитца    return np.exp( predict(beta, xspitz) )  

84.20713139038605

Этот пример вернет число 84.21, которое соответствует ожидаемому весу 84.21 кг. Это намного тяжелее зарегистрированного веса Марка Шпитца 73 кг. Наша модель, похоже, не сработала как надо.

Интервал уверенности для конкретного предсказания

Ранее мы рассчитали интервалы уверенности для параметров популяции. Мы можем также построить интервал уверенности для отдельно взятого предсказания, именуемый интервалом предсказания, или интервалом изменений предсказанной величины. Интервал предсказания количественно определяет размер неопределенности в предсказании, предлагая минимальное и максимальное значение, между которыми истинное значение ожидаемо попадает с определенной вероятностью. Интервал предсказания для yшире интервала уверенности для параметра популяции, такого как например, среднее значение . Это объясняется тем, что в интервале уверенности учитывается только неопределенность при оценке среднего значения, в то время как интервал предсказания также должен принимать в расчет дисперсию в y^{от среднего значения.}

^{Линии интервала предсказания находятся дальше от оптимально подогнанной линии, чем линии интервала уверенности, и тем дальше, чем больше имеется точек данных. 95%-ый интервал предсказания это площадь, в которую ожидаемо попадают 95% всех точек данных. В отличие от него, 95%-й интервал уверенности это площадь, которая с 95%-й вероятностью содержит истинную линию регрессии.}

Приведенный выше график показывает связь между внешним интервалом предсказания и внутренним интервалом уверенности. Мы можем вычислить интервал предсказания, используя следующую ниже формулу:

Здесь y_p это предсказание, плюс или минус интервал. Мы пользуемся t-распределением, где степень свободы равна n - p, т.е. размер выборки минус число параметров. Это та же самая формула, которая ранее применялась при вычислении F-тестов. Хотя указанная формула, возможно, пугает своей сложностью, она относительно прямолинейно транслируется в исходный код, показанный в следующем ниже примере, который вычисляет 95%-ый интервал предсказания.

def prediction_interval(x, y, xp):    '''Вычисление интервала предсказания'''    xtx    = np.matmul(x.T, np.asarray(x))    xtxi   = np.linalg.inv(xtx)      xty    = np.matmul(x.T, np.asarray(y))     coefs  = linear_model(x, y)     fitted = np.matmul(x, coefs)    resid  = y - fitted    rss    = resid.dot(resid)      n      = y.shape[0]  # строки    p      = x.shape[1]  # столбцы    dfe    = n - p     mse    = rss / dfe    se_y   = np.matmul(np.matmul(xp.T, xtxi), xp)    t_stat = np.sqrt(mse * (1 + se_y))         # t-статистика    intl   = stats.t.ppf(0.975, dfe) * t_stat       yp     = np.matmul(coefs.T, xp)    return np.array([yp - intl, yp + intl])

Поскольку t-статистика параметризуется степенью свободы ошибки, она принимает в расчет присутствующую в модели неопределенность.

Если вместо интервала предсказания потребуется рассчитать интервал уверенности для среднего значения, мы попросту можем опустить прибавление единицы к se_y при вычислении t-статистики t_stat.

Приведенный выше исходный код можно использовать для генерирования следующего ниже графика, который показывает, как интервал предсказания изменяется вместе со значением независимой переменной:

В приведенном выше графике модель тренируется на выборке из 5 элементов и показывает, как 95%-й интервал предсказания увеличивается по мере нашего движения дальше от среднего роста. Применив приведенную выше формулу к Марку Шпитцу, в результате получим следующее:

def ex_3_30():    '''Интервал предсказания       применительно к данным о Марке Шпитце'''    df = swimmer_data()    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int)     df['Год рождения'] = df['Дата рождения'].map(str_to_year)    X = df[['Рост, см', 'бин_Пол', 'Год рождения']]     X.insert(0, 'константа', 1.0)    y = df['Вес'].apply(np.log)     xspitz = pd.Series([1.0, 183, 1, 1950])  # данные М.Шпитца    return np.exp( prediction_interval(X, y, xspitz) )

array([72.74964444, 97.46908087])

Этот пример возвращает диапазон между 72.7 и 97.4 кг., который как раз включает в себя вес Марка 73 кг., поэтому наше предсказание находится в пределах 95%-ого интервала предсказания. Правда оно лежит неудобно близко к границам диапазона.

Границы действия модели

Марк Шпитц родился в 1950 г., за несколько десятилетий до рождения самых возрастных пловцов на Олимпийских играх 2012 г. Пытаясь предсказать вес Марка, используя его год рождения, мы виноваты в том, что пытаемся экстраполировать слишком далеко за пределы наших тренировочных данных. Мы превысили границы действия нашей модели.

Есть и вторая причина, по которой этот результат сопряжен с проблемами. Наши данные полностью основаны на пловцах, которые соревнуются в наши дни на уровне международных стандартов, тогда как Марк не участвовал в соревнованиях уже много лет. Другими словами, Марк теперь не является частью популяции, на которой мы натренировали нашу модель. В целях исправления обеих этих проблем нам нужно поискать подробную информацию о Марке в 1979 г., когда он еще участвовал в соревнованиях по плаванию.

Согласно данным, в 1972 г. 22-летний Марк Шпитц имел рост 185 см. и весил 79 кг.

^{Отбор верных признаков это одна из самых важных предпосылок для получения хороших результатов от любого алгоритма предсказания.}

Следует отбирать признаки не только, опираясь на их предсказательную силу, но и на их актуальность для моделируемой области.

Окончательная модель

Хотя наша модель имеет чуть более низкий R², давайте натренируем ее повторно, взяв в качестве признака возраст вместо года рождения. Это позволит нам легко предсказывать веса для прошлых и будущих ранее не встречавшихся данных, т.к. она ближе моделирует переменную, которая, как мы полагаем, имеет причинно-следственную связь с весом.

Модель произведет приблизительно с такими значениями:

Наши признаки для Марка на играх 1972 г. таковы:

Эти признаки можно использовать для предсказания его соревновательного веса при помощи приведенного ниже примера:

def ex_3_32():    '''Окончательная модель для предсказания        соревновательного веса'''    df = swimmer_data()    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int)     X = df[['Рост, см', 'бин_Пол', 'Возраст']]     X.insert(0, 'константа', 1.0)    y = df['Вес'].apply(np.log)     beta = linear_model(X, y)    # предсказать вес Марка Шпитца    xspitz = pd.Series([1.0, 185, 1, 22])     return np.exp( predict(beta, xspitz) )

78.46882772630318

Пример возвращает число 78.47, т.е. предсказывает вес 78.47 кг. Теперь результат находится очень близко к истинному соревновательному весу Марка, равному 79 кг.

Примеры исходного кода для этого поста находятся в моемрепона Github. Все исходные данные взяты врепозиторииавтора книги.

Резюме

В этой серии постов мы узнали о том, как определять наличие общей линейной связи между двумя и более переменными. Мы увидели, как выражать силу их корреляции при помощи корреляции Пирсона rи насколько хорошо линейная модель объясняет вариацию при помощи коэффициентов детерминации R²и R². Мы также выполнили проверку статистических гипотез и рассчитали интервалы уверенности, чтобы вывести диапазон истинного популяционного параметра для корреляции.

Определив корреляцию между переменными, мы смогли построить предсказательную модель с использованием регрессии методом обычных наименьших квадратов и простых функций Python. Затем мы обобщили наш подход, воспользовавшись для этого матричной функциональностью библиотек pandas и numpy и нормальным уравнением. Эта простая модель продемонстрировала принципы машинного усвоения закономерностей путем определения модельных параметров , извлеченных из выборочных данных, которые могут использоваться для выполнения предсказаний. Наша модель смогла предсказать ожидаемый вес нового атлета, который полностью попадал в пределы интервала предсказания для истинного значения.

Проблемы, страхи и надежды традиционных банков отлично резюмировал Google. В своём тактичном предложении стать партнёрами, которое составил для них в прошлом году:

"Youre lagging in technology. Your current vendors are years behind. Consumers think youre irrelevant. Were hip, were cool, we have all the latest technologies, and boy have weve got data! Come partner with us on our new checking account!"

Перевод: Вы технологически отстаёте. Ваши вендоры застряли в прошлом. Клиенты считают вас устаревшими. Мы модные, мы классные, у нас все последние технологии и боже мой, сколько же у нас данных! Становитесь нашими партнёрами в новом банковском приложении!

Забавно, что 75% первых партнёров, похоже, даже не знали, о каком приложении речь! Это был Google Plex часть обновленного Google Pay с фичами вроде кэшбэков, специальных предложений и персональных финансовых советов.

Пример ярко показывает кризис привычной модели взаимодействия банк-клиент. Обе стороны ощущают нехватку в этой цепочке ещё одного звена. И это финтех.

Я уже несколько лет работаю с банками, ниже моё субъективное представление о важности и возможностях их коллаборации с ИТ-шниками + статистика и ссылки на источники.

Банки вынуждены открывать доступ для финтеха

Как говорится из каждого утюга после пандемии мир никогда не будет прежним. А вот в нашем случае вирус выступил катализатором уже происходящих процессов.

Ожидания потребителей по отношению к банкам изменились ещё в доковидную эпоху. Прошлогоднее исследование Ernst&Young показало, что клиенты ждут от банкинга:

1. Прозрачности и возможности доверять

2. Максимальной персонализации

3. Защиты данных

4. Мгновенного удовлетворения запросов

5. Помощи с достижением личных целей

Иными словами, того, что умеют делать agile финтехи и не умеют мастодонты старой закалки.

Accenture, в свою очередь, проанализировали изменение поведения людей в связи с локдаунами. Тенденции подтвердились. Клиентов банков поделили на четыре группы пользователей:

1. Прагматики ценят полезность услуг вне зависимости от используемого канала

2. Традиционалисты предпочитают общаться с людьми, а не онлайн

3. Пионеры любители инноваций, первыми пробуют новинки

4. Скептики разбираются в технологиях, но часто недовольны качеством сервисов

И вот несколько выводов:

• доверие к банкам со стороны всех групп клиентов падает (в том числе, с точки зрения безопасности данных)

• использование всех цифровых каналов связи с банками растёт (включая мессенджеры, соцсети, видеозвонки и т.п.)

• основной критерий при выборе финансовых услуг соотношение цена/качество

• пионеры и прагматики чаще других обращаются за финансовыми рекомендациями

• пионеры и прагматики до 44 лет переходят в необанки ради новых услуг

Вспомним и о директиве PSD2, обязывающей европейские банки открыть безопасный доступ к данным своих клиентов для провайдеров финансовых услуг. Американское Бюро по финансовой защите потребителей тоже задумывается о похожей инициативе учитывая размах общего тренда, скорее всего, это лишь вопрос времени.

Таким образом, традиционные банки вынуждены открываться внешним партнёрам, превращаясь по сути в платформу для сторонних разработок. Причём давление происходит одновременно с двух сторон: со стороны спроса и со стороны регуляторов.

А раз эффективной альтернативы нет (сказал он, категорично сдвинув брови), состоявшимся банковским брендам лучше искать в сложившейся ситуации возможности, а не уязвимости.

Банкам приходится впускать финтех в свою жизнь. Однако, это вовсе не вторжение Uber в мир таксопарков. Обе стороны нуждаются друг в друге, чтобы оставаться на рынке. И вот почему...

У банков есть проблемы

1. Безопасность

У банков есть серьёзные проблемы с безопасностью. К такому наблюдению (помимо меня) пришла ImmuniWeb международная компания, занимающаяся безопасностью приложений. Анализ сайтов и приложений 100 банков с разных континентов показал, что:

• 85% веб-приложений не соответствует GDPR (General Data Protection Regulation)

• 49% веб-приложений не соответствует PCI DSS (Payment Card Industry Data Security Standard)

• 92% приложений имеет хотя бы одну уязвимость среднего уровня риска

• 100% сайтов имеет уязвимости, связанные с забытыми субдоменами (читай сами создают комфортные условия для фишинга)

А это только отчёт о проверке крупных банков. Можно даже не сомневаться, что у банков поменьше недостаточно ресурсов, чтобы обеспечить даже такой уровень безопасности.

Чтобы сделать безопасный IT-продукт, в первую очередь нужны толковые люди с соответствующими хард скилами.

Вот только IT-специалистов сейчас не найти. Душная корпоративная культура не позволяет конкурировать с более привлекательными для талантливой молодёжи стартапами (например, финтех-стартапами) и гигантами вроде Google, Apple и Facebook. А опытные айтишники хотят чудовищно много денег.

2. Технологичность

Банки тратят уйму ресурсов на внедрение новой технологии. Ведь нужно, чтобы она как-то встраивалась в бизнес-процессы и дружила с легаси кодом, которые накопились в огромных масштабах. Неподъёмных для полноценного осмысления одним человеком, а иногда даже небольшой командой.

Финтех-компании лишены подобного наследия, они более гибки при выборе пути и внедрении инноваций. К тому же, партнёрство с ними дешевле и внутренних разработок, и покупки целых стартапов.

3. Зарегулированность

Деятельность кредитных организаций (aka банков) очень сложная система законов, правил и норм, которым нужно соответствовать.

Требования сокращают желание рисковать и ухудшают возможности для инноваций. Драйвером развития в таких условиях может стать кто-то, чьи решения не ограничивает банковская лицензия.

К финтех-компаниям без банковской лицензии требования тоже есть, но они гораздо менее суровые.

И у финтеха тоже есть проблемы

1. Доверие

Каким бы модным и классным не выглядел финтех, ему хватает своих проблем. Одна из них доверие клиентов. Всё-таки традиционалисты и скептики преобладают на большинстве рынков, включая европейские и российский.

По данным Accenture, 20% потребителей не готовы доверить свое финансовое благополучие необанкам. А 45% уверены, что подобные стартапы не продержатся на рынке и года! (отчасти это правда, что уж)

2. Прибыльность

Несмотря на взрывные темпы прироста клиентов, необанки далеко не всегда способны зарабатывать. Так, британские звезды Monzo и Starling столкнулись с удвоением убытков в 2020 г.

Один из первых необанков в мире Moven закрылся по этой причине. Самое показательное в кейсе то, что с рынка он не ушёл трансформировался в IT-компанию для обслуживания кредитных организаций. И сразу же заключил контракт с STC Pay, цифровым кошельком из Саудовской Аравии. Ещё раз подтверждая, что финтех эффективнее в союзе с банком.

Starling в итоге опять вышел в плюс. Но общее правило сохраняется: стабильная доходность не про инновационные финтехи. Даже сами представители брендов не упоминают рентабельность среди главных приоритетов.

3. Отсутствие лицензии

Отсутствие банковской лицензии дар и проклятие финтеха (той его части, у которой нет лицензии). С одной стороны, без жёсткого контроля проще использовать новые технологии и разрабатывать новые типы услуг.

С другой, без разрешения на ведение банковской деятельности большинство из этих услуг просто не получается оказывать. Финтех (в основном), как ни крути, это придаток банка, а не полноценный банк.

А вместе у них есть успешные коллабы

Итак, традиционным банкам и финансовым IT-компаниям необходима правильная модель кооперации. Модель, при которой банки отвечают за операции со счетами и предоставляют партнёрам свой бренд, авторитет и доступ к клиентской базе.

IT-шники, в свою очередь, разрабатывают технологии для сформированного спроса, получают необходимые ресурсы и инвестиции, сокращая свои риски. А заодно экономя бюджет и нервы больших братьев. Подобные схемы партнёрства уже встречаются на рынке.

1. Westpac, второй по капитализации банк Австралии, сотрудничает с британским провайдером облачных технологий 10x Future Technologies. Вместе они создают BaaS-платформу (Banking as a Service), которая позволит другим брендам предлагать услуги банка своим клиентам.

2. В России Модульбанк выступил платформой для финтех-сервиса Хайс, а Открытие и QIWI платформой для Точки. Работа на платформе кредитной организации позволила Точке и Хайсу тоже называть себя банками, хотя фактически они ими не являются.

3. Немецкий Commerzbank работает с одной из ведущих европейских систем удалённой идентификации и электронных подписей IDnow. Она позволяет клиентам подключаться к банку откуда угодно и когда угодно с удобным способом верификации например, по видео.

4. Ак Барс Банк скооперировался с технологической платформой APIBank для выпуска нового типа карт LetyBank.

5. Бывает и наоборот: в США финтех LendingClub купил себе бостонский Radius Bank. Основная причина сделки получение доступа к более дешевому и удобному фондированию.

Банки потеряют своих клиентов

Открытость для внешних партнёров означает, что доступ к данным клиентов больше не будет прерогативой банка. И это главное, что останавливает корпорации от новых коллабораций (корпорации коллаборировались, коллаборировались, да не выколлаборировались)

Статистика PriceWaterhouseCoopers показывает, что 88% игроков рынка боится потерять деньги из-за партнёрств с финтехом (по сравнению с 83% пять лет назад). Они полагают, что рискуют, в среднем, 24% выручки.

С другой стороны, это также плохо для банка, как хорошо для клиента.

Если тренд будет усиливаться, традиционным банкам рано или поздно придётся отказаться от стратегий развития, основанных на эксклюзивности данных (допродажи, агрессивное кредитование и т.п.) и нужно будет искать другие точки роста. Финтехам было бы полезно тоже задуматься об этом, чтобы помочь потенциальным партнёрам.

Я бы сказал, что роль банка, как процессингового центра для ИТ-компаний сейчас выглядит как вполне реалистичная модель. С ней обе стороны могут заниматься тем, что у них получается лучше всего.

А это конец статьи.

Автор идеи и главный редактор Денис Элиановский. Статью написал Станислав Лушин. Редактор Татьяна Китаева

Проблемы, страхи и надежды традиционных банков попробовал резюмировать Google. В своём предложении стать партнёрами, которое составил для них в прошлом году:

"Youre lagging in technology. Your current vendors are years behind. Consumers think youre irrelevant. Were hip, were cool, we have all the latest technologies, and boy have weve got data! Come partner with us on our new checking account!"

Перевод: Вы опаздываете с технологиями. Ваши вендоры отстали на годы. Клиенты считают вас устаревшими. Мы модные, мы классные, у нас все последние технологии и боже мой, сколько же у нас данных! Становитесь нашими партнёрами в новом банковском приложении!

Забавно, что 75% первых партнёров, похоже, даже не знали, о каком приложении речь! Это был Google Plex часть обновленного Google Pay с фичами вроде кэшбэков, специальных предложений и персональных финансовых советов.

Пример ярко показывает кризис привычной модели взаимодействия банк-клиент. Обе стороны ощущают нехватку в этой цепочке ещё одного звена. И это финтех.

Я уже несколько лет работаю с банками, ниже моё субъективное представление о важности и возможностях их коллаборации с ИТ-шниками, как со внешними независимыми разработчиками + статистика и ссылки на источники.

Мы поговорим о том кейсе, когда данные клиента передаются не с целью их продажи (чего порядочные банки и так не делают), а с целью улучшения клиентского сервиса. За что банки иногда даже сами приплачивают.

Банки вынуждены открывать доступ для финтеха

Как говорится из каждого утюга после пандемии мир никогда не будет прежним. А вот в нашем случае вирус выступил катализатором уже происходящих процессов.

Ожидания потребителей по отношению к банкам изменились ещё в доковидную эпоху. Прошлогоднее исследование Ernst&Young показало, что клиенты ждут от банкинга:

1. Прозрачности и возможности доверять

2. Максимальной персонализации

3. Защиты данных

4. Мгновенного удовлетворения запросов

5. Помощи с достижением личных целей

Иными словами, того, что умеют делать agile финтехи и не умеют мастодонты старой закалки.

Accenture, в свою очередь, проанализировали изменение поведения людей в связи с локдаунами. Тенденции подтвердились. Клиентов банков поделили на четыре группы пользователей:

1. Прагматики ценят полезность услуг вне зависимости от используемого канала

2. Традиционалисты предпочитают общаться с людьми, а не онлайн

3. Пионеры любители инноваций, первыми пробуют новинки

4. Скептики разбираются в технологиях, но часто недовольны качеством сервисов

И вот несколько выводов:

• доверие к банкам со стороны всех групп клиентов падает (в том числе, с точки зрения безопасности данных)

• использование всех цифровых каналов связи с банками растёт (включая мессенджеры, соцсети, видеозвонки и т.п.)

• основной критерий при выборе финансовых услуг соотношение цена/качество

• пионеры и прагматики чаще других обращаются за финансовыми рекомендациями

• пионеры и прагматики до 44 лет переходят в необанки ради новых услуг

Вспомним (или только что узнаем) и о директиве PSD2, обязывающей европейские банки открыть безопасный доступ к данным своих клиентов для провайдеров финансовых услуг.

Американское Бюро по финансовой защите потребителей тоже задумывается о похожей инициативе. Учитывая размах общего тренда, скорее всего, это лишь вопрос времени.

Таким образом, традиционные банки вынуждены открываться внешним партнёрам, превращаясь по сути в платформу для сторонних разработок. Причём давление происходит одновременно с двух сторон: со стороны спроса и со стороны регуляторов.

И раз уж такая ситуация сложилась, банкам эффективнее искать в ней возможности, а не уязвимости.

Приходится впускать финтех в свою жизнь. Однако, это совсем не вторжение Uber в мир таксопарков. Обе стороны нуждаются друг в друге, чтобы оставаться на рынке. И вот почему...

У банков есть проблемы

1. Безопасность

У банков есть серьёзные проблемы с безопасностью. К такому наблюдению пришла ImmuniWeb международная компания, занимающаяся безопасностью приложений. Анализ сайтов и приложений 100 банков с разных континентов показал, что:

• 85% веб-приложений не соответствует GDPR (General Data Protection Regulation)

• 49% веб-приложений не соответствует PCI DSS (Payment Card Industry Data Security Standard)

• 92% приложений имеет хотя бы одну уязвимость среднего уровня риска

• 100% сайтов имеет уязвимости, связанные с забытыми субдоменами (читай сами создают комфортные условия для фишинга)

Это отчёт о проверке крупных банков. Можно даже не сомневаться, что у банков поменьше недостаточно ресурсов, чтобы обеспечить даже такой уровень безопасности.

Чтобы сделать безопасный IT-продукт, в первую очередь нужны толковые люди с соответствующими хард скилами.

Вот только IT-специалистов сейчас не найти. Сложившаяся корпоративная культура не позволяет конкурировать с более привлекательными для талантливой молодёжи стартапами (например, финтех-стартапами) и гигантами вроде Google, Apple и Facebook. А опытные айтишники хотят чудовищно много денег.

2. Технологичность

Банки тратят уйму ресурсов на внедрение новой технологии. Ведь нужно, чтобы она как-то встраивалась в бизнес-процессы и дружила с легаси кодом, которые накопились в огромных масштабах. Неподъёмных для полноценного осмысления одним человеком, а иногда даже небольшой командой.

Финтех-компании лишены подобного наследия, они более гибки при выборе пути и внедрении инноваций. К тому же, партнёрство с ними дешевле и внутренних разработок, и покупки целых стартапов.

3. Зарегулированность

Деятельность кредитных организаций (aka банков) очень сложная система законов, правил и норм, которым нужно соответствовать.

Требования сокращают желание рисковать и ухудшают возможности для инноваций. Драйвером развития в таких условиях может стать кто-то, чьи решения не ограничивает банковская лицензия.

К финтех-компаниям без банковской лицензии требования тоже есть, но они гораздо менее суровые.

И у финтеха тоже есть проблемы

1. Доверие

Каким бы модным и классным не выглядел финтех, ему хватает своих проблем. Одна из них доверие клиентов. Всё-таки традиционалисты и скептики преобладают на большинстве рынков, включая европейские и российский.

По данным Accenture, 20% потребителей не готовы доверить свое финансовое благополучие необанкам. А 45% уверены, что подобные стартапы не продержатся на рынке и года! (отчасти это правда, что уж)

2. Прибыльность

Несмотря на взрывные темпы прироста клиентов, необанки далеко не всегда способны зарабатывать. Так, британские звезды Monzo и Starling столкнулись с удвоением убытков в 2020 г.

Один из первых необанков в мире Moven закрылся по этой причине. Самое показательное в кейсе то, что с рынка он не ушёл трансформировался в IT-компанию для обслуживания кредитных организаций. И сразу же заключил контракт с STC Pay, цифровым кошельком из Саудовской Аравии. Ещё раз подтверждая, что финтех эффективнее в союзе с банком.

Starling в итоге опять вышел в плюс. Но общее правило сохраняется: стабильная доходность не про инновационные финтехи. Даже сами представители брендов не упоминают рентабельность среди главных приоритетов.

3. Отсутствие лицензии

Отсутствие банковской лицензии дар и проклятие финтеха (той его части, у которой нет лицензии). С одной стороны, без жёсткого контроля проще использовать новые технологии и разрабатывать новые типы услуг.

С другой, без разрешения на ведение банковской деятельности большинство из этих услуг просто не получается оказывать. Финтех (в основном), как ни крути, это придаток банка, а не полноценный банк.

А вместе у них есть успешные коллабы

Итак, традиционным банкам и финансовым IT-компаниям необходима правильная модель кооперации. Модель, при которой банки отвечают за операции со счетами и предоставляют партнёрам свой бренд, авторитет и доступ к клиентской базе.

IT-шники, в свою очередь, разрабатывают технологии для сформированного спроса, получают необходимые ресурсы и инвестиции, сокращая свои риски. А заодно экономя бюджет и нервы больших братьев. Подобные схемы партнёрства уже встречаются на рынке.

1. Westpac, второй по капитализации банк Австралии, сотрудничает с британским провайдером облачных технологий 10x Future Technologies. Вместе они создают BaaS-платформу (Banking as a Service), которая позволит другим брендам предлагать услуги банка своим клиентам.

2. В России Модульбанк выступил платформой для финтех-сервиса Хайс, а Открытие и QIWI платформой для Точки. Работа на платформе кредитной организации позволила Точке и Хайсу тоже называть себя банками, хотя фактически они ими не являются.

3. Немецкий Commerzbank работает с одной из ведущих европейских систем удалённой идентификации и электронных подписей IDnow. Она позволяет клиентам подключаться к банку откуда угодно и когда угодно с удобным способом верификации например, по видео.

4. Ак Барс Банк скооперировался с технологической платформой APIBank для выпуска нового типа карт LetyBank.

5. Бывает и наоборот: в США финтех LendingClub купил себе бостонский Radius Bank. Основная причина сделки получение доступа к более дешевому и удобному фондированию.

Банки потеряют своих клиентов?

Открытость для внешних партнёров означает, что доступ к данным клиентов больше не будет прерогативой банка. И это главное, что останавливает корпорации от новых коллабораций (корпорации коллаборировались, коллаборировались, да не выколлаборировались)

Статистика PriceWaterhouseCoopers показывает, что 88% игроков рынка боится потерять деньги из-за партнёрств с финтехом (по сравнению с 83% пять лет назад). Они полагают, что рискуют, в среднем, 24% выручки.

С другой стороны, это также плохо для банка, как хорошо для клиента.

Если тренд будет усиливаться, традиционным банкам рано или поздно придётся отказаться от стратегий развития, основанных на эксклюзивности данных (допродажи, агрессивное кредитование и т.п.) и нужно будет искать другие точки роста. Финтехам было бы полезно тоже задуматься об этом, чтобы помочь потенциальным партнёрам.

Я бы сказал, что роль банка, как процессингового центра для ИТ-компаний сейчас выглядит как вполне реалистичная модель. С ней обе стороны могут заниматься тем, что у них получается лучше всего.

Но для этого придётся шэрить данные, и в гораздо больших объёмах, чем это делается сейчас.

Было бы здорово послушать ваше мнение, стоит ли банкам открывать часть данных внешним компаниям в обмен на лучший сервис для своих клиентов?

Автор идеи и главный редактор Денис Элиановский. Статью написал Станислав Лушин. Редактор Татьяна Китаева

Гистограммам не чужды систематические ошибки. Дело в том, что они достаточно условны и могут привести к неправильным выводам о данных. Если вы хотите визуализировать переменную, лучше выбрать другой график.

Независимо от того, находитесь ли вы на встрече с высшим руководящим составом или со специалистами по обработке данных, в одном вы можете быть уверены: в какой-то момент появится гистограмма.

И нетрудно догадаться почему. Гистограммы весьма интуитивно наглядны: любой поймет их с первого взгляда. Более того, они объективно представляют реальность, не так ли? А вот и нет.

Гистограмма может ввести в заблуждение и привести к ошибочным выводам даже на простейшем наборе данных!

В этой статье мы на примерах рассмотрим 6 причин, почему, когда дело доходит до визуализации данных, гистограммы точно не является лучшим выбором:

1. Они слишком сильно зависят от количества интервалов.

2. Они слишком сильно зависят от максимума и минимума переменной.

3. Они не дают возможности заметить значимые значения переменной.

4. Они не позволяют отличить непрерывные переменные от дискретных.

5. Они делают сравнение распределений сложным.

6. Их построение затруднено, если в памяти находятся не все данные.

Ладно, я понял: гистограммы не идеальны. Но есть ли у меня выбор? Конечно есть!

В конце статьи я порекомендую другой график, называемый CDP, который минует эти недостатки.

Итак, что же не так с гистограммой?

1. Она слишком сильно зависит от количества интервалов.

Чтобы построить гистограмму, вы должны сначала определить количество интервалов, также называемых корзинами (bins). Для этого существует множество различных практических методов (вы можете ознакомиться с их обзором на этой странице). Но насколько критичен этот выбор? Давайте возьмем реальные данные и посмотрим, как меняется гистограмма в зависимости от количества интервалов.

Переменная представляет собой максимальную частоту сердечных сокращений (ударов в минуту), полученную у 303 людей во время некоторой физической активности (данные взяты из набора данных UCI по сердечным заболеваниям: источник).

Глядя на верхний левый график (который мы получим по умолчанию в Python и R), у нас сложится впечатление хорошего распределения с одним пиком (модой). Однако если бы мы рассмотрели бы другие варианты гистограммы, мы получили бы совершенно другую картину. Разные гистограммы одних и тех же данных могут привести к противоречивым выводам.

2. Она слишком сильно зависит от максимума и минимума переменной.

Даже после того, как количество интервалов установлено, интервалы зависят от положения минимума и максимума переменной. Достаточно, чтобы один из них немного изменился, и все интервалы также изменятся. Другими словами, гистограммы не являются надежными.

Например, давайте попробуем изменить максимум переменной, не меняя количество интервалов.

Отличается только одно значение, а весь график получается другим. Это нежелательное свойство, потому что нас интересует общее распределение: одно значение не должно так влиять на график!

3. Не дает возможности заметить значимые значения переменной.

Если в общем, то когда переменная содержит некоторые часто повторяющиеся значения, нам конечно нужно об этом знать. Однако гистограммы этому препятствуют, поскольку они основаны на интервалах, а интервалы скрывают отдельные значения.

Классическим примером является случай, когда отсутствующим значениям массово присваивается 0. В качестве примера давайте рассмотрим набор данных переменной, состоящий из 10 тысяч значений, 26% из которых нули.

График слева это то, что вы получаете по умолчанию в Python. Глядя на него, вы не заметите скопление нулей, и вы даже можете подумать, что эта переменная имеет плавную динамику.

График справа получен путем сужения интервалов и дает более четкое представление о реальности. Но дело в том, что как бы вы ни сужали интервалы, вы никогда не будете уверены, содержит ли первый интервал только 0 или какие-то другие значения.

4. Не позволяет отличить непрерывные переменные от дискретных.

Зачастую мы бы хотели знать, является ли числовая переменная непрерывной или дискретной. По гистограмме это практически невозможно сказать.

Возьмем переменную Возраст (Age). Вы можете получить Возраст = 49 лет (когда возраст округлен) или Возраст = 49,828884325804246 лет (когда возраст рассчитывается как количество дней с момента рождения, деленное на 365,25). Первая дискретная переменная, вторая непрерывная.

Тот, что слева, непрерывен, а тот, что справа, дискретен. Однако на верхних графиках (по умолчанию в Python) вы не увидите никакой разницы между ними: они выглядят совершенно одинаково.

5. Сложно сравнивать распределения.

Часто бывает необходимо сравнить одну и ту же переменную в разных кластерах. Например, в отношении данных UCI о сердечных заболеваниях, приведенных выше, мы можем сравнить:

• все население (для справки)

• люди моложе 50 страдающие сердечными заболеваниями

• люди моложе 50 НЕ страдающие сердечными заболеваниями

• люди старше 60 лет страдающие сердечными заболеваниями

• люди старше 60 и НЕ страдающие сердечными заболеваниями.

Вот что мы получили бы в итоге:

Гистограммы основаны на областях, и, когда мы пытаемся провести сравнение, области в конечном итоге перекрываются, что делает эту задачу практически невыполнимой.

6. Сложно построить, если в памяти находятся не все данные.

Если все ваши данные находятся в Excel, R или Python, построить гистограмму легко: в Excel вам просто нужно кликнуть по иконке гистограммы, в R выполнить команду hist(x), а в Python plt.hist(х).

Но предположим, что ваши данные хранятся в базе данных. Вы же не хотите выгружать все данные только для того, чтобы построить гистограмму, верно? По сути, все, что вам нужно, это таблица, содержащая для каждого интервала крайние значения и количество наблюдений. Примерно такая:

| INTERVAL_LEFT | INTERVAL_RIGHT | COUNT |

|---------------|----------------|---------------|

| 75.0 | 87.0 | 31 |

| 87.0 | 99.0 | 52 |

| 99.0 | 111.0 | 76 |

| ... | ... | ... |

Но получить ее с помощью SQL-запроса не так просто, как кажется. Например, в Google Big Query код будет выглядеть так:

WITHSTATS AS (  SELECT     COUNT(*) AS N,    APPROX_QUANTILES(VARIABLE_NAME, 4) AS QUARTILES  FROM    TABLE_NAME),BIN_WIDTH AS (  SELECT    -- freedman-diaconis formula for calculating the bin width    (QUARTILES[OFFSET(4)]  QUARTILES[OFFSET(0)]) / ROUND((QUARTILES[OFFSET(4)]  QUARTILES[OFFSET(0)]) / (2 * (QUARTILES[OFFSET(3)]  QUARTILES[OFFSET(1)]) / POW(N, 1/3)) + .5) AS FD  FROM     STATS),HIST AS (  SELECT     FLOOR((TABLE_NAME.VARIABLE_NAME  STATS.QUARTILES[OFFSET(0)]) / BIN_WIDTH.FD) AS INTERVAL_ID,    COUNT(*) AS COUNT  FROM     TABLE_NAME,    STATS,    BIN_WIDTH  GROUP BY     1)SELECT   STATS.QUARTILES[OFFSET(0)] + BIN_WIDTH.FD * HIST.INTERVAL_ID AS INTERVAL_LEFT,  STATS.QUARTILES[OFFSET(0)] + BIN_WIDTH.FD * (HIST.INTERVAL_ID + 1) AS INTERVAL_RIGHT,  HIST.COUNTFROM   HIST,   STATS,   BIN_WIDTH

Немного громоздко, не правда ли?

Альтернатива: график кумулятивного распределения.

Узнав 6 причин, по которым гистограмма не является идеальным выбором, возникает естественный вопрос: Есть ли у меня альтернатива? Хорошие новости: существует лучшая альтернатива, которая называется График кумулятивного распределения (Cumulative Distribution Plot - CDP). Я знаю, что это название не такое запоминающееся, но гарантирую, оно того стоит.

График кумулятивного распределения это график квантилей переменной. Другими словами, каждая точка CDP показывает:

• по оси x: исходное значение переменной (как в гистограмме);

• по оси y: сколько наблюдений имеют такое же или меньшее значение.

Давайте посмотрим на пример с переменной максимальной частотой пульса.

Возьмем точку с координатами x = 140 и y = 90 (30%). По горизонтальной оси вы видите значение переменной: 140 ударов сердца в минуту. По вертикальной оси вы видите количество наблюдений, у которых частота сердцебиение равна или ниже 140 (в данном случае 90 человек, что означает 30% выборки). Следовательно, у 30% нашей выборки максимальная частота сердцебиения составляет 140 или менее ударов в минуту.

Какой смысл в графике, показывающем, сколько наблюдений равно или ниже заданного уровня? Почему не просто равно? Потому что в противном случае результат зависел бы от отдельных значений переменной. И это не сработает, потому что каждое значение имеет очень мало наблюдений (обычно только одно, если переменная непрерывна). Напротив, CDP полагаются на квантили, которые более стабильны, выразительны и легко читаются.

Вдобавок CDP намного полезнее. Если задуматься, вам часто приходится отвечать на такие вопросы, как у скольких из них от 140 до 160? Или у скольких из них больше 180?. Имея перед глазами CDP, вы можете дать немедленный ответ. С гистограммой это было бы невозможно.

CDP решает все проблемы, которые мы видели выше. Фактически, по сравнению с гистограммой:

1. Не требует пользовательского выбора. Для одного набора данных, существует только один возможный CDP.

2. Не страдает от выпадающих значений. Экстремальные значения не влияют на CDP, поскольку квантили не меняются.

3. Позволяет определять значимые значения. Если существует концентрация точек данных на каком-то конкретном значении, это сразу видно, поскольку будет вертикальный сегмент, соответствующий значению.

4. Позволяет с первого взгляда распознать дискретную переменную. Если существует только конкретный набор возможных значений (т.е. переменная дискретна), это сразу видно, поскольку кривая примет форму лестницы.

5. Упрощает сравнение распределений. На одном графике легко сравнить два или более распределения, поскольку это просто кривые, а не области. Кроме того, ось y всегда находится в диапазоне от 0 до 100%, что делает сравнение еще более простым. Для сравнения, это пример, который мы видели выше:

6. Его легко построить, даже если у вас нет всех данных в памяти. Все, что вам нужно, это квантили, которые можно легко получить с помощью SQL:

SELECT   COUNT(*) AS N,  APPROX_QUANTILES(VARIABLE_NAME, 100) AS PERCENTILESFROM  TABLE_NAME

Как построить график кумулятивного распределения в Excel, R, Python

В Excel вам нужно построить два столбца. Первый с 101 числом, равномерно распределенными от 0 до 1. Второй столбец должен содержать процентили, которые могут быть получены по формуле: =PERCENTILE(DATA, FRAC), где DATA - это вектор, содержащий данные, а FRAC - это первый столбец: 0,00, 0,01, 0,02, 0,03,, 0,98, 0,99, 1. Затем вам просто нужно построить график по этим двум столбцам, разместив значения переменной на оси x.

В R это делается в одну строчку:

plot(ecdf(data))

В Python:

from statsmodels.distributions.empirical_distribution import ECDFimport matplotlib.pyplot as pltecdf = ECDF(data)plt.plot(ecdf.x, ecdf.y)

Спасибо за внимание! Надеюсь, эта статья оказалась для вас полезной.

Я ценю отзывы и конструктивную критику. Если вы хотите поговорить об этой статье или других связанных темах, вы можете написать мне в Linkedin.

Перевод материала подготовлен в рамках онлайн-курса "Machine Learning. Basic". Всех заинтересованных приглашаем на день открытых дверей курса, где можно будет узнать все подробности об обучении и пообщаться с преподавателем.

- Узнать подробнее о курсе "Machine Learning. Basic"

- Смотреть онлайн-встречу "День открытых дверей"

Зарабатывать продажей лекарств, которые заведомо не работают, не только аморально, но и не особо легко. Люди всё-таки обычно не хотят покупать препараты, неэффективность которых была доказана. А вот если вы сумели выдавить заветное p < 0.05 в пользу того, что акупунктура таки работает из данных, которые явно утверждают обратное, то серия публикаций, успех в карьере и вечная благодарность всех акупунктурщиков вам гарантированы.

Со ставками на спорт история такая же: чтобы выигрывать деньги у букмекера, нужно уметь считать коэффициенты лучше него. А вот чтобы заработать на продаже советов, на кого ставить, достаточно считать лучше своих читателей. Например, Винс Экинс давал своим читателям следующий совет:

Philadelphia Eagles выиграли 10 предыдущих игр, если они не были фаворитом, играли на собственном поле и в предыдущем матче прошли более 150 ярдов в наступлении.

10 игр! Это подразумевает, что вероятность поражения в аналогичных условиях 1/(2^10) = 1/1024 = 0.0009. Такая степень уверенности не каждый день встречается даже в серьёзных исследованиях. И это, разумеется, не потому что Eagles действительно так уж хороши, а потому что анализ был проведён неправильно. Иногда так получается просто потому, что авторы не умеют нормально анализировать данные, и в результате целые области науки оказываются под вопросом. А иногда так поступают откровенные мошенники (что в академической науке, что за её пределами).

Как получить p<0.05

По определению, p-значение в 0.05 это вероятность в 5%, что результат будет получен, если верна нулевая гипотеза, то есть если лекарство не работает, исследуемые феномены не взаимосвязаны и вообще ничего интересного тут нет. Так что основная идея проста: если взять два десятка гипотез, то, скорее всего, одна из них будет иметь p0.05.

Получить много гипотез для проверки можно двумя основными способами. Если изначально разбить рассматриваемую выборку на множество подгрупп, то, скорее всего, хоть в одной из них да получится статистически значимый результат. Этот метод прекрасно иллюстрирует xkcd: мармеладки вообще прыщи не вызывают (p>0.05), и красные мармеладки, в частности, тоже не вызывают (p>0.05). И жёлтые не вызывают, и сиреневые, и оранжевые, и коричневые, и ещё два десятка цветов не вызывают а вот для зелёных p<0.05.

Второй метод (он же метод Латиноамериканской Бабушки) подразумевает деление исходной выборки на любые произвольные подгруппы до тех пор, пока не найдётся комбинация условий, при которой p<0.05. Допустим, в вымышленной стране существует вымышленное заболевание, от которого сама собой излечивается ровно половина больных. Вторая половина умирает. На первый взгляд, ваше чудо-лекарство, повышающее долю выживших аж до 50%, выглядит так себе даже по меркам British Journal of General Practice. Но всегда можно взглянуть поподробнее.

Допустим, по чистой случайности удачно вылечившиеся распределены по полам слегка неравномерно: выздоровело 49% мужчин и 51% женщин. А среди женщин старше 60 лет препарат помог аж 55%. И, допустим, 13 пожилых женщин, участвовавших в исследовании, родом из Мексики. Вполне может оказаться, что препарат помог 10 из них. Это уже не 50%, а 77%, и к тому же вполне приличное p-значение в 0.046. Можно наслаждаться репутацией спасителя латиноамериканских бабушек (а если бы с ними не прокатило можно было бы проверить еврейских мальчиков, девочек-негритянок, белых среднего возраста и все остальные комбинации).

Надеюсь, всем и так очевидно, что 10 человек из 13 это несерьёзно. Но тем не менее я повторю: если выборка слишком маленькая, то любой полученный на ней результат почти наверняка ничего не стоит.

Но что, если бабушек было не 13, а 90, и препарат помог аж 61? Даже если всем остальным он не помог, две трети выздоровевших и p-значение в 0.0005 выглядят впечатляюще, а 90 человек это уже вполне приличная выборка. Давайте посчитаем, но для начала немного теории.

P-значения в общем-то довольно плохой инструмент. Сама идея обнаруживать то, чего нет, не чаще, чем в 1 эксперименте из 20 звучит не очень впечатляюще, а к тому же даже в теории 5%-ный порог p-значения обманывает экспериментатора в 30% случаев. Но вот что они делают хорошо так это конвертируют любое распределение в равномерное. Например, если взять несколько значений из нормального распределения, то они в основном лягут примерно по центру. А вот их p-значения равномерно распределятся между 0 и 1.

В статье вы чаще всего увидите p-значения только для лучшей гипотезы (что прыщи возникают от зелёных мармеладок, а лекарство прекрасно помогает пожилым женщинам латиноамериканского происхождения). Первый шаг к нормальному анализу это поправка Бонферрони:

Правило Бонферрони: Порог P-значения в для одной гипотезы эквивалентен порогу в /N для лучшей из N гипотез.

Обычно её интерпретируют как верхнюю границу: порог в /N для N гипотез заставляет принять нулевую гипотезу не чаще, чем порог в для единственной протестированной гипотезы. Но, на самом деле, это неплохая аппроксимация: пусть h1,,hN это N p-значений для N независимых нулевых гипотез, и все они находятся в диапазоне от 0 до 1. Тогда вероятность того, что хотя бы одно из них ниже /N = P(min(h1,,hN) < /N) = 1 (1 /N)^N 1 e^- 1 (1-) = . Последние шаги основываются на линейной аппроксимации e^x 1+x, которая работает при близких к нулю x.

Но это доказательство работает для независимых гипотез. В случае мармеладок так оно и есть, а вот для мексиканских бабушек гипотезы взаимосвязаны. Шанс выздоровления, например, белых мужчин положительно коррелирует и с шансом выздоровления белых людей вообще (т.е. более широкой категории), и с шансом выздоровления белых пожилых мужчин (подмножества белых мужчин). Чтобы проверить применимость правила Бонферрони в такой ситуации, я прогнал симуляцию (код доступен на Github).

Пусть в нашей стране живут люди трёх возрастов (молодые, среднего возраста и старые), двух полов и четырёх рас. В каждой из 2*3*4=24 подгрупп по пятьсот человек, общее население 12 000. По условиям задачи болезнь убивает 50% больных, так что в среднем ожидается 12 000/2=6 000 выживших. Для всей выборки мы получим p=0.05, если выздоровеет 50.75% (90 дополнительных выздоровевших) и p=0.0005, если вылечится 51.5%.

Правило больших выборок: при достаточно большой выборке статистическая значимость не обязательно обозначает практическую значимость. Можно получить копеечный эффект с очень низким p-значением.

P-значения не особенно помогают добраться до истины. Но они повсюду, с ними легко работать, и они неплохо отсекают откровенную чепуху. Поскольку статья именно про это, я не буду влезать в байесовскую статистику и ограничусь обсуждением возможных манипуляций с p-значениями.

Вернёмся к нашей симуляции. Я прогнал её 1 000 раз для трёх лекарств: плацебо, вылечивающее 50%; статистически значимое лекарство, вылечивающее 50.75%; и хорошее лекарство, вылечивающее 51.5% (да, вот такие вот у нас критерии хорошего). Для каждого из лекарств я искал подгруппу, в которой оно выдаст лучшее p-значение:

13 hispanic 1    0.12253041651147314 female hispanic 2    0.18079730402678315 young hispanic 2    0.2517223358154316 young female hispanic 3    0.17187517 white 1    0.046230490536462118 female white 2    0.57223222404718419 young white 2    0.2517223358154320 young female white 3   0.945312521 adult 1    0.36877715449216222 female adult 2    0.78520474607830623 asian 1    0.95376950946353824 female asian 2    0.819202695973217

Второе число это глубина выбранной подгруппы (вся выборка 0, азиаты 1, азиатские женщины среднего возраста 3). В нашем случае возможно 60 групп: 1 полная выборка, 9 групп глубины 1, 26 глубины 2, 24 - глубины 3. Так что поправка Бонферрони требует порога p-значения в 0.05/60=0.00083

В каждой из 1000 симуляций я выбрал самую удачную группу и построил график. Цветами показана глубина подгруппы, вертикальные линии соответствуют нескорректированному значению в 0.05 и скорректированному 0.00083. По горизонтальной оси логарифм p-значения, по вертикальной сколько симуляций (из 1000) имеют значение не ниже данного.

Безо всякого мошенничества плацебо получает p<0.05 в 5% случаев (что очевидно из определения), значимое лекарство в 50% случаев, а хорошее в 95. Но если мы применим поправку Бонферрони, то работающие лекарства пройдут проверку всего в 23% и 72% случаев соответственно. У плацебо дела ещё хуже, но всё-таки получается, что в таких случаях поправка оказывается чересчур агрессивной.

Как заметить подвох

Что всё это даёт на практике? Давайте соберём воедино все имеющиеся у нас советы.

1. Мощность превыше всего: если результат получен на крошечной выборке (и особенно если в исходных данных много шума) дальше можно не читать.

2. Считайте категории: если вы видите результат только для одной подгруппы, посчитайте, сколько аналогичных категорий можно выделить из тех же данных.

3. Применяйте поправку: разделите исходный порог p-значения на то, что получилось в предыдущем пункте (или, что эквивалентно, умножьте на него само p-значение). Если полученный результат вас устраивает, хорошо.

4. Сохраняйте скептицизм: даже со всеми поправками сомнительный результат есть сомнительный результат, а статистическая значимость необязательно подразумевает реальную значимость. Сколь угодно низкое p-значение ещё ничего не гарантирует.

Существует мета-анализ, подтверждающий с p=0.00000000012 способность некоторых людей предсказывать будущее. Цифра потрясающая, но она имеет смысл, только если нет ни малейших сомнений в том, что это исследование (и все предыдущие работы, на которые оно опирается) было проведено безупречно. Если есть причины считать, что это не так, то на самом деле p-значение намного выше полученного.

Можно предположить, что как минимум одна психологическая статья из тысячи выполнена некорректно, а то и вовсе написана по сфабрикованным результатам. Соответственно, любое полученное в психологии p-значение ниже 1/1000 ничем не лучше p-значения ровно в 1e-3. Вероятность получить результат в отсутствие изучаемого феномена равна собственно p-значению плюс вероятность того, что исследование в целом некорректно.

Посмотрим, что этот метод говорит нам насчёт ставки на Eagles. Во-первых, раз речь идёт о 10 выигранных играх, то 11-ю игру (в смысле 11-ю с конца, предшествующую этим 10) в аналогичных условиях они проиграли, иначе обсуждалась бы серия из 11 побед. Во-вторых, 10 игр это не так чтоб очень большая выборка, но зато нет никакой погрешности измерения. Мы знаем со стопроцентной вероятностью, выиграли ли они ту или иную игру или проиграли. К мощности эксперимента формальных претензий нет.

А вот со вторым пунктом некоторые проблемы . Даже если Eagles действительно неплохо играют как не-фаворит на своём поле после удачного наступления в предыдущем матче, тот же самый матч может быть описан как Игра сиэтлской команды после победы на выезде против команды, которая плохо пасует (прим. пер.: речь об описании матча с точки зрения их противников, Seattle Seahawks), или как Матч команды восточного дивизиона на западе против команды, выигравшей предыдущий матч, или ещё несколькими тысячами способов. Число возможных параметров сложно даже посчитать, но попробуем прикинуть:

1. Описание одной команды: дивизион, родной город, история в этом сезоне, результат предыдущей игры, качество атаки и защиты, статистика отдельных игроков итого не меньше 20 параметров.

2. Столько же для другой команды.

3. Обстоятельства игры: домашний / выездной матч, погода, время года, состояние поля, история игр между конкретными двумя командами итого не меньше 10 параметров.

Даже если выбрать в каждой категории всего по 1 параметру, получится 4000 моделей. Что это означает для Eagles? Вероятность поражения, согласно данной модели, 1/1024, но поправка Бонферрони для 4000 гипотез говорит нам, что примерно 4 гипотезы аналогичной сложности должны оказаться верными по чистому совпадению. Разумеется, Экинс не перебрал их все; он просто порылся в данных, нашёл интересное совпадение и опубликовал его. Но при таком количестве возможных гипотез единственное совпадение ничего не стоит. В обсуждаемом матче Eagles проиграли со счётом 15:26.

Тот же критерий нужно применять каждый раз, когда кто-то хвастается превосходством в неожиданно узкой категории. На этой рекламе коллектив Тринити-Университета признан лучшим в Техасе (и 6-м в стране) среди преподавателей, занимающихся исследованиями и читающих лекции. При виде этой фразы сразу нужно задуматься, сколько ещё категорий было в цитируемом исследовании, сколько вообще существует рейтингов, и в скольких из них Тринити никаких мест не занял. Университет 1 в общем рейтинге, скорее всего, действительно хорош в преподавании; про университет 1 в данной узкой категории мы знаем только то, что они умеют перебирать все возможные рейтинги и подкатегории.

1. Введение

При выполнении инженерно-геологических изысканий может возникнуть задача, связанная с сопоставлением данных полевых и лабораторных исследований на одних и тех же грунтах, с целью подтверждения корректной транспортировки проб от объекта изысканий до лаборатории (образцы не были деформированы и/или разрушены в ходе перевозки).

При данной постановке задачи можно применить методику A/B-тестирования со следующими параметрами:

1. Измеряемой метрикой будет среднее значение плотности скелета грунта (p_d, г/см³), характеризующее сложение проб. Данная величина имеет нормальный закон распределения;

2. Критерием проверки гипотезы будет служить t-критерий (критерий Стьюдента):для двух независимых выборок, если сопоставляемые полевые (до транспортировки) и лабораторные (после транспортировки) данные проводились на разных пробах грунта;для двух зависимых выборок, если исследования выполнены на одних и тех же пробах.

В рамках данной темы мы сгенерируем две случайные выборки, которые будем сопоставлять, сформулируем статистические гипотезы, проверим их и сделаем выводы.

2. Генерация выборок

2.1 Оценка объема выборок

В рамках дизайна эксперимента, перед генерацией выборок плотностей, прикинем их необходимый объем при заданномразмере эффекта (ES - effect size),мощности (power)идопустимой ошибке I рода ()(определения данных терминов приведено ниже). Расчет произведем с привлечением пакетаstatsmodels.

Размер эффекта (стандартизированный) величина, характеризующая различие, которое мы хотим выявить, равная отношению разности средних значений по выборкам к взвешенному стандартному отклонению. В нашем случае:

Взвешенное стандартное отклонениеS_pooledдля выборок одинакового размера можно расcчитать по формуле:

Существует условная классификация размера эффекта (Cohen, 1988) ES = 0.2 - маленький; 0.5 - средний; 0.8 - большой.

Мощность вероятность не совершить ошибку II рода (обычно принимается равной 80%).

Пояснения по ошибкам I и II рода приведены в таблице ниже:

Для описанных выше величин примем следующие значения:

• = 0.05 (вероятность выявить различия между средними при их отсутствии)

• ES = 0.5 (размер эффекта составит половину от дисперсии измеряемых величин плотности).

• Power = 0.8 (вероятность выявления установленного различия между средними значениями).

Теперь к коду:

#Импорт библиотекimport numpy as npfrom statsmodels.stats.power import TTestIndPowerfrom matplotlib.pyplot import figureimport matplotlib.pyplot as pltimport scipyfrom statsmodels.stats.weightstats import *

#Задаем параметрыeffect = 0.5alpha = 0.05power = 0.8analysis = TTestIndPower()#Оценка размера выборкиsize = analysis.solve_power(effect, power=power, alpha=alpha)print(f'Размер выборки, шт.: {int(size)}')

Размер выборки, шт.: 63

При заданных значения мощности, размера эффекта и уровня значимости минимальный размер каждой выборки должен составлять63пробы. Для красоты округлим полученное число до65шт.

Давайте построим график зависимости необходимого размера выборок от размера эффекта при заданной мощности и уровне значимости.

plt.figure(figsize=(10, 7), dpi=80)results = dict((i/10, analysis.solve_power(i/10, power=power, alpha=alpha))                for i in range(2, 16, 1))plt.plot(list(results.keys()), list(results.values()), 'bo-')plt.grid()plt.title('График зависимости необходимого объема выборки \n от размера эффекта')plt.ylabel('Размер выборки n, шт.')plt.xlabel('Размер эффекта ES, д.е.')for x,y in zip(list(results.keys()),list(results.values())):    label = "{:.0f}".format(y)    plt.annotate(label,                  (x,y),                  textcoords="offset points",                  xytext=(0,10),                  ha='center')plt.show()

Данный график позволяет увидеть, как быстро изменяется необходимый объем выборок при уменьшении фиксируемого размера эффекта ES. Например: при выявлении различия в плотности проб грунта до и после их транспортировки в 0,03г/см³при стандартном отклонении в 0,1г/cм³(ES = 0,03г/см³/ 0,1г/см³= 0,3 д.е.), необходимый объем проб по каждой выборке должен составить не менее 175 проб для заданной мощности и уровня значимости (power=0.80,=0.05).

2.2 Генерация выборок

Теперь зная необходимый минимальный размер выборок, сгенерируем их с помощью библиотекиnumpy.

Измеряемая физическая характеристика грунта (плотность скелета) имеет нормальный закон распределения. В рамках данного примера зададим генератору следующие значения среднего (X) и стандартного отклонения (S):

• для первой выборки X₁= 1,65г/см³,S₁= 0.15г/см³;

• для второй X₂= 1,60г/см³,S₂= 0.15г/см³.

loc_1 = 1.65sigma_1 = 0.15loc_2 = 1.60sigma_2 = 0.15sample_size = 65#Генерируем выборки с заданными параметрамиsample_1 = np.random.normal(loc=loc_1, scale=sigma_1, size=sample_size)sample_2 = np.random.normal(loc=loc_2, scale=sigma_2, size=sample_size)

Постоим гистограммы и "ящик с усами" по полученным выборкам.

fig, axes = plt.subplots(ncols=2, figsize=(18, 5))max_y = np.max(np.hstack([sample_1,sample_2]))#Гистрограмма по выборке 1count_1, bins_1, ignored_1 = axes[0].hist(sample_1, 10, density=True,                                           label="Выборка 1", edgecolor='black',                                          linewidth=1.2)axes[0].plot(bins_1, 1/(sigma_1 * np.sqrt(2 * np.pi)) *               np.exp( - (bins_1 - loc_1)2 / (2 * sigma_12)),         linewidth=2, color='r', label='плотность вероятности')axes[0].legend()axes[0].set_xlabel(u'Длина сессии, с')axes[0].set_ylabel(u'Количество сессий, шт.')axes[0].set_ylim([0, 5])axes[0].set_xlim([1.1, 2.2])#Гистрограмма по выборке 2count_2, bins_2, ignored_2 = axes[1].hist(sample_2, 10, density=True,                                           label="Выборка 2", edgecolor='black',                                           linewidth=1.2, color="green")axes[1].plot(bins_2, 1/(sigma_2 * np.sqrt(2 * np.pi)) *               np.exp( - (bins_2 - loc_2)2 / (2 * sigma_22)),         linewidth=2, color='r', label='плотность вероятности')axes[1].legend()axes[1].set_xlabel(u'Длина сессии, с')axes[1].set_ylabel(u'Количество сессий, шт.')axes[1].set_ylim([0, 5])axes[1].set_xlim([1.1, 2.2])plt.show()

#Ящик с усамиfig, ax = plt.subplots(figsize=(8, 8))axis = ax.boxplot([sample_1, sample_2], labels=['Выборка 1', 'Выборка 2'])data = np.array([sample_1, sample_2])means = np.mean(data, axis = 1)stds = np.std(data, axis = 1)for i, line in enumerate(axis['medians']):    x, y = line.get_xydata()[1]    text = ' ={:.2f}\n ={:.2f}'.format(means[i], stds[i])    ax.annotate(text, xy=(x, y))plt.ylabel('Плотность скелета грунта, г/см3')plt.show()

3. Формулировка гипотез

Пришло время для формулировки гипотез. У нас могут быть два случая:

• Случай 1. Сопоставляемые полевые и лабораторные данные по определению плотности скелета грунта относятся к разным пробам, тогда t-критерий будет рассчитываться для двух независимых выборок;

• Случай 2. Исследования в поле и лаборатории выполнены на одних и тех же пробах, тогда t-критерий будет рассчитываться для двух зависимых выборок.

Начнем с первого варианта.

Вариант 1. Для двух независимых выборок

С помощью двухвыборочного критерия Стьюдента проверим гипотезу о равенстве средних выборок.

Нулевая гипотезаH₀:средние значения равны1=2.

Альтернативная гипотезаH₁:средние не равны₁₂.

Статистика:

Нулевое распределение:T(X₁ⁿ¹,X₂ⁿ²)~St(), где степень свободывычисляется по следующей формуле

Для расчета достигаемого уровня значимости воспользуемся методомttest_indмодуляstats.

t_st, p_val = scipy.stats.ttest_ind(sample_1, sample_2, equal_var = False)print(f't-критерий составил {round(t_st, 2)}')print(f'Рассчитанный t-критерий дает достигаемый \уровень значимости (p-value) равный {round(p_val, 3)}')

t-критерий составил 2.92

Рассчитанный t-критерий дает достигаемый уровень значимости (p-value) равный 0.004

Вывод для варианта 1

Нулевая гипотезаH₀о том, что средняя плотность скелета грунта не изменилась после транспортировки,отвергаетсяна уровне значимости 0,05 (достигаемый уровень значимостиp-valueдля сгенерированных выборок составил 0.004) в пользу альтернативной.

Давайте интервально оценим разность средних по данным выборкам.

c_m = CompareMeans(DescrStatsW(sample_1), DescrStatsW(sample_2))print("95%% доверительный интервал: \[%.4f, %.4f]" % c_m.tconfint_diff(usevar='unequal'))

95% доверительный интервал: [0.0235, 0.1228]

Так как ноль не попадает в рассматриваемый 95% доверительный интервал, мы можем сделать вывод, что средние значения рассматриваемых выборок отличаются на уровне значимости в 5%.

Вариант 2. Для двух связанных выборок

Допустим, что оценка плотности скелета грунта в полевых (до транспортировки) и лабораторных (после транспортировки) условиях проводилась для каждого образца. Тем самым выборки будут является зависимыми, а проверка нулевой гипотезы об отсутствии изменений в плотности грунта при транспортировке будет осуществляться с помощью двухвыборочного критерия Стьюдента для связанных выборок.

Нулевая гипотезаH₀:средние значения равны ₁=₂.

Альтернативная гипотезаH₁:средние не равны₁₂.

Статистика:

Нулевое распределение: T(X₁ⁿ, X₂ⁿ)~St(n-1)

Для расчета достигаемого уровня значимости воспользуемся методомttest_relмодуляstats.

t_st, p_val = stats.ttest_rel(sample_1, sample_2)print(f't-критерий составил {round(t_st, 2)}')print(f'Рассчитанный t-критерий дает достигаемый \уровень значимости (p-value) равный {round(p_val, 3)}')

t-критерий составил 2.79

Рассчитанный t-критерий дает достигаемый уровень значимости (p-value) равный 0.007

Вывод для варианта 2

Нулевая гипотезаH₀о том, что средняя плотность скелета грунта не изменилась после транспортировки,отвергаетсяна уровне значимости 0,05 (достигаемый уровень значимостиp-valueдля сгенерированных выборок составил 0.007).

Для наглядности также давайте интервально оценим разность средних по данным выборкам

print("95%% confidence interval: [%.4f, %.4f]"      % DescrStatsW(sample_1 - sample_2).tconfint_mean())

95% confidence interval: [0.0208, 0.1255]

Так как ноль не попадает в рассматриваемый 95% доверительный интервал, мы можем сделать вывод, что средние значения рассматриваемых выборок отличаются.

5. Итог

В данной статье мы рассмотрели возможность применения языка Python при решении практической задачи в инженерной геологии, с попутным исследованием вопроса о необходимом объеме выборки для проверки гипотез.

Карантинные ограничения продолжают сниматься, и компании готовятся к отмене удаленной работы. Однако сотрудники не торопятся возвращаться.

Профессионалы в области технологий не исключение. По данным опроса Hackajob, эксклюзивно поделившейся его результатами с City A. M., 86% из них хотят продолжить работать на дому после пандемии.

Только 14% из 1700 опрошенных технических специалистов хотят вернуться в офис компании на полный рабочий день, в то время как примерно каждый четвертый хотел бы работать удаленно на постоянной основе.

Шестьдесят процентов с удовольствием работают из офиса время от времени и проводят остаток недели, работая дома.

"Гибридная форма работы это новый прорыв для IT специалистов",-подчеркнул Марк Чаффи, соучредитель и генеральный директор Hackajob.

"Хотя работа из дома, возможно, не была самой легкой для людей в прошлом году, специалисты явно высоко ценят возможность не появляться в офисе каждый день. Сотрудники чувствуют себя более комфортно и счастливо, работая из дома и соблюдая баланс между работой и личной жизнью", - сказал Чаффи в интервью City AM.

Он добавил, что технические специалисты не теряют в продуктивности, работая из дома, даже более того, становятся эффективнее благодаря меньшему количеству отвлекающих факторов и отсутствию необходимости в поездках до офиса.

Двигаясь дальше

Также весьма любопытно, что 80% технических специалистов не думают, что продолжат работать в той же компании через два года.

Это говорит о том, что сотрудники стали больше заботиться о гибкости работы и личном удовлетворении, чем о стабильности, льготах и пенсиях, сказал Чаффи.

"Для беби-бумеров и поколения X было вполне нормальным оставаться в одной и той же организации большую часть своей карьеры. Миллениалы и поколение Z, наоборот, предпочитают концепцию частой смены места работы", - продолжил он.

"Тенденция смены работы каждый год или два может рассматриваться работодателями как тревожный звонок, но, на самом деле, это показывает амбициозность, адаптивность и хорошую осведомлённость человека характеристики, которые ищет каждая компания."

Процессы найма

Опрос также подчеркивает необходимость для потенциальных работодателей отказаться от резюме и устранить предвзятость, которая все еще существует в процессе найма.

Более половины (57%) технических специалистов хотят, чтобы их оценивали по их конкретным навыкам, а не по их резюме или по их полу, этнической принадлежности, образованию, сексуальности, инвалидности или социально-экономическому статусу.

Кроме того, многие говорят, что процессы собеседования вслепую имеют решающее значение, если организации стремятся нанимать людей на основе их потенциала и стремиться к созданию разнообразной и инклюзивной культуры на рабочем месте.

В современном мире мы обладаем все большим и большим объемом данных о событиях, происходящих вокруг. Зачастую у нас появляются вопросы, на которые хотелось бы быстро ответить на основе имеющейся информации, для этого как нельзя лучше подходит процесс, связанный с проверкой статистических гипотез. Однако, многие считают, что это занятие подразумевает под собой большое число вычислений и в принципе довольно сложно для понимания. На самом деле, алгоритм проверки гипотез достаточно прост, а для осуществления расчетов с каждым годом появляется все больше и больше готовых инструментальных средств, не требующих от человека глубоких познаний в области. Далее я попытаюсь показать, что мало того, что процесс проверки гипотез может быть полезным, так и осуществляется достаточно быстро и без серьезных усилий.

Статистические гипотезы и области их применения

Статистическая гипотеза - это предположение о каких-либо характеристиках случайной величины. Например: существенно ли изменение числа AI-стартапов в Европе в два разных года и т. д.

Проверка статистических гипотез является важнейшим классом задач математической статистики. С помощью данного инструмента можно подтвердить или отвергнуть предположение о свойствах случайной величины путем применения методов статистического анализа для элементов выборки. Если в предыдущем предложении какие-либо термины являются не совсем понятными, ниже можно найти пояснение на простом языке.

Случайная величина - это величина, которая в зависимости от той или иной ситуации принимает конкретные значения с определенными вероятностями. Примеры: отметка на экзамене; результат игры в кости; количество AI-стартапов по странам Европы. В общем, почти все что угодно!

Генеральная совокупность - совокупность всех объектов для анализа. Например: все AI-стартапы в Европе в 2019-м году.

Выборка - часть данных из генеральной совокупности. Например: официально зарегистрированные AI-стартапы в некоторых странах Европы в 2019-м году.

Статистический анализ - использование различных методов для того, чтобы определить свойства генеральной совокупности по выборке.

Для проверки статистических гипотез зачастую применяются статистические тесты, о которых будет рассказано далее.

Алгоритм проверки статистической гипотезы

В обобщенном виде алгоритм выглядит таким образом:

1. Формулировка основной (H0) и альтернативной (H1) гипотез

2. Выбор уровня значимости

3. Выбор статистического критерия

4. Определения правила принятия решения

5. Итоговое принятие решения на основе исходной выборки данных

Данные шаги являются унифицированными и схему можно использовать почти во всех случаях. Далее подробнее рассмотрим пример работы данного алгоритма на конкретных данных.

Пример проверки статистической гипотезы

Итак, как вы, наверное, догадались по вышеприведенным примерам, будем проверять гипотезу о том, что имеется существенное различие между числом созданных европейских AI-стартапов в 2019-м и 2020-м годах. Пример достаточно простой, чтобы было проще разобраться в ходе работы алгоритма.

Сначала обратим внимание на исходную выборку (рис. 1): датасет представлен для 30-ти Европейских стран, внесены только официально зарегистрированные в стране стартапы. Данные количественные по двум годам. Стоит отметить, что выборки - парные, то есть мы наблюдаем один и тот же показатель для одних и тех же стран с разницей в год.

Сразу стоит отметить, что будут проверены две статистические гипотезы подряд. Для того, чтобы применять критерий для сравнения средних выборок двух лет нужно сначала определить закон распределения данных. Таким образом, шаг 1 - проверка статистической гипотезы о законе распределения данных. Шаг 2 - проверка статистической гипотезы о равенстве между средними.

Проверка гипотезы о законе распределения

Для данных 2019-го года проверим нормальность распределения.

1. H0: случайная величина распределена нормально

   H1: случайная величина не распределена нормально

2. Пусть уровень значимости alpha = 0.05 (как и в 95-ти процентах статистических тестов). Определение уровня значимости достойно отдельного поста, так что не будем заострять на нем внимание.

3. Будет использован критерий Шапиро-Уилка.

4. На этом шаге необходимо разобраться, как работает критерий. В данном случае рассчитывается следующая статистика - функция от нашей выборки:

   , , , ;

   Как видно, формула не слишком простая, плюс существует непростой механизм определения параметра a, поэтому в таких случаях проще пользоваться онлайн-калькуляторами для расчета статистики. Я, например, воспользуюсь хорошим статистическим онлайн-ресурсом - https://www.statskingdom.com/320ShapiroWilk.html.

   Итак, калькулятор показал нам, что p-value =1.20005e-9 , W = 0.435974; Что же делать дальше? Есть два варианта:

   Можно сравнить статистику W с критическим значением Wкрит. Критическое значение чаще всего приведено в готовых таблицах (по строкам/столбцам там отмечен объем выборки и уровень значимости, а на пересечении как раз-таки и лежит Wкрит.). Если W>Wкрит., то не отвергаем H0 и наоборот. Но это не очень удобно, поэтому чаще используется второй способ.

   Можно сравнить p-value с alpha (выбран на 2-ом шаге). Если p-value < alpha, то отвергаем H0. Если нет, то НЕ отвергаем H0. В нашем случае p-value < alpha, следовательно с 95%-ой уверенностью отвергаем H0.

5. H0 отклонена, распределение выборочных данных за 2019-й год не подчинено нормальному закону распределения.

Для данных 2020-го года проверим нормальность распределения. Здесь шаги абсолютно те же самые. Получилось, что p-value = 3.41343e-9. Значение p-value < alpha, следовательно отвергаем H0.

Таким образом, значения в обеих выборках распределены не нормально. Для сравнения средних в двух годах будем использовать критерий Вилкоксона.

Проверка гипотезы о различии в числе AI-стартапов в европейских странах для 2019-го и 2020-го годов

1. H0: отсутствует статистически значимое различие между числом AI-стартапов в Европе в двух годах.

   H1: признается статистическая значимость изменения показателя числа AI-стартапов в Европе между 2019-м и 2020-м годами.

2. Пусть уровень значимости alpha = 0.05.

3. Будет использован критерий Вилкоксона.

4. На этом шаге необходимо разобраться, как работает критерий. Безусловно, для данного критерия также существуют онлайн-калькуляторы, но его достаточно просто посчитать и вручную. Алгоритм очень прост:

   Шаг 1 - Для каждой страны нужно вычислить разность между значениями двух лет.

   Шаг 2 - Далее понять, какие из разностей являются типичными, то есть соответствуют преобладающему по частоте направлению изменения показателя.

   Шаг 3 - Далее в порядке возрастания проранжировать разности пар по их абсолютным значениям. Меньшему абсолютному значению разности приписывается меньший ранг.

   Шаг 4 - Рассчитать сумму рангов, соответствующих нетипичным сдвигам. Это и будет значением T-критерия.

   Пример расчета для двенадцати стран приведен на рисунке ниже (рис. 2). Не пугайтесь, приведенные ранги рассчитаны по всем 30-ти элементам выборки, двенадцать стран приведены лишь для иллюстрации. Проведя такой расчет по всем 30-ти странам и сложив ранги для стран с нетипичными изменениями, получилось, что T = 28.

   Сравним T и Tкрит.=163. T < Tкрит, значит с 95-ой уверенностью изменение числа стартапов статистически значимо.

5. H0 отвергается, различия между числом европейских AI-стартапов в 2019-м и 2020-м годах существенны.

Разнообразие статистических критериев

Как мы увидели на примере, важным шагом в проверке статистической гипотезы является выбор критерия. В примере выше я использовала лишь два статистических критерия, но по факту их гораздо больше, так сказать, на все случаи жизни. Данные критерии важно знать и четко нужно осознавать, когда и какой можно применить. Многие из них направлены на сравнение центров распределений случайных величин, например, сравнение средних, медиан, равенство параметра распределения какому-либо числу и т. д. В основном они делятся на параметрические (знаем закон распределения случайной величины) и непараметрические.

Для вашего удобства внизу (рис. 3) приведена таблица с основными, с моей точки зрения, критериями сравнения центров распределения и их классификацией. Надеюсь, она будет вам полезна, ее можно дополнять и расширять по вашему желанию.