Статистика

Серия из 5 постов для начинающих представляет собой ремикс первой главы книги 2015 года под названием Clojure для науки о данных (Clojure for Data Science). Автор книги, Генри Гарнер, любезно дал согласие на использование материалов книги для данного ремикса с использованием языка Python.

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

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

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

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

Пост 1 посвящен подготовке среды и данных.

Статистика

Важно не кто голосует, а кто подсчитывает голоса

Иосиф Сталин

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

При изучении распределений чрезвычайную важность играет наглядная и удобная визуализация данных, и для этого мы воспользуемся Python-овской библиотекой pandas. Мы покажем, как пользоваться ею для загрузки, преобразования и разведывательного анализа реальных данных, а также начнем работать с фундаментальной библиотекой numpy для научных вычислений. Мы проведем сопоставительный анализ результатов двух общенациональных выборов всеобщих выборов в Великобритании 2010 г. и российских выборов депутатов Государственной Думы Федерального Собрания РФ шестого созыва 2011 г. и увидим, каким образом даже элементарный анализ может предъявить подтверждающие данные о потенциальных фальсификациях.

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

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

Кроме того, мы будем пользоваться встроенными в Python модулями. Так, например, модуль random позволяет генерировать случайные числа и извлекать выборки, и модуль collections содержит дополнительные структуры данных, из которых мы воспользуемся специальным словарем Counter.

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

• Если данные представлены текстовым файлом с разделением полей данных запятыми (.csv) или символами табуляции (.tsv), то мы будем использовать функцию чтения данных read_csv

• Если данные представлены файлом Excel (например, файл .xls или .xlsx), то мы воспользуемся функцией чтения данных read_excel

• Для любого другого источника данных (внешняя база данных, веб-сайт, буфер обмена данными, JSON-файлы, HTML-файлы и т. д.) предусмотрен ряд других функций

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

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

pd.read_excel('data/ch01/UK2010.xls')

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

def load_uk(): '''Загрузить данные по Великобритании''' return pd.read_excel('data/ch01/UK2010.xls') 

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

Первая строка электронной таблицы UK2010.xls содержит имена столбцов. Функция библиотеки pandas read_excel резервирует их в качестве имен столбцов возвращаемого кадра данных. Начнем обследование данных с их проверки атрибут кадра данных columns возвращает имена столбцов в виде списка, при этом адресация атрибутов осуществляется при помощи оператора точки (.):

def ex_1_1(): '''Получить имена полей кадра данных''' return load_uk().columns

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

Index(['Press Association Reference', 'Constituency Name', 'Region', 'Election Year', 'Electorate', 'Votes', 'AC', 'AD', 'AGS', 'APNI', ... 'UKIP', 'UPS', 'UV', 'VCCA', 'Vote', 'Wessex Reg', 'WRP', 'You', 'Youth', 'YRDPL'], dtype='object', length=144)

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

• Информация для Ассоциации прессы: число, идентифицирующее избирательный округ (представленный одним депутатом)

• Название избирательного округа: стандартное название, данное избирательному округу

• Регион: географический район Великобритании, где округ расположен

• Год выборов: год, в котором выборы состоялись

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

• Голосование: общее число проголосовавших

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

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

def ex_1_2(): '''Получить значения поля "Год выборов"''' return load_uk()['Election Year']

В результате будет выведен следующий список:

0 2010.01 2010.02 2010.0...646 2010.0647 2010.0648 2010.0649 2010.0650 NaNName: Election Year, dtype: float64

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

def ex_1_3(): '''Получить значения в поле "Год выборов" без дубликатов''' return load_uk()['Election Year'].unique()

[ 2010. nan]

Значение 2010 еще больше подкрепляет наши ожидания в отношении того, что эти данные относятся к 2010 году. Впрочем, наличие специального значения nan, от англ. not a number, т.е. не число, которое сигнализирует о пропущенных данных, является неожиданным и может свидетельствовать о проблеме с данными.

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

def ex_1_4(): '''Рассчитать частоты в поле "Год выборов"  (количества появлений разных значений)''' return Counter( load_uk()['Election Year'] )

Counter({nan: 1, 2010.0: 650})

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

Исправление данных

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

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

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

def ex_1_5(): '''Вернуть отфильтрованную по полю "Год выборов"  запись в кадре данных (в виде словаря)''' df = load_uk() return df[ df['Election Year'].isnull() ]

Выражение dt['Election Year'].isnull() вернет булеву последовательность, в которой все элементы, кроме последнего, равны False, в результате чего будет возвращена последняя запись кадра данных. Если Вы знаете язык запросов SQL, то отметите, что этот метод очень похож на условный оператор WHERE.

Присмотревшись к результатам примера ex_1_5, можно заметить, что в полученной записи все поля (кроме одного) имеют значение NaN. Дальнейший анализ данных подтверждает, что строка с непустым полем на самом деле является строкой итоговой суммы в листе файла Excel. Эту строку следует из набора данных удалить. Мы можем удалять проблемные строки путем обновления коллекции предикативной функцией notnull(), которая в данном случае вернет только те строки, в которых год выборов не равен NaN:

 df = load_uk() return df[ df[ 'Election Year' ].notnull() ]

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

def load_uk_scrubbed(): '''Загрузить и отфильтровать данные по Великобритании''' df = load_uk() return df[ df[ 'Election Year' ].notnull() ]

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

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

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

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

Булочник и Пуанкаре

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

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

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

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

Генерирование распределений

В целях развития нашего интуитивного понимания относительно нормального распределения и дисперсии, давайте смоделируем честного и нечестного булочников, и для этого воспользуемся функцией генерирования нормально распределенных случайных величин stats.norm.rvs. (rvs от англ. normal variates, т.е. случайные величины). Честного булочника можно смоделировать в виде нормального распределения со средним значением 1000, что соответствует справедливой буханке хлеба весом 1 кг. При этом мы допустим наличие дисперсии в процессе выпекания, которая приводит к стандартному отклонению в 30г.

def honest_baker(mu, sigma): '''Модель честного булочника''' return pd.Series( stats.norm.rvs(loc, scale, size=10000) )def ex_1_18(): '''Смоделировать честного булочника на гистограмме''' honest_baker(1000, 30).hist(bins=25) plt.xlabel('Честный булочник')  plt.ylabel('Частота') plt.show()

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

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

def dishonest_baker(mu, sigma): '''Модель нечестного булочника''' xs = stats.norm.rvs(loc, scale, size=10000)  return pd.Series( map(max, bootstrap(xs, 13)) ) def ex_1_19(): '''Смоделировать нечестного булочника на гистограмме''' dishonest_baker(950, 30).hist(bins=25) plt.xlabel('Нечестный булочник')  plt.ylabel('Частота') plt.show()

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

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

Асимметрия

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

Положительная и отрицательная асимметрии

Библиотека pandas располагает функцией skew для измерения асимметрии:

def ex_1_20(): '''Получить коэффициент асимметрии нормального распределения''' s = dishonest_baker(950, 30) return { 'среднее' : s.mean(),  'медиана' : s.median(),  'асимметрия': s.skew() }

{'асимметрия': 0.4202176889083849, 'медиана': 998.7670301469957, 'среднее': 1000.059263920949}

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

Графики нормального распределения

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

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

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

def qqplot( xs ): '''Квантильный график (график квантиль-квантиль, Q-Q plot)''' d = {0:sorted(stats.norm.rvs(loc=0, scale=1, size=len(xs))), 1:sorted(xs)} pd.DataFrame(d).plot.scatter(0, 1, s=5, grid=True) df.plot.scatter(0, 1, s=5, grid=True) plt.xlabel('Квантили теоретического нормального распределения') plt.ylabel('Квантили данных') plt.title ('Квантильный график', fontweight='semibold')def ex_1_21(): '''Показать квантильные графики  для честного и нечестного булочников''' qqplot( honest_baker(1000, 30) ) plt.show() qqplot( dishonest_baker(950, 30) ) plt.show()

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

Выше показан квантильный график для честного булочника. Далее идет квантильный график для нечестного булочника:

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

^{Надписи: нормально распределенные, тяжелые хвосты, легкие хвосты, скошенность влево, скошенность вправо, раздельные кластеры}

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

Технические приемы сопоставительной визуализации

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

Коробчатые диаграммы

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

def ex_1_22(): '''Показать коробчатую диаграмму с данными честного и нечестного булочников''' d = {'Честный булочник' :honest_baker(1000, 30), 'Нечестный булочник':dishonest_baker(950, 30)}  pd.DataFrame(d).boxplot(sym='o', whis=1.95, showmeans=True) plt.ylabel('Вес буханки (гр.)') plt.show()

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

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

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

Интегральные функции распределения

Интегральные функции распределения (ИФР), также именуемые кумулятивными функциями распределения, от англ. Cumulative Distribution Function (CDF), описывают вероятность, что значение, взятое из распределения, будет меньше x. Как и все распределения вероятностей, их значения лежат в диапазоне между 0 и 1, где 0 это невозможность, а 1 полная определенность. Например, представьте, что я собираюсь бросить шестигранный кубик. Какова вероятность, что выпадет значение меньше 6?

Для уравновешенного кубика вероятность выпадения пятерки или меньшего значения равна 5/6. И наоборот, вероятность, что выпадет единица, равна всего1/6. Тройка или меньше соответствуют равным шансам то есть вероятности 50%.

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

^{ИФР и квантили тесно друг с другом связаны ИФР является инверсией квантильной функции. Если 0.5-квантиль соответствует значению 1000, тогда ИФР для 1000 составляет 0.5.}

Подобно тому, как функция pandas quantile позволяет нам отбирать значения из распределения в конкретных точках, эмпирическая ИФР empirical_cdf позволяет нам внести значение из последовательности и вернуть значение в диапазоне между 0 и 1. Это функция более высокого порядка, т.е. она принимает значение (в данном случае последовательность значений) и возвращает функцию, которую потом можно вызывать, сколько угодно, с различными значениями на входе, и возвращая ИФР для каждого из них.

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

Построим график ИФР одновременно для честного и нечестного булочников. Для этих целей можно воспользоваться функцией библиотеки pandas построения двумерного графика plot для визуализации ИФР, изобразив на графике исходные данные то есть выборки из распределений честного и нечестного булочников в сопоставлении с вероятностями, вычисленными относительно эмпирической ИФР. Функция plot ожидает, что значения x и значения y будут переданы в виде двух раздельных последовательностей значений. Для этих целей мы воспользуемся конструктором кадра данных pandas DataFrame.

Чтобы изобразить оба распределения на одном графике, мы должны передать функции plot несколько серий. Для многих своих графиков pandas предоставляет функции, которые позволяют добавлять дополнительные серии. В случае с функцией plot мы можем присвоить указатель на создаваемый график, присвоив временной переменной (ax) результат первого вызова функции plot, и затем при повторных вызовах указывать эту переменную в именованном аргументе функции (ax=ax). Можно также передать необязательную метку серии. Мы выполним это в следующем ниже примере, чтобы на готовом графике отличить две серии друг от друга. Сначала определим универсальную функцию построения эмпирической ИФР против теоретической, которая получает на вход кортеж из двух серий (tp[1] и tp[3]) и их названий и метки осей, и затем вызовем ее:

def empirical_cdf(x): """Вернуть эмпирическую ИФР для x""" sx = sorted(x) return pd.DataFrame( {0: sx, 1:sp.arange(len(sx))/len(sx)} )def ex_1_23(): '''Показать графики эмпирической ИФР честного булочника в сопоставлении с нечестным''' df = empirical_cdf(honest_baker(1000, 30)) df2 = empirical_cdf(dishonest_baker(950, 30)) 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()

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

Несмотря на то, что этот график выглядит совсем по-другому, он в сущности показывает ту же самую информацию, что и коробчатая диаграмма. Мы видим, что две линии пересекаются примерно в медиане 0.5, соответствующей 1000 гр. Линия нечестного булочника обрезается в нижнем хвосте и удлиняется на верхнем хвосте, что соответствует асимметричному распределению.

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

Пост 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, т.е. значимый результат. Разные варианты веб-сайта действительно имеют разные дисперсии, которые нельзя просто объяснить одной лишь случайной ошибкой в выборке.

F-распределение со степенями свободы 19 и 980

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

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 = 20 и степенью свободы df = 5

Приведенный выше рисунок показывает, как 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.

t-распределение, степень свободы = 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-го рода. Правильнее было бы допустить возможность того, что обновленный веб-сайт может в действительности оказаться хуже существующего. Этот подход распределяет нашу одинаково по обоим хвостам распределения и обеспечивает значимый результат, не искаженный под воздействием априорного допущения об улучшении работы обновленного веб-сайта.

Надписи: t-распределение, степень свободы = 298

В действительности в модуле 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, будут матричные операции, нормальное уравнение и коллинеарность.

---

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

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

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

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

В этой статье мы на примерах рассмотрим 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%, что делает сравнение еще более простым. Для сравнения, это пример, который мы видели выше:

Сравнение распределений в CDP. [Рисунок автора]

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"

- Смотреть онлайн-встречу "День открытых дверей"

Вячеслав Ермолин 2 мая 2021 года
Орбитальные запуски 2021 года. Январь-апрель.

Текущая статистика орбитальных запусков 2021 года. На 30 апреля 2021 года..

Апрель США, Китай, Россия и Европа. 11 запусков. Без аварий. 7 человек на орбиту. Выведено 176 спутников и кораблей.

Орбитальные запуски:
Месяц прошел под знаком человек на орбите
два запуска пилотируемых кораблей для смены экипажа МКС
Союз МС-18 и Dragon Crew-2. Семь космонавтов и астронавтов. Запущен на орбиту базовый модуль китайской национальной станции. Первый европейский запуск в этом году. Легкая РН Vega с миссией Rideshare.
Два запуска Starlink от SpaceX и запуск OneWeb.
8 запусков из 11 сделаны в последнюю неделю апреля.

2021-027: 7 апреля. США. Falcon 9 Block 5 | Starlink 23.
Starlink-23 запуск 60 спутников связи. Starlink v1.0. Для системы низкоорбитального коммерческого интернета Starlink от SpaceX.

2021-028: 9 апреля. Китай. CZ-4B | Shiyan-6-03
Запуск Shiyan-6-03, на солнечно-синхронную орбиту (SSO). Возможно технологический спутник.

2021-029: 9 апреля. Россия. Союз 2.1а|Союз МС-18
Пилотируемый запуск экипажа МКС-65 на Международную Космическую Станцию (МКС) на корабле Союз. Сверхбыстрая, 3-х часовая схема полета. Работа на орбите 191 суток. Два члена этой команды (русский и американец) останутся на МКС до года.

2021-030: 22 апреля. США. Falcon 9 | Crew-2.
Пилотируемый запуск корабля Crew Dragon SpaceX. Четыре астронавта МКС-65 для работы на Международной космической станции (МКС). Второй штатный полет к МКС. Корабль Crew Dragon C206-2. Время работы экипажа на орбите планируется до шести месяцев. Возвращение на Землю и приводнение в Атлантическом океане.

2021-031: 25 апреля. Россия. Soyuz 2.1b/Fregat-M | OneWeb 6.
Запуск 36 спутников связи Oneweb Launch#6. Спутники для будущей группировки низкоорбитального интернета Oneweb. Запуск с космодрома Восточный.

2021-032: 26 апреля. США. Delta IV-H | NROL-82.
Спутник Национального управления военно-космической разведки США. Предположительно спутник оптической разведки серии KH-11/Kennen/Crystal.

2021-033: 27 апреля. Китай. CZ-6 | Девять коммерческих спутников.
Запуск 9 коммерческих спутников в первой китайской миссии Rideshare от компании Great Wall Company. Три основных спутника ДЗЗ весом до 100 кг каждый и шесть дополнительных. Выведены на орбиту ССО высотой 500 км.

2021-034: 29 апреля. Европа. Vega VV18 | Rideshare.
Пуск ракеты-носителя Vega. Миссия Rideshare с основной нагрузкой в виде французского спутника-шпиона и пяти дополнительных нагрузок CubSat из трех стран.

2021-035: 29 апреля. Китай. CZ-5B | Tianhe.
Пуск ракеты-носителя Чанчжэн-5B с основном модулем будущей китайской орбитальной станции Тяньгун (Небесный дворец).

2021-036: 29 апреля. США. Falcon 9 Block 5 | Starlink 24.
Starlink-24 запуск 60 спутников связи. Starlink v1.0. Для системы низкоорбитального коммерческого интернета Starlink от SpaceX.

2021-037: 30 апреля. Китай. Long March 4C | Yaogan 34.
Яогань-34. Пуск успешный. Предположительно новый спутник оптической разведки.

Планы орбитальных запусков на 2021 год и мой прогноз

Максимальный прогноз (неофициальный) на 2021 годЛегенда к статистике

Зарабатывать продажей лекарств, которые заведомо не работают, не только аморально, но и не особо легко. Люди всё-таки обычно не хотят покупать препараты, неэффективность которых была доказана. А вот если вы сумели выдавить заветное 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%) технических специалистов хотят, чтобы их оценивали по их конкретным навыкам, а не по их резюме или по их полу, этнической принадлежности, образованию, сексуальности, инвалидности или социально-экономическому статусу.

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

Перевод подготовлен в рамках курса "Machine Learning. Professional".

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

---

Ошибки наблюдения и различия в подгруппах вызывают статистические парадоксы

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

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

1. Парадокс Берксона

Первым ярчайшим примером является обратная корреляция между степенью тяжести заболевания COVID-19 и курением сигарет (см., например, обзор Европейской комиссии Wenzel 2020). Курение сигарет широко известный фактор риска респираторных заболеваний, так как же объяснить это противоречие?

Работа Griffith 2020, недавно опубликованная в Nature, предполагает, что это может быть случай ошибки коллайдера (Collider Bias), также называемой парадоксом Берксона. Чтобы понять этот парадокс, давайте рассмотрим следующую графическую модель, в которую мы включили третью случайную переменную: госпитализация.

Парадокс Берксона: госпитализация это переменная-коллайдер для курения сигарет, и для тяжести течения COVID-19. (Изображение автора)

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

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

Парадокс Берксона: если мы добавим условие в соответствии с коллайдером госпитализация, мы увидим обратную связь между курением и COVID-19! (Изображение автора)

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

Но на правом рисунке где мы рассматриваем только пациентов больниц мы видим противоположную тенденцию! Чтобы понять это, обратите внимание на следующие моменты.

1.Тяжелая форма COVID-19 увеличивает шансы на госпитализацию. То есть, если степень тяжести заболевания выше 1, то требуется госпитализация.

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

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

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

2. Скрытые (латентные) переменные

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

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

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

Парадокс скрытой переменной: степень тяжести пожара это скрытая переменная для n задействованных пожарных и для n пострадавших. (Изображение автора)

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

Давайте рассмотрим следующий пример с набором данных. На левом рисунке у нас отражены общие данные по всем видам пожаров, а на правом рисунке мы рассматриваем только сведения, соответствующие трем фиксированным степеням тяжести пожара (т.е. мы обусловливаем наши данные наблюдений скрытой переменной).

Скрытые переменные: если мы примем за условие скрытую переменную степень тяжести пожара, мы увидим обратную корреляцию между количеством задействованных пожарных и количеством раненых! (Изображение автора)

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

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

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

3. Парадокс Симпсона

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

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

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

• С другой стороны, общий процент приема среди абитуриентов женского пола был ниже, чем у абитуриентов мужского пола.

Чтобы понять, как как такое может быть, давайте рассмотрим следующий набор данных с двумя факультетами: Факультет A и Факультет B.

• Из 100 абитуриентов мужского пола: 80 подали заявки на Факультет A, из которых 68 были приняты (85%), а 20 подали заявки на Факультет В, из которых приняты были 12 человек (60%).

• Из 100 абитуриентов женского пола: 30 подали заявки на Факультет А, из которых 28 были приняты (93%), в то время как 70 подали заявки на Факультет B, из которых были приняты 46 (66%).

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

Парадокс выражается следующими неравенствами.

Парадокс Симпсона: неравенство, лежащее в основе очевидного противоречия. (Изображение автора)

Теперь мы можем понять происхождение наших, казалось бы, противоречивых наблюдений. Дело в том, что существует ощутимый классовый гендерный дисбаланс среди абитуриентов на каждом из двух факультетов (Факультет A: 8030, Факультет B: 2070). Действительно, большинство студентов женского пола подали заявку на более конкурентный Факультет B (который имеет низкие показатели приема), в то время как большинство студентов мужского пола подали документы на менее конкурентный Факультет А (который имеет более высокие показатели приема). Это обусловливает противоречивые данные, которые мы получили.

Заключение

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

---

Узнать подробнее о курсе "Machine Learning. Professional"

Участвовать в онлайн-интенсиве Деплой ML модели: от грязного кода в ноутбуке к рабочему сервису

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

Что изменилось за эти годы? Стало ли лучше или хуже? Работают ли порталы открытых данных? Публикуются ли данные?

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

• за 2020 год на федеральном портале открытых данных (data.gov.ru) было опубликовано 223 набора данных, за 5 месяцев 2021 года - только 2 набора данных

• всего с 2020 года объём этих 225 наборов данных - 405 мегабайт из которых более 390 мегабайт - это данные Минкультуры России и ФНС России (и то есть подозрение что цифры завышены потому что в реестре наборов данных есть дублирующиеся записи. Скорее всего реально данных значительно меньше)

• лишь 9 178 наборов данных из 24 002 опубликованы федеральными органами власти, остальные региональными и муниципальными

• 10 ФОИВов не опубликовали ни одного нового набора данных с 2013 года (за 8 лет)

• 20 ФОИВов не опубликовали ни одного нового набора данных с 2015 года (за 6 лет)

• 42 ФОИВа не опубликовали ни одного нового набора данных с 2017 года (за 4 лет)

• 68 ФОИВов не опубликовали ни одного нового набора данных с 2019 года (за 2 года)

• иначе говоря в 2020 и 2021 года лишь 6 ФОИВов разместили хотя бы один новый набор данных на портале открытых данных

• некоторые ФОИВы, при этом, кое что опубликовали на своих сайтах, но куда меньше чем раньше и чем могли бы

• общий объём опубликованных данных на портале data.gov.ru оценить сложно, сайт не даёт статистики, API сайта очень куцое, требуется очень много запросов сделать чтобы подсчитать хоть самые приблизительные цифры, но они будут невелики.

• параллельно этому на сайтах и FTP серверах органов власти опубликовано открытых данных, оценочно, на 20 терабайт в форме архивов. Количественно - это сотни наборов данных, качественно - это данные большого объёма.

• безусловные лидеры по масштабам раскрытия данных - Минкультуры, ФНС России, Федеральное казначейство, Минфин России. Даже при том что тенденции там не только к раскрытию, текущие объёмы доступных данных очень велики.

• источники наиболее крупных наборов данных:

  • Сайт ФНС России (nalog.ru)

  • Портал открытых данных Минкультуры России (opendata.mkrf.ru)

  • Портал госзакупок (zakupki.gov.ru)

  • Единый портал бюджетной системы (budget.gov.ru)

  • ФИАС (fias.nalog.ru)

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

• в ряде субъектов федерации закрывают порталы открытых данных. Например, его закрыли в Московской области

• в других субъектах федерации перестали публиковать новые наборы данных, это касается и портала открытых данных Москвы и десятков других субъектов федерации (чуть ли не всех)

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

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