In [65]:
from __future__ import division
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.proportion import proportion_confint
from statsmodels.stats.weightstats import CompareMeans, DescrStatsW, ztest
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
В одном из выпусков программы "Разрушители легенд" проверялось, действительно ли заразительна зевота. В эксперименте участвовало 50 испытуемых, проходивших собеседование на программу. Каждый из них разговаривал с рекрутером; в конце 34 из 50 бесед рекрутер зевал. Затем испытуемых просили подождать решения рекрутера в соседней пустой комнате.
Во время ожидания 10 из 34 испытуемых экспериментальной группы и 4 из 16 испытуемых контрольной начали зевать. Таким образом, разница в доле зевающих людей в этих двух группах составила примерно 4.4%. Ведущие заключили, что миф о заразительности зевоты подтверждён.
Можно ли утверждать, что доли зевающих в контрольной и экспериментальной группах отличаются статистически значимо? Посчитайте достигаемый уровень значимости при альтернативе заразительности зевоты, округлите до четырёх знаков после десятичной точки.
In [2]:
data_exp = np.array( [1 if i < 10 else 0 for i in range(34)] )
data_ctrl = np.array( [1 if i < 4 else 0 for i in range(16)] )
In [10]:
print('Mean experimental value: %.4f' % data_exp.mean())
print('Mean control value: %.4f' % data_ctrl.mean())
In [12]:
conf_interval_banner_exp = proportion_confint(np.sum(data_exp), len(data_exp), method = 'wilson')
conf_interval_banner_ctrl = proportion_confint(np.sum(data_ctrl), len(data_ctrl), method = 'wilson')
In [13]:
print('95%% confidence interval for exp group: [%f, %f]' % conf_interval_banner_exp)
print('95%% confidence interval for a ctrl group: [%f, %f]' % conf_interval_banner_ctrl)
In [18]:
def proportions_diff_confint_ind(sample1, sample2, alpha = 0.05):
z = stats.norm.ppf(1 - alpha / 2.)
p1 = float(sum(sample1)) / len(sample1)
p2 = float(sum(sample2)) / len(sample2)
left_boundary = (p1 - p2) - z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
right_boundary = (p1 - p2) + z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
return (left_boundary, right_boundary)
In [19]:
def proportions_diff_z_stat_ind(sample1, sample2):
n1 = len(sample1)
n2 = len(sample2)
p1 = float(sum(sample1)) / n1
p2 = float(sum(sample2)) / n2
P = float(p1*n1 + p2*n2) / (n1 + n2)
return (p1 - p2) / np.sqrt(P * (1 - P) * (1. / n1 + 1. / n2))
In [27]:
def proportions_diff_z_test(z_stat, alternative='two-sided'):
if alternative not in ('two-sided', 'less', 'greater'):
raise ValueError("alternative not recognized\n"
"should be 'two-sided', 'less' or 'greater'")
if alternative == 'two-sided':
return 2 * (1 - stats.norm.cdf(np.abs(z_stat)))
if alternative == 'less':
return stats.norm.cdf(z_stat)
if alternative == 'greater':
return 1 - stats.norm.cdf(z_stat)
In [28]:
print('95%% confidence interval for a difference: [%.4f, %.4f]' % proportions_diff_confint_ind(data_exp, data_ctrl))
In [29]:
print('p-value: %.4f' % proportions_diff_z_test(proportions_diff_z_stat_ind(data_exp, data_ctrl), 'greater'))
Имеются данные измерений двухсот швейцарских тысячефранковых банкнот, бывших в обращении в первой половине XX века. Сто из банкнот были настоящими, и сто — поддельными. На рисунке ниже показаны измеренные признаки:
Загрузите данные:
banknotes.txt
Отделите 50 случайных наблюдений в тестовую выборку с помощью функции sklearn.cross_validation.train_test_split (зафиксируйте random state = 1). На оставшихся 150 настройте два классификатора поддельности банкнот:
логистическая регрессия по признакам X1,X2,X3; логистическая регрессия по признакам X4,X5,X6. Каждым из классификаторов сделайте предсказания меток классов на тестовой выборке. Одинаковы ли доли ошибочных предсказаний двух классификаторов? Проверьте гипотезу, вычислите достигаемый уровень значимости. Введите номер первой значащей цифры (например, если вы получили 5.5×10−8, нужно ввести 8).
In [32]:
banknotes_data = pd.read_table('banknotes.txt')
banknotes_data.describe()
banknotes_data.info()
banknotes_data.head()
Out[32]:
Out[32]:
In [33]:
banknotes_train, banknotes_test = train_test_split(banknotes_data, test_size=50, random_state=1)
banknotes_train.shape
banknotes_test.shape
Out[33]:
Out[33]:
In [35]:
banknotes_data.columns
Out[35]:
In [36]:
features_1 = [u'X1', u'X2', u'X3']
clf_logreg_1 = LogisticRegression()
clf_logreg_1.fit(banknotes_train[features_1].values, banknotes_train[u'real'].values)
Out[36]:
In [58]:
pred_1 = clf_logreg_1.predict(banknotes_test[features_1].values)
print('Error rate pred1: %f' % (1 - accuracy_score(banknotes_test[u'real'].values, pred_1)))
err_1 = np.array( [1 if banknotes_test[u'real'].values[i] == pred_1[i] else 0 for i in range(len(pred_1))] )
err_1
Out[58]:
In [43]:
features_2 = [u'X4', u'X5', u'X6']
clf_logreg_2 = LogisticRegression()
clf_logreg_2.fit(banknotes_train[features_2].values, banknotes_train[u'real'].values)
Out[43]:
In [59]:
pred_2 = clf_logreg_2.predict(banknotes_test[features_2].values)
print('Error rate pred2: %f' % (1 - accuracy_score(banknotes_test[u'real'].values, pred_2)))
err_2 = np.array( [1 if banknotes_test[u'real'].values[i] == pred_2[i] else 0 for i in range(len(pred_2))] )
err_2
Out[59]:
In [49]:
def proportions_diff_confint_rel(sample1, sample2, alpha = 0.05):
z = stats.norm.ppf(1 - alpha / 2.)
sample = zip(sample1, sample2)
n = len(sample)
f = sum([1 if (x[0] == 1 and x[1] == 0) else 0 for x in sample])
g = sum([1 if (x[0] == 0 and x[1] == 1) else 0 for x in sample])
left_boundary = float(f - g) / n - z * np.sqrt(float((f + g)) / n**2 - float((f - g)**2) / n**3)
right_boundary = float(f - g) / n + z * np.sqrt(float((f + g)) / n**2 - float((f - g)**2) / n**3)
return (left_boundary, right_boundary)
In [50]:
def proportions_diff_z_stat_rel(sample1, sample2):
sample = zip(sample1, sample2)
n = len(sample)
f = sum([1 if (x[0] == 1 and x[1] == 0) else 0 for x in sample])
g = sum([1 if (x[0] == 0 and x[1] == 1) else 0 for x in sample])
return float(f - g) / np.sqrt(f + g - float((f - g)**2) / n )
In [64]:
print('95%% confidence interval for a difference between predictions: [%.4f, %.4f]' %
proportions_diff_confint_rel(err_1, err_2))
In [63]:
print('p-value: %f' % proportions_diff_z_test(proportions_diff_z_stat_rel(err_1, err_2)))
Ежегодно более 200000 людей по всему миру сдают стандартизированный экзамен GMAT при поступлении на программы MBA. Средний результат составляет 525 баллов, стандартное отклонение — 100 баллов.
Сто студентов закончили специальные подготовительные курсы и сдали экзамен. Средний полученный ими балл — 541.4. Проверьте гипотезу о неэффективности программы против односторонней альтернативы о том, что программа работает. Отвергается ли на уровне значимости 0.05 нулевая гипотеза? Введите достигаемый уровень значимости, округлённый до 4 знаков после десятичной точки.
In [70]:
avg_mark = 541.4
avg_mark_data = np.ones(100) * avg_mark
In [78]:
exp_mark = 525.
st_dev = 100.
In [81]:
def z_test(mean_val, exp_val, st_dev, num):
standard_error = st_dev / np.sqrt(num)
return (mean_val - exp_val) / standard_error
In [82]:
def p_val_greater(z_stat):
return 1 - stats.norm.cdf(z_stat)
In [90]:
z_stat = z_test(avg_mark, exp_mark, st_dev, (len(avg_mark_data)))
z_stat
Out[90]:
In [95]:
p_val = p_val_greater(z_stat)
p_val
Out[95]:
In [96]:
avg_mark2 = 541.5
In [97]:
z_stat2 = z_test(avg_mark2, exp_mark, st_dev, (len(avg_mark_data)))
z_stat2
Out[97]:
In [98]:
p_val2 = p_val_greater(z_stat2)
p_val2
Out[98]: