Когнитивные технологии

Алла Тамбовцева

Библиотека scipy для статистики

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

Импортируем библиотеку pandas и модуль stats из библиотеки scipy.

In [5]:
import scipy.stats as st
import pandas as pd

Загрузим базу данных (датафрейм) из csv-файла. Описание данных см.здесь.

In [6]:
df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')
In [7]:
df.head()
Out[7]:
Unnamed: 0 Fertility Agriculture Examination Education Catholic Infant.Mortality
0 Courtelary 80.2 17.0 15 12 9.96 22.2
1 Delemont 83.1 45.1 6 9 84.84 22.2
2 Franches-Mnt 92.5 39.7 5 5 93.40 20.2
3 Moutier 85.8 36.5 12 7 33.77 20.3
4 Neuveville 76.9 43.5 17 15 5.16 20.6

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

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

In [8]:
df[df.Catholic > 50] # для иллюстрации
Out[8]:
Unnamed: 0 Fertility Agriculture Examination Education Catholic Infant.Mortality
1 Delemont 83.1 45.1 6 9 84.84 22.2
2 Franches-Mnt 92.5 39.7 5 5 93.40 20.2
5 Porrentruy 76.1 35.3 9 7 90.57 26.6
6 Broye 83.8 70.2 16 7 92.85 23.6
7 Glane 92.4 67.8 14 8 97.16 24.9
8 Gruyere 82.4 53.3 12 7 97.67 21.0
9 Sarine 82.9 45.2 16 13 91.38 24.4
10 Veveyse 87.1 64.5 14 6 98.61 24.5
30 Conthey 75.5 85.9 3 2 99.71 15.1
31 Entremont 69.3 84.9 7 6 99.68 19.8
32 Herens 77.3 89.7 5 2 100.00 18.3
33 Martigwy 70.5 78.2 12 6 98.96 19.4
34 Monthey 79.4 64.9 7 3 98.22 20.2
35 St Maurice 65.0 75.9 9 9 99.06 17.8
36 Sierre 92.2 84.6 3 3 99.46 16.3
37 Sion 79.3 63.1 13 13 96.83 18.1
45 Rive Droite 44.7 46.6 16 29 50.43 18.2
46 Rive Gauche 42.8 27.7 22 29 58.33 19.3
In [9]:
sample1 = df[df.Catholic > 50]["Infant.Mortality"] # выборка 1
sample1
Out[9]:
1     22.2
2     20.2
5     26.6
6     23.6
7     24.9
8     21.0
9     24.4
10    24.5
30    15.1
31    19.8
32    18.3
33    19.4
34    20.2
35    17.8
36    16.3
37    18.1
45    18.2
46    19.3
Name: Infant.Mortality, dtype: float64
In [10]:
sample2 = df[df.Catholic <= 50]["Infant.Mortality"] # выборка 2
sample2
Out[10]:
0     22.2
3     20.3
4     20.6
11    16.5
12    19.1
13    22.7
14    18.7
15    21.2
16    20.0
17    20.2
18    10.8
19    20.0
20    18.0
21    22.4
22    16.7
23    15.3
24    21.0
25    23.8
26    18.0
27    16.3
28    20.9
29    22.5
38    20.3
39    20.5
40    18.9
41    23.0
42    20.0
43    19.5
44    18.0
Name: Infant.Mortality, dtype: float64

Теперь приступим к формальной проверке гипотез.

T-test (критерий Стьюдента для двух выборок)

T-test используется для сравнения средних значений двух генеральных совокупностей в предположении, что обе выборки взяты из нормального распределения. Существует две разновидности t-теста: t-тест для независимых выборок и t-тест для связных (парных) выборок. В связных выборках объекты связаны друг с другом. Пример связных выборок: значения уровня смертности в одних и тех же кантонах до и после какой-нибудь реформы)

Нулевая гипотеза: средние значения двух генеральных совокупностей, откуда взяты выборки, равны, то есть $H_0: a_1=a_2$.

Альтернативная гипотеза: средние значения двух генеральных совокупностей не равны, то есть: $H_1: a_1 \ne a_2$

В stats в t-тесте в качестве альтернативной гипотезы используется двусторонняя альтернатива (средние не равны) и всегда выводится соответствующее p-value (two-tailed). То же будет характерно для всех последующих тестов. Так как наши выборки независимы, нам нужна функция ttest_ind(), от independent.

In [11]:
st.ttest_ind(sample1, sample2)
Out[11]:
Ttest_indResult(statistic=1.1297965130690064, pvalue=0.26454837328688746)

Что возвращает эта функция? Наблюдаемое значение t-статистики и p-value. Результат понятный, только более лакончиный по сравнению с тем, что выводит R и другие статистические пакеты.

Выводы: так как p-value больше любого конвенционального уровня значимости (1%, 5%, 10%), на имеющихся данных на любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу. Средний уровень детской смертности в католических и протестантских районах можно считать одинаковым.

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

In [12]:
st.ttest_ind(sample1, sample2, equal_var = False)
Out[12]:
Ttest_indResult(statistic=1.0863471703503398, pvalue=0.28551301767919196)

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

Нулевая гипотеза: дисперсии двух генеральных совокупностей равны, то есть $H_0: \sigma_1^2 = \sigma_2^2$

Альтернативная гипотеза: дисперсии двух генеральных совокупностей не равны, то есть $H_1: \sigma_1^2 \ne \sigma_2^2$.

Реализовать F-тест, который нам нужен именно для этого случая, в Python сразу не получится: встроенная функция f_oneway() используется для однофакторного дисперсионного анализа (ANOVA), речь о котором пойдёт далее. Можно попробовать реализовать этот тест «вручную», рассчитав частное выборочных дисперсий и поработав с F-распределением, но давайте пойдём другим путём и воспользуемся другими критериями для сравнения дисперсий, которые явно встроены в Python.

Например, тест Левена:

In [14]:
st.levene(sample1, sample2)
Out[14]:
LeveneResult(statistic=1.0380526498527436, pvalue=0.3137214012581385)

На любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу о равенстве дисперсий.

Wilcoxon test & Mann-Whitney test (две выборки)

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

In [15]:
st.wilcoxon(sample1, sample2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-6bfff17a9e4f> in <module>()
----> 1 st.wilcoxon(sample1, sample2)

/usr/local/lib/python3.5/dist-packages/scipy/stats/morestats.py in wilcoxon(x, y, zero_method, correction)
   2374         x, y = map(asarray, (x, y))
   2375         if len(x) != len(y):
-> 2376             raise ValueError('Unequal N in wilcoxon.  Aborting.')
   2377         d = x - y
   2378 

ValueError: Unequal N in wilcoxon.  Aborting.

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

Нулевая гипотеза: выборки взяты из одного и того же распределения, то есть $H_0: F(x) = G(x)$

Альтернативная гипотеза: выборки взяты из разных распределений, то есть $H_1: F(x) \ne G(x)$.

In [16]:
st.mannwhitneyu(sample1, sample2)
Out[16]:
MannwhitneyuResult(statistic=235.5, pvalue=0.2920645577220585)

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

ANOVA

Если выборок больше, чем две, то использовать указанные выше критерии нельзя. В предположении, что все выборки взяты из нормального распределения, для сравнения средних значений более чем в двух группах используется однофакторный дисперсионный анализ (ANOVA, analysis of variance).

Нулевая гипотеза: средние значения по всем $k$ группам (во всех генеральных совокупностях) равны, то есть $H_0: a_1 = a_2 = \dots = a_k$

Альтернативная гипотеза: средние значения по всем группам (во всех генеральных совокупностях) не равны.

Чтобы не создавать искусственные группы на основе данных в swiss, загрузим таблицу с весами цыплят, которых кормили разным кормом :) Описание см. здесь.

In [7]:
dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')
In [8]:
dat.head()
Out[8]:
Unnamed: 0 weight feed
0 1 179 horsebean
1 2 160 horsebean
2 3 136 horsebean
3 4 227 horsebean
4 5 217 horsebean

Задание: разбить датафрейм на группы с помощью groupby по переменной feed и сохранить значения weight в словарь.

Решение:

In [9]:
wgt = {}
for name, d in dat.groupby('feed'):
    wgt[name] = d.weight
In [10]:
wgt
Out[10]:
{'casein': 59    368
 60    390
 61    379
 62    260
 63    404
 64    318
 65    352
 66    359
 67    216
 68    222
 69    283
 70    332
 Name: weight, dtype: int64, 'horsebean': 0    179
 1    160
 2    136
 3    227
 4    217
 5    168
 6    108
 7    124
 8    143
 9    140
 Name: weight, dtype: int64, 'linseed': 10    309
 11    229
 12    181
 13    141
 14    260
 15    203
 16    148
 17    169
 18    213
 19    257
 20    244
 21    271
 Name: weight, dtype: int64, 'meatmeal': 48    325
 49    257
 50    303
 51    315
 52    380
 53    153
 54    263
 55    242
 56    206
 57    344
 58    258
 Name: weight, dtype: int64, 'soybean': 22    243
 23    230
 24    248
 25    327
 26    329
 27    250
 28    193
 29    271
 30    316
 31    267
 32    199
 33    171
 34    158
 35    248
 Name: weight, dtype: int64, 'sunflower': 36    423
 37    340
 38    392
 39    339
 40    341
 41    226
 42    320
 43    295
 44    334
 45    322
 46    297
 47    318
 Name: weight, dtype: int64}

Теперь ANOVA (f_oneway от One-Way ANOVA):

In [11]:
st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
Out[11]:
F_onewayResult(statistic=15.364799774712534, pvalue=5.936419853471331e-10)

Функция возвращает наблюдаемое значение F-статистики и p-value. В данном случае p-value близко к 0, поэтому гипотезу о равенстве средних генеральных совокупностей по группам можно отвергнуть на 1% уровне значимости. Средний вес цыплят, которых кормили разным кормом, отличается (ещё бы, horsebean или sunflower!).

Kruskal-Wallis test

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

Нулевая гипотеза: выборки взяты из одного и того же распределения, то есть: $H_0: F(x) = G(x) = \dots = H(x)$

Альтернативная гипотеза: выборки взяты из разных распределений.

In [12]:
st.kruskal(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
Out[12]:
KruskalResult(statistic=37.34271769425624, pvalue=5.112829511937094e-07)