In [ ]:
import pandas as pd
import scipy.stats as sps
import matplotlib.pyplot as plt 
import numpy as np
import seaborn as sns

import random
# зафиксируем сид для воспроизводимости генерации
SEED = 42
np.random.seed(SEED)
random.seed(SEED)

Критерии вида t-test¶

Одновыборочный¶

Дана одна нормальная выборка $X_1, ..., X_n \sim \mathcal{N}(a, \sigma^2)$.

Критерий проверяет гипотезы

$\mathsf{H}_0\colon a = a_0$

$\mathsf{H}_1\colon a \not= a_0$

ttest_1samp(a, popmean): statistic, pvalue

  • a — выборка
  • popmean — равно $a_0$

Сгенерируйте выборку $X_1, ..., X_n \sim \mathcal{N}(0, 1)$ и с помощью критерия проверьте:

  • равенство среднего нулю
  • равенство среднего 0.5
In [ ]:
size=100

sample = <...>
<...>

Двухвыборочный¶

Независимые выборки¶

Даны две независимые нормальные выборки

  • $X_1, ..., X_n \sim \mathcal{N}(a_1, \sigma_1^2)$,

  • $Y_1, ..., Y_m \sim \mathcal{N}(a_2, \sigma_2^2)$.

Критерий проверяет для их гипотезы о равенстве среднего:

$\mathsf{H}_0\colon a_1 = a_2$

$\mathsf{H}_1\colon a_1 \not= a_2$

ttest_ind(a, b, equal_var=True): statistic, pvalue

a, b — выборка

equal_var — известно ли равенство дисперсий

Сгенерируйте выборки $X_1, ..., X_n \sim \mathcal{N}(0, 1)$ и $X_1, ..., X_m \sim \mathcal{N}(1, 1)$. Используя критерий, проверьте равенство средних двух выборок.

In [ ]:
sample_1 = <...>
sample_2 = <...>
<...>

Сгенерируйте выборки $X_1, ..., X_n \sim \mathcal{N}(0, 1)$ и $X_1, ..., X_m \sim \mathcal{N}(1, 7)$. Используя критерий, проверьте равенство средних двух выборок.

In [ ]:
sample_1 = <...>
sample_2 = <...>
<...>

Связные выборки¶

Даны две связные нормальные выборки

  • $X_1, ..., X_n \sim \mathcal{N}(a_1, \sigma_1^2)$,

  • $Y_1, ..., Y_n \sim \mathcal{N}(a_2, \sigma_2^2)$.

Критерий проверяет для их гипотезы о равенстве среднего:

$\mathsf{H}_0\colon a_1 = a_2$

$\mathsf{H}_1\colon a_1 \not= a_2$

ttest_rel(a, b): statistic, pvalue

a, b — выборка

Сгенерируйте выборку $X_1, ..., X_n \sim \mathcal{N}(0, 1)$. Вторую сгенерируйте по формуле выборка_1 + случайный шум, случайный шум из $\mathcal{N}(0, 0.5)$. Используя критерий, проверьте гипотезу о равенстве среднего.

In [ ]:
sample_1 = <...>
sample_2 = <...>
<...>

Сгенерируйте выборку $X_1, ..., X_n \sim \mathcal{N}(0, 1)$. Вторую сгенерируйте по формуле выборка_1 + случайный шум, случайный шум из $\mathcal{N}(0.5, 0.5)$. Используя критерий, проверьте гипотезу о равенстве среднего.

In [ ]:
sample_1 = <...>
sample_2 = <...>
<...>

Пример: ирисы Фишера¶

Визуализация данных

In [ ]:
df = sns.load_dataset("iris")

g = sns.PairGrid(df, hue='species')
g.map_lower(sns.kdeplot, cmap ="Blues_d")
g.map_upper(plt.scatter)
g.map_diag(sns.kdeplot, lw=3);

Как выглядят данные

In [ ]:
df.head()

Виды ирисов

In [ ]:
np.unique(df.species)
In [ ]:
sps.ttest_ind(df[df.species == 'setosa'].sepal_length, 
              df[df.species == 'versicolor'].sepal_length,
              equal_var=False)
In [ ]:
sps.ttest_ind(df[df.species == 'virginica'].sepal_length, 
              df[df.species == 'versicolor'].sepal_length,
              equal_var=False)
In [ ]:
sps.ttest_ind(df[df.species == 'virginica'].sepal_width, 
              df[df.species == 'versicolor'].sepal_width,
              equal_var=False)

Замечание. Строго говоря, неоходима поправка на множественное тестирование гипотез.

Вывод

Множественная проверка гипотез¶

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

Сформулируем задачу на статистическом языке:

$X_1...X_{100}$ — выборка из распределения $Bern(p)$

$p=0.5$ отвечает случайному угадыванию.

Проверьте гипотезу: $\mathsf{H}_0 \colon p=0.5$ vs $\mathsf{H}_1 \colon p \neq 0.5$. Используйте критерий Вальда.

В качестве асимптотически нормальной оценки можно использовать $\widehat{p} = \overline{X}$ с асимптотической дисперсией $\sigma^2(p) = p (1 - p)$.

Выпишем состоятельную оценку дисперсии и статистику критерия Вальда:

$\widehat{\sigma} = \sqrt{\overline{X} (1 - \overline{X})}$, $W = \sqrt{n} \frac{\overline{X} - 0.5}{\sqrt{\overline{X} (1 - \overline{X})}}$

Оценим реальный уровень значимости для этого критерия при размере выборки равном 100. К чему он должен быть близок? Для скорости вычислений используйте количество выборок равное $10^3$.

In [ ]:
sample_size = 100
sample_count = 1000

theta = 0.5
In [ ]:
def wald_test(sample, theta, estimation_theta, estimation_sigma, alternative='two_sided'):
    """
    param sample: реализация выборки
    param theta: истинное значение параметра
    param estimation_theta: оценка параметра
    param estimation_sigma: оценка асимптотической дисперсии оценки estimation_sigma
    param alternative: вид альтернативной гипотезы, может принимать одно из значений 'two_sided', 'less', 'greater'

    return statistic
    return p_value
    """

    alpha = 0.05
    z = sps.norm.ppf(1 - alpha/2)
    n = len(sample)
    statistic = np.sqrt(n) * (estimation_theta - theta) / estimation_sigma

    if alternative == 'two_sided':
        p_value = sps.norm.sf(np.abs(statistic)) + sps.norm.cdf(-np.abs(statistic))
        conf_int = round(estimation_theta - z*estimation_sigma/np.sqrt(n), 4), round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4)


    elif alternative == 'less':
        p_value = sps.norm.cdf(statistic)
        conf_int = (-np.inf, round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4))

    
    elif alternative == 'greater':
        p_value = sps.norm.sf(statistic)
        conf_int = (round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4), np.inf)

    else:
        raise ValueError('alternative name is wrong')

    return statistic, p_value, conf_int

Оценим реальный уровень значимости

In [ ]:
sample = <...>

estimation_theta = <...>
estimation_sigma = <...>


counter = 0

for i in range(sample_count):
    _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])
    is_rejected = <...>
    if is_rejected:
        counter += 1
        
counter / sample_count

Теперь представим, что мы хотим проверить большое количество людей на экстрасенсорные способности с помощью данного критерия.

Проведите аналогичный эксперимент: сгенерируйте $10^3$ выборок размера $100$ для $100$ людей. Посчитайте, сколько раз из 1000 в вашем наборе из 100 выборок хотя бы для одной гипотеза будет отвергнута.

In [ ]:
sample_all = <...>

counter = 0
for sample in sample_all:
    estimation_theta = <...>
    estimation_sigma = <...>

    for i in range(100):
        _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])
        is_rejected = <...>
        if is_rejected:
            counter += 1
            break
In [ ]:
counter / 1000

Вывод: <...>

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

Чему равен этот уровень значимости, если одновременно проверяются n гипотез?

Ответ: <...>

Проведите предыдущий эксперимент с использованием корректной процедуры. Поскольку в реализованной выше функции $\alpha$ зафиксировано, используйте критерий отвержения гипотезы с помощью p-value.

In [ ]:
sample_all = <...>

counter = 0
for sample in sample_all:
    estimation_theta = <...>
    estimation_sigma = <...>

    for i in range(100):
        _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])
        is_rejected = <...>
        if is_rejected:
            counter += 1
            break
In [ ]:
counter / 1000

Вывод: <...>