В одном из выпусков программы "Разрушители легенд" проверялось, действительно ли заразительна зевота. В эксперименте участвовало 50 испытуемых, проходивших собеседование на программу. Каждый из них разговаривал с рекрутером; в конце 34 из 50 бесед рекрутер зевал. Затем испытуемых просили подождать решения рекрутера в соседней пустой комнате.
Во время ожидания 10 из 34 испытуемых экспериментальной группы и 4 из 16 испытуемых контрольной начали зевать. Таким образом, разница в доле зевающих людей в этих двух группах составила примерно 4.4%. Ведущие заключили, что миф о заразительности зевоты подтверждён.
Можно ли утверждать, что доли зевающих в контрольной и экспериментальной группах отличаются статистически значимо? Посчитайте достигаемый уровень значимости при альтернативе заразительности зевоты, округлите до четырёх знаков после десятичной точки.
In [1]:
import numpy as np
# import pandas as pd
# import math
# from scipy.stats import chisquare
# from statsmodels.stats.descriptivestats import sign_test
# from statsmodels.sandbox.stats.multicomp import multipletests
# import scipy
# import scipy as sc
# from statsmodels.stats.weightstats import *
# import pandas as pd
# import statsmodels.formula.api as smf
# import statsmodels.stats.api as sms
# import statsmodels.stats.multitest as smm
# import matplotlib.pyplot as plt
# %pylab inline
In [3]:
n = 50
n1 = 34
n2 = 16
p1 = 10.0/34
p2 = 4.0/16
print p1-p2
In [15]:
def proportions_diff_confint_ind(sample1, sample2, alpha = 0.05):
z = scipy.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 [16]:
alpha = 4.4/100
z = scipy.stats.norm.ppf(1 - alpha / 2.)
left_boundary = (p1 - p2) - z * np.sqrt(p1 * (1 - p1)/ 34 + p2 * (1 - p2)/ 16)
right_boundary = (p1 - p2) + z * np.sqrt(p1 * (1 - p1)/ 34 + p2 * (1 - p2)/ 16)
print (left_boundary, right_boundary)
In [18]:
P = float(p1*n1 + p2*n2) / (n1 + n2)
z_stat = (p1 - p2) / np.sqrt(P * (1 - P) * (1. / n1 + 1. / n2))
print 'greater', 1 - scipy.stats.norm.cdf(z_stat)
Ежегодно более 200000 людей по всему миру сдают стандартизированный экзамен GMAT при поступлении на программы MBA. Средний результат составляет 525 баллов, стандартное отклонение — 100 баллов.
Сто студентов закончили специальные подготовительные курсы и сдали экзамен. Средний полученный ими балл — 541.4. Проверьте гипотезу о неэффективности программы против односторонней альтернативы о том, что программа работает. Посчитайте достигаемый уровень значимости, округлите до 4 знаков после десятичной точки. Отвергается ли на уровне значимости 0.05 нулевая гипотеза?
In [28]:
import math
mu = 525
sigma = 100
n1 = 100
mu1 = 541.4
mu2 = 541.5
p_val_1 = 1 - stats.norm.cdf((mu1 - mu)/(sigma/math.sqrt(n1)))
p_val_2 = 1 - stats.norm.cdf((mu2 - mu)/(sigma/math.sqrt(n1)))
print p_val_1, p_val_2
In [2]:
import pandas as pd
df = pd.read_csv('banknotes.csv', delimiter='\t')
In [3]:
df.head()
Out[3]:
In [4]:
# # Split the Learning Set
from sklearn.cross_validation import train_test_split
target = df.real
df.drop(['real'], axis = 1, inplace = True)
X_fit, X_eval, y_fit, y_eval= train_test_split(
df, target, test_size=0.25, random_state=1
)
In [5]:
from sklearn import linear_model
from sklearn.metrics import accuracy_score
logreg = linear_model.LogisticRegression()
mod1 = logreg.fit(X_fit[['X1', 'X2', 'X3']], y_fit)
y1 = mod1.predict(X_eval[['X1', 'X2', 'X3']])
mod2 = logreg.fit(X_fit[['X4', 'X5', 'X6']], y_fit)
y2 = mod2.predict(X_eval[['X4', 'X5', 'X6']])
p1 = accuracy_score(y_eval, y1)
p2 = accuracy_score(y_eval, y2)
print p1, p2
In [6]:
delta1 = y_eval == y1
delta2 = y_eval == y2
In [9]:
def proportions_diff_confint_rel(sample1, sample2, alpha = 0.05):
z = scipy.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)
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 )
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 - scipy.stats.norm.cdf(np.abs(z_stat)))
if alternative == 'less':
return scipy.stats.norm.cdf(z_stat)
if alternative == 'greater':
return 1 - scipy.stats.norm.cdf(z_stat)
print proportions_diff_z_test(proportions_diff_z_stat_rel(delta1, delta2))
print proportions_diff_confint_rel(delta2, delta1)
In [71]:
print y1.size, X_fit.shape
In [19]:
long_life = np.array([49,58,75,110,112,132,151,276,281,362])
mu = 200
print stats.wilcoxon(long_life - mu), long_life.mean()
In [27]:
trees_no_cut = np.array([22,22,15,13,19,19,18,20,21,13,13,15])
trees_with_cut = np.array([17,18,18,15,12,4,14,15,10])
print stats.mannwhitneyu(trees_no_cut, trees_with_cut)
In [49]:
df = pd.read_csv('challenger.csv', delimiter='\t')
print df.head()
no_incident = df[df.Incident == 0].Temperature.as_matrix()
with_incident = df[df.Incident == 1].Temperature.as_matrix()
In [83]:
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
In [85]:
np.random.seed(0)
#get_bootstrap_samples(no_incident, 1000)
no_incident_scores = map(np.mean, get_bootstrap_samples(no_incident, 1000))
with_incident_scores = map(np.mean, get_bootstrap_samples(with_incident, 1000))
print 'no_incident_scores ', sum(no_incident_scores) / float(len(no_incident_scores))
print 'with_incident_scores ', sum(with_incident_scores) / float(len(with_incident_scores))
delta_median_scores = map(lambda x: x[0] - x[1], zip(no_incident_scores, with_incident_scores))
print 'delta ', sum(delta_median_scores) / float(len(delta_median_scores))
print "95% confidence interval for the difference between medians", stat_intervals(delta_median_scores, 0.05)
На данных предыдущей задачи проверьте гипотезу об одинаковой средней температуре воздуха в дни, когда уплотнительный кольца повреждались, и дни, когда повреждений не было. Используйте двустороннюю альтернативу. Чему равен достигаемый уровень значимости? Округлите до четырёх знаков после десятичной точки.
Чтобы получить такое же значение, как мы:
установите random seed = 0; возьмите 10000 перестановок.
In [93]:
df = pd.read_csv('challenger.csv', delimiter='\t')
print df.head()
no_incident = df[df.Incident == 0].Temperature.as_matrix()
with_incident = df[df.Incident == 1].Temperature.as_matrix()
def permutation_t_stat_ind(sample1, sample2):
return np.mean(sample1) - np.mean(sample2)
def get_random_combinations(n1, n2, max_combinations):
index = range(n1 + n2)
indices = set([tuple(index)])
for i in range(max_combinations - 1):
np.random.shuffle(index)
indices.add(tuple(index))
return [(index[:n1], index[n1:]) for index in indices]
def permutation_zero_dist_ind(sample1, sample2, max_combinations = None):
joined_sample = np.hstack((sample1, sample2))
n1 = len(sample1)
n = len(joined_sample)
if max_combinations:
indices = get_random_combinations(n1, len(sample2), max_combinations)
else:
indices = [(list(index), filter(lambda i: i not in index, range(n))) \
for index in itertools.combinations(range(n), n1)]
distr = [joined_sample[list(i[0])].mean() - joined_sample[list(i[1])].mean() \
for i in indices]
return distr
pylab.hist(permutation_zero_dist_ind(no_incident, with_incident, max_combinations = 10000))
pylab.show()
In [96]:
def permutation_test(sample, mean, max_permutations = None, alternative = 'two-sided'):
if alternative not in ('two-sided', 'less', 'greater'):
raise ValueError("alternative not recognized\n"
"should be 'two-sided', 'less' or 'greater'")
t_stat = permutation_t_stat_ind(sample, mean)
zero_distr = permutation_zero_dist_ind(sample, mean, max_permutations)
if alternative == 'two-sided':
return sum([1. if abs(x) >= abs(t_stat) else 0. for x in zero_distr]) / len(zero_distr)
if alternative == 'less':
return sum([1. if x <= t_stat else 0. for x in zero_distr]) / len(zero_distr)
if alternative == 'greater':
return sum([1. if x >= t_stat else 0. for x in zero_distr]) / len(zero_distr)
np.random.seed(0)
print "p-value: %f" % permutation_test(no_incident, with_incident, max_permutations = 10000)
In [47]:
df = pd.read_csv('water.csv', delimiter='\t')
print df.head()
In [48]:
#1 -0.6548
df.corr(method='pearson')
Out[48]:
In [49]:
#2 -0.6317
df.corr(method='spearman')
Out[49]:
In [9]:
df_north = df[df.location == 'North']
df_south = df[df.location == 'South']
In [50]:
df_north.corr(method='pearson')
Out[50]:
In [51]:
df_south.corr(method='pearson')
Out[51]:
In [ ]:
#3 min = -0.3686
In [18]:
#4 корреляция метьюса = 0.109
def mcc(a,b,c,d):
return (a*d-b*c)/(math.sqrt((a+b)*(a+c)*(b+d)*(c+d)))
print mcc(239, 515, 203, 718)
In [3]:
p_man1 = 239./(239+515)
p_man2 = 1 - p_man1
p_woman1 = 203./(203+718)
p_woman2 = 1 - p_woman1
p = [[p_man1, p_man2], [p_woman1, p_woman2]]
num1 = [[239,515],[203, 718]]
num2 = [[239,203],[515, 718]]
print scipy.stats.chi2_contingency(p,0), scipy.stats.chi2_contingency(num1),scipy.stats.chi2_contingency(num2)
Returns:
chi2 : float
The test statistic.
p : float
The p-value of the test
dof : int
Degrees of freedom
expected : ndarray, same shape as observed
The expected frequencies, based on the marginal sums of the table.
In [12]:
#6 0.0952
alpha = 0.05
z = scipy.stats.norm.ppf(1 - alpha / 2.)
p1 = p_man1
p2 = p_woman1
n1 = 239+515
n2 = 203+718
left_boundary = (p1 - p2) - z * np.sqrt(p1 * (1 - p1)/ n1 + p2 * (1 - p2)/ n2)
right_boundary = (p1 - p2) + z * np.sqrt(p1 * (1 - p1)/ n1 + p2 * (1 - p2)/ n2)
print (left_boundary, right_boundary)
In [5]:
P = float(p1*n1 + p2*n2) / (n1 + n2)
z_stat = (p1 - p2) / np.sqrt(P * (1 - P) * (1. / n1 + 1. / n2))
print 2 * (1 - scipy.stats.norm.cdf(np.abs(z_stat)))
In [6]:
# Не доволен Более или менее Доволен
# Не очень счастлив 197 111 33
# Достаточно счастлив 382 685 331
# Очень счастлив 110 342 333
pd_happiness = [[197,111,33],[382,685,331],[110,342,333]]
scipy.stats.chi2_contingency(pd_happiness)
Out[6]:
In [10]:
print math.sqrt((293.68311039689746)/((n1+n2)*(3-1)))
In [64]:
df = pd.read_csv('Aucs.csv', delimiter='\t')
df
Out[64]:
In [65]:
df['C4.5'].as_matrix()
Out[65]:
In [68]:
#2 C4.5 C4.5+m+cf
#3 2
#4 m
corr_data = []
for i, lhs_column in enumerate(df.columns):
for j, rhs_column in enumerate(df.columns):
if i >= j or (i == 0 or j ==0):
continue
print lhs_column
corr_data.append([lhs_column, rhs_column, stats.wilcoxon(df[lhs_column] - df[rhs_column]).pvalue])
corr_data
Out[68]:
In [69]:
#5 метод холма ans = 0
reject, p_corrected, a1, a2 = multipletests([row[2] for row in corr_data], alpha = 0.05, method = 'holm')
reject, p_corrected
Out[69]:
In [70]:
#5 метод Бенджамини-Хохберга ans = 0
reject, p_corrected, a1, a2 = multipletests([row[2] for row in corr_data], alpha = 0.05, method = 'fdr_bh')
reject, p_corrected
Out[70]:
In [175]:
df = pd.read_csv('botswana.tsv', delimiter='\t')
df.head()
Out[175]:
In [176]:
# 5ая - ответ 4
print df.religion.unique(), len(df.religion.unique()), df.shape
In [177]:
df.dropna().shape
Out[177]:
In [178]:
# 5ая - ответ 1834
df.isnull().sum()
Out[178]:
In [179]:
df['nevermarr'] = 0
df.nevermarr.mask(df.agefm.isnull(), 1, inplace = True)
df.drop(['evermarr'], axis = 1, inplace=True)
df.agefm.fillna(value = 0, inplace = True)
df.heduc.mask(df.nevermarr == 1,-1, inplace = True)
#7ая 123
df[df.heduc.isnull() == True].shape
Out[179]:
In [180]:
df['nevermarr'].sum()
Out[180]:
In [181]:
df['idlnchld_noans'] = 0
df.idlnchld_noans.mask(df.idlnchld.isnull(), 1, inplace = True)
df.idlnchld.fillna(value = -1, inplace = True)
df['heduc_noans'] = 0
df.heduc_noans.mask(df.heduc.isnull(), 1, inplace = True)
df.heduc.fillna(value = -2, inplace = True)
df['usemeth_noans'] = 0
df.usemeth_noans.mask(df.usemeth.isnull(), 1, inplace = True)
df.usemeth.fillna(value = -1, inplace = True)
In [182]:
#8 78264
print df.dropna().shape, df.dropna().shape[0]*df.dropna().shape[1]
df.dropna(inplace = True)
print df.shape
In [183]:
df.columns
Out[183]:
In [184]:
m1 = smf.ols('ceb ~ age + educ + religion + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + radio + tv + bicycle + nevermarr + idlnchld_noans + heduc_noans + usemeth_noans',
data=df)
fitted = m1.fit()
print fitted.summary()
#9ая 0.644
In [185]:
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
In [186]:
m4 = smf.ols('ceb ~ age + educ + religion + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + radio + tv + bicycle + nevermarr + idlnchld_noans + heduc_noans + usemeth_noans', data=df)
fitted = m4.fit(cov_type='HC1')
print fitted.summary()
plt.figure()
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
pylab.show()
In [187]:
m5 = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans + heduc_noans + usemeth_noans',
data=df)
fitted = m5.fit(cov_type='HC1')
print fitted.summary()
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
In [188]:
#12ая ответ 0.4672
print "F=%f, p=%f, k1=%f" % m4.fit(cov_type='HC1').compare_f_test(m5.fit(cov_type='HC1'))
In [190]:
m_meh = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + '\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans + heduc_noans',
data=df)
fitted = m_meh.fit(cov_type='HC1')
print fitted.summary()
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
print "F=%f, p=%f, k1=%f" % m4.fit(cov_type='HC1').compare_f_test(m_meh.fit(cov_type='HC1'))
print m4.fit(cov_type='HC1').compare_f_test(m_meh.fit(cov_type='HC1'))
#13ая 36
In [2]:
df = pd.read_csv('gene_high_throughput_sequencing.csv')
df.head()
Out[2]:
In [3]:
df.Diagnosis.unique()
Out[3]:
In [4]:
df.drop(['Patient_id'], axis = 1, inplace = True)
In [5]:
df_norm = df[df.Diagnosis == 'normal']
df_norm.drop(['Diagnosis'], axis = 1, inplace = True)
df_neop = df[df.Diagnosis == 'early neoplasia']
df_neop.drop(['Diagnosis'], axis = 1, inplace = True)
df_cans = df[df.Diagnosis == 'cancer']
df_cans.drop(['Diagnosis'], axis = 1, inplace = True)
In [6]:
df_neop.head()
Out[6]:
In [7]:
pvalues1 = scipy.stats.ttest_ind(df_norm, df_neop, equal_var = False).pvalue
pvalues2 = scipy.stats.ttest_ind(df_neop, df_cans, equal_var = False).pvalue
In [8]:
#1
alpha = 0.05
num1 = 0
num2 = 0
for i in pvalues1:
if i < alpha:
num1 = num1 + 1
for i in pvalues2:
if i < alpha:
num2 = num2 + 1
print num1, num2, num1 + num2
In [9]:
#2 fold change > 1.5.
import statsmodels.stats.multitest as smm
alpha = 0.05/2
reject1, p_corrected1, a1, a2 = smm.multipletests(pvalues1, alpha = alpha, method = 'holm')
reject2, p_corrected2, a1_2, a2_2 = smm.multipletests(pvalues2, alpha = alpha, method = 'holm')
# reject, p_corrected, a1, a2 = multipletests(sales_correlation.p,
# alpha = 0.05,
# method = 'holm')
In [15]:
def fold_change(C, T):
if T > C:
return 1.*T/C
if C > T:
return 1.*C/T
else:
return 0
fold_change(1,6)
Out[15]:
In [ ]:
In [34]:
num1 = 0
for id, p_corrected in enumerate(p_corrected1):
if (p_corrected < 0.025) and (fold_change(df_norm[[id]].mean(axis = 0)[0], df_neop[[id]].mean(axis = 0)[0]) > 1.5):
num1 = num1 + 1
print num1
In [35]:
num2 = 0
for id, val in enumerate(p_corrected2):
if (val < 0.025) and (fold_change(df_neop[[id]].mean(axis = 0)[0], df_cans[[id]].mean(axis = 0)[0]) > 1.5):
num2 = num2 + 1
print num2
In [20]:
fold_change(df_norm[[0]].mean(axis = 0)[0], df_neop[[0]].mean(axis = 0)[0])
Out[20]:
In [23]:
df_cans.shape
Out[23]:
In [36]:
alpha = 0.05/2
reject1, p_corrected1, a1, a2 = smm.multipletests(pvalues1, alpha = alpha, method = 'fdr_bh')
reject2, p_corrected2, a1_2, a2_2 = smm.multipletests(pvalues2, alpha = alpha, method = 'fdr_bh')
num1 = 0
for id, p_corrected in enumerate(p_corrected1):
if (p_corrected < 0.025) and (fold_change(df_norm[[id]].mean(axis = 0)[0], df_neop[[id]].mean(axis = 0)[0]) > 1.5):
num1 = num1 + 1
print num1
num2 = 0
for id, val in enumerate(p_corrected2):
if (val < 0.025) and (fold_change(df_neop[[id]].mean(axis = 0)[0], df_cans[[id]].mean(axis = 0)[0]) > 1.5):
num2 = num2 + 1
print num2
In [ ]:
state — штат США
account_length — длительность использования аккаунта
area_code — деление пользователей на псевдорегионы, использующееся в телекоме
intl_plan — подключена ли у пользователя услуга международного общения
vmail_plan — подключена ли у пользователя услуга голосовых сообщений
vmail_message — количество голосых сообщений, который пользователь отправил / принял
day_calls — сколько пользователь совершил дневных звонков
day_mins — сколько пользователь проговорил минут в течение дня
day_charge — сколько пользователь заплатил за свою дневную активность
eve_calls, eve_mins, eve_charge — аналогичные метрики относительно вечерней активности
night_calls, night_mins, night_charge — аналогичные метрики относительно ночной активности
intl_calls, intl_mins, intl_charge — аналогичные метрики относительно международного общения
custserv_calls — сколько раз пользователь позвонил в службу поддержки
treatment — номер стратегии, которая применялись для удержания абонентов (0, 2 = два разных типа воздействия, 1 = контрольная группа)
mes_estim — оценка интенсивности пользования интернет мессенджерами
churn — результат оттока: перестал ли абонент пользоваться услугами оператора
In [3]:
df = pd.read_csv('churn_analysis.csv')
df.head()
In [ ]: