In [1]:
from __future__ import division
import numpy as np
import pandas as pd
from scipy import stats
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
Давайте вернёмся к данным выживаемости пациентов с лейкоцитарной лимфомой из видео про критерий знаков:
49,58,75,110,112,132,151,276,281,362∗
Измерено остаточное время жизни с момента начала наблюдения (в неделях); звёздочка обозначает цензурирование сверху — исследование длилось 7 лет, и остаточное время жизни одного пациента, который дожил до конца наблюдения, неизвестно.
Поскольку цензурировано только одно наблюдение, для проверки гипотезы H0:medX=200 на этих данных можно использовать критерий знаковых рангов — можно считать, что время дожития последнего пациента в точности равно 362, на ранг этого наблюдения это никак не повлияет.
Критерием знаковых рангов проверьте эту гипотезу против двусторонней альтернативы, введите достигаемый уровень значимости, округлённый до четырёх знаков после десятичной точки.
In [2]:
surv_data = np.array( [49, 58, 75, 110, 112, 132, 151, 276, 281, 362] )
surv_data
surv_data.mean()
Out[2]:
Out[2]:
In [3]:
exp_surv = 200
stats.wilcoxon(surv_data - exp_surv)
Out[3]:
В ходе исследования влияния лесозаготовки на биоразнообразие лесов острова Борнео собраны данные о количестве видов деревьев в 12 лесах, где вырубка не ведётся:
22,22,15,13,19,19,18,20,21,13,13,15,
и в 9 лесах, где идёт вырубка:
17,18,18,15,12,4,14,15,10.
Проверьте гипотезу о равенстве среднего количества видов в двух типах лесов против односторонней альтернативы о снижении биоразнообразия в вырубаемых лесах. Используйте ранговый критерий. Чему равен достигаемый уровень значимости? Округлите до четырёх знаков после десятичной точки.
In [4]:
forest_not_cut = np.array( [22, 22, 15, 13, 19, 19, 18, 20, 21, 13, 13, 15] )
forest_not_cut.mean()
forest_cut = np.array( [17, 18, 18, 15, 12, 4, 14, 15, 10] )
forest_cut.mean()
Out[4]:
Out[4]:
In [5]:
stats.mannwhitneyu(forest_not_cut, forest_cut, alternative='greater')
Out[5]:
28 января 1986 года космический шаттл "Челленджер" взорвался при взлёте. Семь астронавтов, находившихся на борту, погибли. В ходе расследования причин катастрофы основной версией была неполадка с резиновыми уплотнительными кольцами в соединении с ракетными ускорителями. Для 23 предшествовавших катастрофе полётов "Челленджера" известны температура воздуха и появление повреждений хотя бы у одного из уплотнительных колец.
challenger.txt
С помощью бутстрепа постройте 95% доверительный интервал для разности средних температур воздуха при запусках, когда уплотнительные кольца повреждались, и запусках, когда повреждений не было. Чему равна его ближайшая к нулю граница? Округлите до четырёх знаков после запятой.
Чтобы получить в точности такой же доверительный интервал, как у нас:
установите random seed = 0 перед первым вызовом функции get_bootstrap_samples, один раз сделайте по 1000 псевдовыборок из каждой выборки.
In [14]:
challenger = pd.read_csv('challenger.txt', delimiter='\t')
challenger.describe()
challenger.info()
challenger.head()
Out[14]:
Out[14]:
In [17]:
challenger.rename(columns={'Unnamed: 0': 'Date'}, inplace=True)
challenger.head()
Out[17]:
In [41]:
challenger_broken = challenger.loc[challenger.loc[:, 'Incident'] == 1, :].drop(['Date', 'Incident'], axis=1)
challenger_broken
challenger_not_broken = challenger.loc[challenger.loc[:, 'Incident'] != 1, :].drop(['Date', 'Incident'], axis=1)
challenger_not_broken
Out[41]:
Out[41]:
In [42]:
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
In [43]:
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
In [71]:
np.random.seed(0)
challenger_broken_bs_mean = np.array(map(np.mean, get_bootstrap_samples(challenger_broken['Temperature'].values, 1000)))
challenger_not_broken_bs_mean = np.array(map(np.mean, get_bootstrap_samples(challenger_not_broken['Temperature'].values, 1000)))
In [77]:
print('95%% confidence interval for times decrease of infarction: %s' %
str(stat_intervals(challenger_broken_bs_mean - challenger_not_broken_bs_mean, 0.05)))
На данных предыдущей задачи проверьте гипотезу об одинаковой средней температуре воздуха в дни, когда уплотнительный кольца повреждались, и дни, когда повреждений не было. Используйте перестановочный критерий и двустороннюю альтернативу. Чему равен достигаемый уровень значимости? Округлите до четырёх знаков после десятичной точки.
Чтобы получить такое же значение, как мы:
установите random seed = 0; возьмите 10000 перестановок.
In [84]:
def permutation_t_stat_ind(sample1, sample2):
return np.mean(sample1) - np.mean(sample2)
In [85]:
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]
In [86]:
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
In [87]:
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)
In [88]:
np.random.seed(0)
print('p-value: %.4f' % permutation_test(challenger_broken['Temperature'].values,
challenger_not_broken['Temperature'].values, max_permutations=10000))