Проверка значимости или гипотез с помощью Python

"Начиная"

Проверка значимости или гипотез с помощью Python

Статистика колледжа с Python

Вступление

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

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

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

Опубликованные статьи:

Как обычно, код доступен на моем GitHub.

Введение в проверку гипотез

Антониу был на 30-дневном отпуске с группой из 4 друзей, и они хотели решить, кто будет готовить ужин каждый вечер. Антониу предложил, чтобы они могли записать свои имена и хранить их в миске. Затем он случайным образом выбирал имя из миски, и имя, которое появилось, должно было приготовить ужин в тот вечер. Теперь предположим, что после 4 раундов случайного выбора имен имя Антониу так и не появилось. Остальные друзья начали подозревать, что что-то происходит. Давайте вычислим эту вероятность.

0.80**4
0.4096000000000001

Мы видим, что вероятность реализации этого сценария составляет около 41%. Это очень большая вероятность, поэтому группе друзей, вероятно, следует доверять Антониу. А теперь представьте, что после 15 раундов Антониу все еще не выбирают.

0.8**15
0.03518437208883203

Если предположить, что процесс полностью случайный, вероятность того, что его не выберут 15 раундов подряд, составляет около 3,5%. Если бы его друзья использовали часто определяемый порог в 5%, они не ожидали бы, что событие, которое произошло чисто случайно, имеет вероятность менее 5%. Тогда они, вероятно, отвергли бы гипотезу о том, что это был действительно случайный процесс, что означало бы, что Антониу жульничал.

Проверка гипотез о доле населения

Однажды на пляже Антониу прочитал статью в журнале, в которой говорилось, что только 30% людей пользуются солнцезащитным кремом на пляже. Антониу подумал, что эта цифра слишком мала для пляжной зоны, где они отдыхали. Он рассказал своим друзьям, и они решили случайным образом выбрать 35 человек с близлежащих пляжей, где, по их оценкам, в тот день на пляже наслаждались не менее 400 человек. Из 35 человек 14 пользовались солнцезащитным кремом.

14/35
0.4

Давайте определим нашу проверку гипотез более формально:

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

Давайте проведем моделирование, чтобы увидеть, насколько вероятно, что выборка, подобная приведенной выше, возникнет случайно. Мы моделируем 40 образцов размером 30, из которых 40% людей использовали солнцезащитный крем.

from scipy.stats import bernoulli, norm, t
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
x = np.zeros((40, 35))
for i in range(40):    
    x[i] = bernoulli.rvs(p=0.3, size=35)

Мы можем построить наше моделирование с помощью гистограммы.

sns.histplot(np.mean(x, axis=1));

p = np.mean(x, axis=1)

Обладая этой информацией, мы можем вычислить наше p-значение. Его можно определить как вероятность того, что p̂ будет больше или равна выборке 0,4, при условии, что H_0 истинно:

Мы можем легко вычислить p-значение.

p[p>=0.40].shape[0]/p.shape[0]
0.125

Если Антониу использует уровень значимости α = 0,05, он не откажется от H_0. Это означает, что при p̂ = 0,4 у них недостаточно доказательств, чтобы отвергнуть возможность того, что только 30% людей используют солнцезащитный крем, а статистика, полученная ими в 40%, была получена исключительно случайно.

if (p[p>=0.4].shape[0]/p.shape[0]<0.05):
    print('Reject H_0')
else:
    print('Do not reject H_0')
Do not reject H_0

Условия Z-теста на пропорцию

В предыдущих статьях мы видели условия для выполнения z-теста. Три условия:

Выборка должна быть случайной.
Нормальное распределение может приблизительно соответствовать выборочному распределению пропорций выборки. Эмпирическое правило гласит, что у вас должно быть не менее 10 успехов и 10 неудач.
Образцы должны быть независимыми. Эмпирическое правило состоит в том, что если вы отбираете выборку без замены, размер вашей выборки должен быть менее 10% от размера генеральной совокупности.

Давайте проработаем наш предыдущий пример и проверим каждое из условий. Начиная со случайной выборки, нам сказали, что группа друзей произвольно выбирала людей, чтобы задать вопросы. Имея эту информацию, мы собираемся предположить, что их процесс был действительно случайным. Нормальное условие говорит нам, что ожидаемое количество успехов и неудач должно быть как минимум равным 10. Мы также смогли соответствовать этим критериям (см. Расчеты ниже). Наконец, у нас есть условие независимости. В нашем случае численность населения составляла около 400 человек, проживающих на пляже. Размер выборки составлял 35, поэтому мы также могли выполнить условие для выполнения z-теста.

n = 35
p = 0.3

print('np = '+ str(n*p))
print('n(1-p) = '+ str(n*(1-p)))
np = 10.5
n(1-p) = 24.5

Вычисление статистики z в тесте на пропорцию

Мы продолжим пользоваться солнцезащитным кремом. В предыдущем разделе мы проверили, что все условия для вывода были выполнены. В этом случае мы собираемся выполнить z-тест.

Прежде чем продолжить, напомним нашу цель. Из всего населения пляжной зоны Антониу и его друзья взяли выборку из 35 человек и вычислили статистику выборки ̂p = 40%. Теперь, если предположить, что нулевая гипотеза верна, какова вероятность получить результат так далеко или дальше от предполагаемой доли населения p = 30%.

Один из способов рассчитать эту вероятность - вычислить, на сколько стандартных отклонений отстоит выбранная пропорция от истинной пропорции (помните, что мы всегда предполагаем, что H_0 истинно, поэтому мы считаем эту пропорцию истинной).

p_0 = 0.3
p_hat = 0.4

Z = (0.4-0.3)/np.sqrt((0.3*(1-0.3)/n))
Z
1.290994448735806

Тогда наше p-значение может быть определено как:

1-norm.cdf(Z)
0.09835280122947343

Сравнивая полученную нами вероятность с уровнем значимости α = 0,05, мы снова видим, что мы не можем отклонить H_0.

Проверка гипотез о среднем населении

Однажды ночью Антониу проверял свой телефон и начал задаваться вопросом, сколько сообщений он и его друзья обменивались в групповом чате WhatsApp, который они создали за 3 месяца до поездки в отпуск. Он подозревал, что в среднем они отправляли друг другу 180 сообщений в день, но один из его друзей утверждал, что они обменивались более чем 180 сообщениями. Они решили проверить теорию.

Давайте определимся с нашим критерием значимости:

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

Как мы видели ранее, в этом случае мы будем выполнять t-тест, а не z-тест. Причина в том, что нам неизвестно стандартное отклонение выборочного распределения выборочного среднего.

Мы можем только оценить это, используя стандартное отклонение выборки

Затем мы можем написать наш тест как:

x = [266, 190, 161, 191, 241, 180, 240, 251]
x = np.asarray(x)
sns.histplot(x, bins=10);

print(x.mean())
print(x.std())
215.0
36.29738282576307
μ = 180
n = 8

t_star = (x.mean()-μ)/(x.std()/np.sqrt(n))
t_star
2.727330227672841
1-t.cdf(t_star, df=n-1)
0.01472601649385541

Наше значение p составляет примерно 0,015, что меньше нашего уровня значимости. Это означает, что мы можем отклонить H_0, выдвинутый Антониу, и принять альтернативную гипотезу. Есть достаточно доказательств, чтобы принять гипотезу друга Антониу о том, что группа отправляла более 180 сообщений в день.

Условия для t-теста на среднее значение по совокупности

В предыдущих статьях мы видели условия для проведения t-теста. Три условия:

Выборка должна быть случайной.
Выборочное распределение выборочного среднего может быть аппроксимировано нормальным распределением. Этого можно добиться тремя способами: размер выборки больше 30 (применяется центральная предельная теорема), исходное распределение является нормальным или распределение выборки является симметричным и не показывает выбросов.
Образцы должны быть независимыми. Эмпирическое правило состоит в том, что если вы отбираете выборку без замены, размер вашей выборки должен быть менее 10% от размера генеральной совокупности.

Еще раз, давайте проверим нашу предыдущую решенную задачу для каждого из условий. Это придаст нам больше уверенности в том, что результаты не имеют какой-либо предвзятости. Нам сказали, что Антониу и его друг выбрали дни случайным образом, что соответствует первому критерию. Во-вторых, наше выборочное распределение примерно симметрично, поэтому мы можем считать, что нормальное условие также выполняется. Наконец, Антониу и его друг использовали размер выборки 8. Нам сказали, что они настроили групповой чат WhatsApp за 3 месяца до отъезда. Размер выборки составляет менее 10%, поэтому можно с уверенностью предположить независимость между выбранными днями, даже если они были выбраны без замены.

Ошибки типа I и типа II

Мы уже видели, что если p-значение меньше уровня значимости α, мы должны отклонить H_0. И наоборот, если значение p больше, чем α, мы не должны отклонять H_0. Когда мы проводим эти сравнения, всегда есть вероятность того, что мы делаем ошибку. Чтобы было понятнее, нарисуем таблицу.

Обратите внимание, что вероятность получения ошибки типа I равна уровню значимости. Представьте, что вы определяете уровень значимости 5%. Это означает, что в 5% случаев, даже если ваша нулевая гипотеза верна, вы получите статистику, которая заставит вас отвергнуть нулевую гипотезу.

Какими будут ошибки типа I и типа II в приведенном выше примере? Давайте еще раз запишем нашу проверку гипотез о доле людей, пользующихся солнцезащитным кремом:

Ошибка типа I произойдет, если группа друзей придет к выводу, что доля людей, пользующихся солнцезащитным кремом, не равна 0,3, хотя это действительно так. С другой стороны, ошибка типа II будет иметь место, если доля людей не равна 0,3, но они не могут ее отвергнуть.

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

Обратите внимание, что

на самом деле является вероятностью совершения ошибки типа II, поэтому вы можете думать о мощности критерия значимости как о вероятности того, что не будет ошибка типа II.

Давайте посмотрим, что мы можем сделать, чтобы повысить эффективность наших проверок гипотез:

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

Заключение

В этой статье мы рассмотрели проверку значимости или проверку гипотез. Они чрезвычайно полезны для использования выборочных данных для оценки вероятности утверждения о значении генеральной совокупности. Проверяемое значение может быть пропорцией или средним значением для генеральной совокупности. Мы рассчитали p-значения для обоих случаев, используя z-тест и t-тест, соответственно.

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

Упражнения

Вы получите решения в статье на следующей неделе.

Согласно большому опросу, проведенному в прошлом году, около 85% домов в Порту имеют доступ к высокоскоростному Интернету. Марко поинтересовался, изменилась ли пропорция, и произвел случайную выборку из 80 домов и обнаружил, что 75 имеют доступ к высокоскоростному Интернету. Он хочет использовать эти образцы данных, чтобы проверить, действительно ли соотношение изменилось. Предполагая, что условия для вывода были выполнены, какой вы можете сделать вывод о доле домов с высокоскоростным Интернетом с учетом уровня значимости 1%?
Марта владеет фруктовым магазином и еженедельно получает арбузы. Поставщик заявляет, что они должны весить 1 кг. Марта решает взвесить случайную выборку из 100 арбузов и находит средний вес 850 г и стандартное отклонение 200 г. Она хочет использовать этот образец данных, чтобы проверить, меньше ли среднее значение, чем заявленное поставщиком, и пересмотреть договор, если это так. Если предположить, что условия для вывода были выполнены, что должна сделать Марта (рассмотрим уровень значимости 5%)?

Ответы за прошлую неделю

Измените функцию confidence_interval_mean_t и постройте 5 различных экспериментов, чтобы вычислить 95% интервал для выборочного среднего, в 3 из них условия для интервала t не выполняются, а в 2 - где они есть. Для трех случаев, когда условия не выполняются, определите следующие сценарии: выборка не случайна, исходное распределение не является приблизительно нормальным и независимость не соблюдается. Для двух случаев, когда выполняются условия, определите один случай, когда исходное распределение является нормальным, и другой, когда исходное распределение искажено, но среднее значение выборки имеет нормальное распределение. Истинное среднее значение генеральной совокупности должно содержаться в рассчитанных доверительных интервалах в 95% случаев для этих двух последних случаев.

Подсказка: вам может быть полезно использовать функцию skewnorm из scipy. Ниже представлена модифицированная версия нормального распределения, искаженная параметром асимметрии, средним значением и стандартным отклонением.

# code adapted from https://stackoverflow.com/questions/49367436/scipy-skewnorm-mean-not-matching-theory

skew = 4.0
mean = 2
stdev = 0.5

delta = skew / math.sqrt(1. + math.pow(skew, 2.))
adjStdev = math.sqrt(math.pow(stdev, 2.) / (1. - 2. * math.pow(delta, 2.) / math.pi))
adjMean = mean - adjStdev * math.sqrt(2. / math.pi) * delta

print('target mean={:.4f} actual mean={:.4f}'.format(mean, float(skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[0])))
print('target stdev={:.4f} actual stdev={:.4f}'.format(stdev, math.sqrt(float(skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[1]))))
target mean=2.0000 actual mean=2.0000
target stdev=0.5000 actual stdev=0.5000
# Original skewed distribution

plt.hist(skewnorm.rvs(a = skew, loc=adjMean, scale=adjStdev, size=2000), bins=50);

# Approximately normal distribution of the sample mean because sample 
# size is bigger than 30 (CTL applies)

plt.hist(np.mean([skewnorm.rvs(a = skew, loc=adjMean, scale=adjStdev, size=35) for _ in range(2000)], axis=1), bins=50);

def confidence_interval_mean_t(μ, σ, n, number_trials, N, ci=0.95, sample='random', dist='normal'):
    
    skew = 4.0
    mean = μ
    stdev = σ

    delta = skew / math.sqrt(1. + math.pow(skew, 2.))
    adjStdev = math.sqrt(math.pow(stdev, 2.) / (1. - 2. * math.pow(delta, 2.) / math.pi))
    adjMean = mean - adjStdev * math.sqrt(2. / math.pi) * delta
    
    if dist!='normal':
        x_ = skewnorm.rvs(a = 4, loc=adjMean, scale=adjStdev, size=N)
    else:
        x_ = norm.rvs(loc=μ, scale=σ, size=N)

    x_hat_list = []
    SE_hat_x_list = []
    
    if sample!='random':
        # Inducing bias on the sampling
        x_.sort()
        x_ = x_[:-int(0.2*N)]
        np.random.shuffle(x_)
    
    for i in range(number_trials):
        s_ = np.random.choice(x_, n, replace=False)
        x_hat = np.mean(s_)
        x_hat_list.append(x_hat)
        SE_hat_x_list.append(t.ppf(ci+(1-ci)/2, df=n-1)*np.std(s_)/np.sqrt(n))
    
    plt.hist(x_hat_list, bins=50)

    j=0
    _, ax = plt.subplots(1, 1, figsize=(6, 8))
    for i in range(len(x_hat_list)):
        if (μ>x_hat_list[i]-SE_hat_x_list[i]) & (μ<x_hat_list[i]+SE_hat_x_list[i]):
            # interval contains p
            if i > len(x_hat_list)-50:
                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='black')
            j +=1
        else:
            # interval does not contain p
            if i > len(x_hat_list)-50:    
                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='red')
    ax.axvline(μ, color='darkorange')
    #plt.xlim(0,1)
    plt.show()
    print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')
# Not random sampling
confidence_interval_mean_t(μ=2, σ=0.5, n=35, number_trials=1000, N=500, sample='not_random')

191/1000=0.19
# not approximately normal dist
confidence_interval_mean_t(μ=2, σ=0.5, n=15, number_trials=1000, N=500, dist='not normal')

944/1000=0.94
# not independent samples
confidence_interval_mean_t(μ=2, σ=0.5, n=150, number_trials=1000, N=600, dist='not normal')

973/1000=0.97
# Not normal distribution but n >= 30, so CTL applies; all other conditions are met
confidence_interval_mean_t(μ=2, σ=0.5, n=35, number_trials=1000, N=500, dist='not normal')

947/1000=0.95
# Original distribution is normal; all other conditions are met
confidence_interval_mean_t(μ=2, σ=0.5, n=35, number_trials=1000, N=500)

930/1000=0.93

Проверка значимости или гипотез с помощью Python

"Начиная"