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

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

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

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

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

import scipy.stats as st
import pandas as pd

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

df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')

df.head()

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

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

df[df.Catholic > 50] # для иллюстрации

sample1 = df[df.Catholic > 50]["Infant.Mortality"] # выборка 1
sample1

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

sample2 = df[df.Catholic <= 50]["Infant.Mortality"] # выборка 2
sample2

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.

st.ttest_ind(sample1, sample2)

Ttest_indResult(statistic=1.1297965130690064, pvalue=0.26454837328688746)

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

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

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

st.ttest_ind(sample1, sample2, equal_var = False)

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.

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

st.levene(sample1, sample2)

LeveneResult(statistic=1.0380526498527436, pvalue=0.3137214012581385)

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

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

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

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)$.

st.mannwhitneyu(sample1, sample2)

MannwhitneyuResult(statistic=235.5, pvalue=0.2920645577220585)

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

ANOVA¶

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

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

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

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

dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')

dat.head()

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

Решение:

wgt = {}
for name, d in dat.groupby('feed'):
    wgt[name] = d.weight

wgt

{'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):

st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])

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)$

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

st.kruskal(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])

KruskalResult(statistic=37.34271769425624, pvalue=5.112829511937094e-07)

	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

	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

	Unnamed: 0	weight	feed
0	1	179	horsebean
1	2	160	horsebean
2	3	136	horsebean
3	4	227	horsebean
4	5	217	horsebean