Python стандартная ошибка среднего

— Я ничего в ней не вижу, — сказал я, возвращая шляпу Шерлоку Холмсу.

— Нет, Уотсон, видите, но не даете себе труда поразмыслить над тем, что видите.

Артур Конан Дойл. Голубой карбункул

В предыдущей серии постов для начинающих (первый пост тут) из ремикса книги Генри Гарнера «Clojure для исследования данных» (Clojure for Data Science) на языке Python было представлено несколько численных и визуальных подходов, чтобы понять, что из себя представляет нормальное распределение. Мы обсудили несколько описательных статистик, таких как среднее значение и стандартное отклонение, и то, как они могут использоваться для краткого резюмирования больших объемов данных.

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

В этой серии из 4-х постов мы рассмотрим статистический вывод: каким образом можно выйти за пределы простого описания выборок и вместо этого описать популяцию, из которой они были отобраны. Мы подробно рассмотрим степени нашей уверенности в выводах, которые мы делаем из выборочных данных. Мы раскроем суть робастного подхода к решению задач в области исследования данных, каким является проверка статистических гипотез, которая как раз и привносит научность в исcледование данных.

Кроме того, по ходу изложения будут выделены болевые точки, связанные с терминологическим дрейфом в отечественной статистике, иногда затуманивающим смысл и подменяющим понятия. В конце заключительного поста можно будет проголосовать за или против размещения следующей серии постов. А пока же…

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

Представляем AcmeContent

Для оказания помощи в иллюстрировании понятий, представленных в этой серии постов, предположим, что в компанию AcmeContent нас недавно назначили в качестве исследователей данных. Компания заведует веб-сайтом, предлагающим своим посетителям возможность делиться между собой понравившимися им видеоклипами по Интернет.

Одна из метрик, которая отслеживается в AcmeContent посредством веб-аналитики — это время пребывания. Указанная метрика служит мерой количества времени, в течение которого посетитель остается на веб-сайте. Безусловно, посетителям, которые проводят на веб-сайте продолжительное время, веб-сайт нравится, и в AcmeContent хотели бы, чтобы посетители оставались на нем максимально долго.

^{Время пребывания (dwell time)— это отрезок времени между временем, прибытия посетителя на веб-сайт и временем, когда он сделал последний запрос.}

^{Отскок (bounce) — это посетитель, который выполняет всего один запрос — его время пребывания равно нулю.}

На вас, как новом исследователе данных, лежит в компании обязанность анализировать время пребывания посетителей на веб-сайте — этот показатель фигурирует в аналитических отчетах посещаемости веб-страниц и разделов веб-сайта — и измерять успех веб-сайта AcmeContent.

Загрузка и обследование данных

Здесь мы будем пользоваться теми же самыми библиотеками, что и ранее: scipy, pandas и matplotlib. В предыдущей серии постов мы использовали библиотеку pandas для загрузки электронных таблиц Excel, задействуя ее функцию read_excel. Здесь мы будем загружать набор данных из текстового файла с разделением значений символом табуляции. Для этого мы воспользуемся функцией pandas read_csv, которая на входе ожидает URL-адрес либо путь к файлу в строковом формате.

Файл был любезно переформатирован веб-командой AcmeContent и содержит всего два столбца — дату запроса и время пребывания на веб-сайте в секундах. Заголовки столбцов расположены в первой строке файла:

^{Названия приводимых примеров имеют формат ex_N_M, где ex — example (пример), N — номер серии постов и M — порядковый номер в посте. Примеры оформлены в виде функций без аргументов и возвращаемых значений. Это сделано намеренно, т.к. задачно-ориентированный стиль изложения требует кратких и четких примеров без отвлекающей внимание информации. К тому же, в таком виде примеры могут быть собраны вместе и исполняться независимо в рамках программной оболочки.}

def load_data( fname ):
    return pd.read_csv('data/ch02/' + fname, 't')

def ex_2_1():
    return load_data('dwell-times.tsv').head()

Если выполнить этот пример (в консоли интерпретатора Python либо в блокноте Jupyter), то можно увидеть результат, который показан ниже:

Посмотрим, как выглядит время пребывания на гистограмме.

Визуализация времени пребывания

Мы можем построить гистограмму времени пребывания, выбрав столбец dwell-time и применив к нему функцию hist:

def ex_2_2():
    load_data('dwell-times.tsv')['dwell-time'].hist(bins=50)
    plt.xlabel('Время пребывания, сек.')
    plt.ylabel('Частота')
    plt.show()   

Приведенный выше пример сгенерирует следующую ниже гистограмму:

Очевидно, что эти данные не являются нормально распределенными; они даже не представляют собой сильно смещенное нормальное распределение. Слева от пика нет никакого хвоста (посетитель явно не может быть на веб-сайте менее 0 сек.). Сначала данные круто убывают вправо и потом тянутся вдоль оси X намного дольше, чем можно было бы ожидать от данных, которые нормально распределены.

Когда встречаются подобные распределения, где большинство значений малы, и иногда случаются предельные значения, может оказаться полезным изобразить ось Y в виде логарифмической шкалы. Логарифмические шкалы используются для того, чтобы представлять события, которые меняются в очень широких пределах. Оси графика обычно линейны, и они делят диапазон на равноразмерные шаги подобно «числовой оси», которую мы изучали в школе. В отличие от них, логарифмические шкалы делят диапазон на шаги, которые становятся все длиннее и длиннее по мере удаления от источника.

Некоторые системы измерения природных явлений, которые изменяются в очень широких пределах, представлены в логарифмической шкале. Например, шкала Рихтера для измерения интенсивности землетрясений является логарифмической шкалой по основанию 10, что означает, что землетрясение магнитудой 5 баллов по шкале Рихтера в 10 раз интенсивнее землетрясения магнитудой 4 балла. Децибельная шкала громкости тоже имеет логарифмическую шкалу, но с другим основанием — звуковая волна магнитудой 30 децибелов в 10 раз больше звуковой волны в 20 децибелов. В каждом случае принцип один и тот же — использование логарифмической шкалы позволяет сжать очень широкий предел значений в диапазон гораздо меньшего размера.

Изобразить на графике ось Y с логарифмической шкалой достаточно легко при помощи именованного аргумента logy=True функции pandas plot.hist:

def ex_2_3():
    load_data('dwell-times.tsv')['dwell-time'].plot.hist(bins=20, logy=True)
    plt.xlabel('Время пребывания, сек.')
    plt.ylabel('Логарифмическая частота')
    plt.show()  

В библиотеке pandas по умолчанию используется десятичная логарифмическая шкала, каждое деление на оси которой представляет собой интервал в 10 раз шире относительно предыдущего. График, в котором только одна ось имеет логарифмическую шкалу, называется полулогарифмическим или лог-линейным. Неудивительно, что график с двумя логарифмическими осями называется просто логарифмическим или логлог-линейным графиком (иногда также двойным логарифмическим, loglog=True).

Нанесение данных времени пребывания на полулогарифмический график показывает, что они обладают скрытой связностью — имеется линейная связь между временем пребывания и логарифмом частоты. На графике ясность связи ухудшается справа, где число посетителей меньше 10, но, помимо этого, связь удивительно постоянная.

Прямая линия на полулогарифмическом графике — это явный индикатор экспоненциального распределения.

Экспоненциальное распределение

Экспоненциальное (или показательное) распределение часто встречается в ситуациях, когда имеется много малых положительных значений и намного меньше более крупных. С учетом того, что мы узнали о шкале Рихтера, не будет никаким секретом, что магнитуда землетрясений подчинена экспоненциальному распределению.

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

Экспоненциальное распределение характерно многими интересными свойствами. Одно из них имеет отношение к среднему значению и стандартному отклонению:

def ex_2_4():
    ts = load_data('dwell-times.tsv')['dwell-time']
    print('Среднее:               ', ts.mean())    
    print('Медиана:               ', ts.median())
    print('Стандартное отклонение:', ts.std())

Среднее:                 93.2014074074074
Медиана:                 64.0
Стандартное отклонение:  93.96972402519819

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

^{Для экспоненциальных распределений средние значения и стандартные отклонения эквивалентны.}

Вторым свойством экспоненциального распределения является отсутствие памяти (или последействия). Это противоречащее интуитивному пониманию свойство лучше всего проиллюстрировать на примере. Мы ожидаем, что пока посетитель продолжает просматривать наш веб-сайт, вероятность, что ему надоест, и он покинет веб-сайт, увеличивается. Поскольку среднее время пребывания составляет 93 сек., то может создаться впечатление, что за пределами этих 93 сек., он будет продолжать просмотр веб-сайта все в меньшей степени.

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

Для распределения без памяти количество истекшего времени не влияет на вероятность продолжения события в течение дополнительных x минут.

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

Поскольку медианное время пребывания составляет 64 сек., то примерно половина наших посетителей остаются на веб-сайте в течение примерно всего одной минуты. Среднее значение в 93 сек. показывает, что некоторые посетители остаются намного дольше медианы. Эти статистики были вычислены по всем посетителям за последние 6 месяцев. Было бы интересно узнать, каким образом эти статистики варьируются в расчете на один день. Давайте их вычислим.

Распределение среднесуточных значений

Предоставленный веб-командой файл включает метку даты посещения. В целях агрегирования данных по дням, необходимо удалить ту часть метки даты, которая соответствует времени. Хотя мы можем проделать это при помощи нативных средств Python, более гибкий подход состоит в использовании функционала pandas для обработки даты и времени — функции to_datetime.

На входе она принимает строковое значение, значение с типом date-time, список, кортеж, 1-мерный массив либо числовой ряд Series библиотеки pandas как в нашем случае, а также набор именованных аргументов. Например, именованный аргумент errors='ignore' позволяет проигнорировать даты, которые неправильно отформатированы либо выходят за допустимые пределы. Отметим также, что в функции mean_dwell_times_by_date использован вспомогательный метод resample для частотных преобразований и перестановки временного ряда. Он выполняет группировку по дате-времени, а после него следует метод свертки по каждой группе. В данном случае аргумент 'D' группирует по будним дням, и затем каждая группа агрегируется методом mean. Таким образом, выражение dt.resample('D').mean() берет средние значения по будним дням:

def with_parsed_date(df):
    '''Привести поле date к типу date-time'''
    df['date'] = pd.to_datetime(df['date'], errors='ignore')
    return df

def filter_weekdays(df): 
    '''Отфильтровать по будним дням'''
    return df[df['date'].index.dayofweek < 5]  # понедельник..пятница

def mean_dwell_times_by_date(df):
    '''Среднесуточные времена пребывания'''
    df.index = with_parsed_date(df)['date']
    return df.resample('D').mean()  # перегруппировать  

def daily_mean_dwell_times(df):
    '''Средние времена пребывания с фильтрацией - только по будним дням'''
    df.index = with_parsed_date(df)['date']
    df = filter_weekdays(df)
    return df.resample('D').mean()

Сочетание приведенных выше функций позволяет вычислить среднее, медиану и стандартное отклонение для времен пребывания по будним дням:

def ex_2_5():
    df  = load_data('dwell-times.tsv')    
    mus = daily_mean_dwell_times(df)
    print('Среднее:                ', float(means.mean()))    
    print('Медиана:                ', float(means.median()))
    print('Стандартное отклонение: ', float(means.std()))

Среднее:                 90.21042865056198
Медиана:                 90.13661202185793
Стандартное отклонение:  3.7223429053200348

По будним дням среднее значение составляет 90.2 сек. Оно близко к среднему значению, которое мы вычислили ранее по всему набору данных, включая выходные дни. А вот стандартное отклонение намного ниже, всего 3.7 сек. Другими словами, распределение значений по будним дням имеет стандартные отклонения намного ниже, чем по всему набору данных. Давайте построим гистограмму значений времен пребывания по будним дням:

def ex_2_6():
    df = load_data('dwell-times.tsv')
    daily_mean_dwell_times(df)['dwell-time'].hist(bins=20)
    plt.xlabel('Время пребывания по будним дням, сек.')
    plt.ylabel('Частота')
    plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

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

Центральная предельная теорема

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

^{Согласно центральной предельной теореме, распределение выборочных средних нормально распределено независимо от распределения, из которого они были вычислены.}

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

В целях выведения нормальной кривой поверх гистограммы последняя должна быть построена как гистограмма плотностей распределения. Она изображает не частоту, а долю всех точек, помещенных в каждую корзину. После чего мы сможем наложить нормальную плотность вероятности с тем же средним значением и стандартным отклонением (Обратите внимание на применение функции dropna, которая удаляет строки, где есть пропущенные значения):

def ex_2_7():
    '''Подгонка нормальной кривой поверх гистограммы'''
    df = load_data('dwell-times.tsv')
    means = daily_mean_dwell_times(df)['dwell-time'].dropna() 
    ax = means.hist(bins=20, density=True)
    xs = sorted(means)    # корзины
    df = pd.DataFrame()
    df[0] = xs
    df[1] = stats.norm.pdf(xs, means.mean(), means.std())
    df.plot(0, 1, linewidth=2, color='r', legend=None, ax=ax)
    plt.xlabel('Время пребывания по будним дням, сек.')
    plt.ylabel('Плотность')    
    plt.show()

Этот пример сгенерирует следующий ниже график:

Нормальная кривая, изображенная поверх гистограммы, имеет стандартное отклонение, равное примерно 3.7 сек. Другими словами, она количественно определяет изменчивость каждого значения в будние дни относительно общего среднего значения, равного 90 сек. Можно представить среднее каждого дня как выборку из общей популяции, где изображенная выше кривая показывает распределение средних по выборкам. Поскольку 3.7 сек. — это величина, на которую среднее значение выборки отличается от популяционного среднего, она называется стандартной ошибкой.

Стандартная ошибка

В то время как стандартное отклонение измеряет величину изменчивости внутри выборки, стандартная ошибка (Standard Error, аббр. SE) среднего, измеряет величину изменчивости между средними значениями выборок, взятыми из той же самой популяции.

^{Стандартная ошибка среднего — это стандартное отклонение распределения средних по выборкам.}

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

Здесь σ_x — это стандартное отклонение, подсчитанное по выборке x, и n — размер выборки. Эта формула не похожа на описательные статистики, которые мы встречали в предыдущей серии постов. Они описывали одиночную выборку. В отличие от них, стандартная ошибка пытается описать свойство выборок в целом — величину изменчивости в средних значениях выборок, которую можно ожидать для выборок заданного размера:

def variance(xs):
    '''Несмещенная дисперсия (варианс) при n <= 30'''
    x_hat = xs.mean() 
    n = len(xs)
    n = n-1 if n in range(1, 30) else n  
    square_deviation = lambda x : (x – x_hat) ** 2 
    return sum( map(square_deviation, xs) ) / n

def standard_deviation(xs):
    return sp.sqrt(variance(xs))

def standard_error(xs):
    return standard_deviation(xs) / sp.sqrt(len(xs))

Стандартная ошибка среднего таким образом связана с двумя факторами:

• Размером выборки

• Популяционным стандартным отклонением

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

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

Примеры исходного кода для этого поста находятся в моем репо на Github. Все исходные данные взяты в репозитории автора книги.

Темой следующего поста, поста №2, будет разница между выборками и популяцией, а также интервал уверенности. Да-да, именно интервал уверенности, а не доверительный интервал.

Python is a great language for doing data analysis, primarily because of the fantastic ecosystem of data-centric python packages. Pandas is one of those packages and makes importing and analyzing data much easier.

Pandas dataframe.sem() function return unbiased standard error of the mean over requested axis. The standard error (SE) of a statistic (usually an estimate of a parameter) is the standard deviation of its sampling distribution[1] or an estimate of that standard deviation. If the parameter or the statistic is the mean, it is called the standard error of the mean (SEM).

Syntax : DataFrame.sem(axis=None, skipna=None, level=None, ddof=1, numeric_only=None, **kwargs)

Parameters :
axis : {index (0), columns (1)}
skipna : Exclude NA/null values. If an entire row/column is NA, the result will be NA
level : If the axis is a MultiIndex (hierarchical), count along a particular level, collapsing into a Series
ddof : Delta Degrees of Freedom. The divisor used in calculations is N – ddof, where N represents the number of elements.
numeric_only : Include only float, int, boolean columns. If None, will attempt to use everything, then use only numeric data. Not implemented for Series

Return : sem : Series or DataFrame (if level specified)

For link to the CSV file used in the code, click here

Example #1: Use sem() function to find the standard error of the mean of the given dataframe over the index axis.

Let’s use the dataframe.sem() function to find the standard error of the mean over the index axis.

Output :

When we include the NaN values then it will cause that particular row or column to be NaN

Python is a great language for doing data analysis, primarily because of the fantastic ecosystem of data-centric python packages. Pandas is one of those packages and makes importing and analyzing data much easier.

Pandas dataframe.sem() function return unbiased standard error of the mean over requested axis. The standard error (SE) of a statistic (usually an estimate of a parameter) is the standard deviation of its sampling distribution[1] or an estimate of that standard deviation. If the parameter or the statistic is the mean, it is called the standard error of the mean (SEM).

Syntax : DataFrame.sem(axis=None, skipna=None, level=None, ddof=1, numeric_only=None, **kwargs)

Parameters :
axis : {index (0), columns (1)}
skipna : Exclude NA/null values. If an entire row/column is NA, the result will be NA
level : If the axis is a MultiIndex (hierarchical), count along a particular level, collapsing into a Series
ddof : Delta Degrees of Freedom. The divisor used in calculations is N – ddof, where N represents the number of elements.
numeric_only : Include only float, int, boolean columns. If None, will attempt to use everything, then use only numeric data. Not implemented for Series

Return : sem : Series or DataFrame (if level specified)

For link to the CSV file used in the code, click here

Example #1: Use sem() function to find the standard error of the mean of the given dataframe over the index axis.

Let’s use the dataframe.sem() function to find the standard error of the mean over the index axis.

Output :

When we include the NaN values then it will cause that particular row or column to be NaN

ch02_Inference

# -*- coding: utf-8 -*-

# Системные библиотеки

import scipy as sp
import pandas as pd    
from scipy import stats

# Графические настройки 

import matplotlib.pyplot as plt
   
from matplotlib import rcParams
rcParams['font.family']     = 'sans-serif'
rcParams['font.sans-serif'] = ['Ubuntu Condensed']
rcParams['figure.figsize']  = (4, 3.3)
rcParams['legend.fontsize'] = 10
rcParams['xtick.labelsize'] = 9
rcParams['ytick.labelsize'] = 9

def saveplot(dest):
    plt.tight_layout()
    plt.savefig('images/' + dest)  

def load_data(fname):
    '''Загрузка данных из файла fname'''
    return pd.read_csv('data/ch02/' + fname, 't')  

def ex_2_1():
    '''Загрузка данных времени задержки на веб-сайте'''
    return load_data('dwell-times.tsv')[0:5:]  # как вариант, head()
    
ex_2_1()

def ex_2_2():
    '''Визуализация времени задержки'''
    load_data('dwell-times.tsv')['dwell-time'].hist(bins=50)
    plt.xlabel('Время задержки, сек.')
    plt.ylabel('Частота')
    #saveplot('ex_2_2.png')  
    plt.show()    
    
ex_2_2()

def ex_2_3():
    '''Визуализация времени задержки
       на полулогарифмическом графике'''
    load_data('dwell-times.tsv')['dwell-time'].plot.hist(bins=20, logy=True, grid=True)
    plt.xlabel('Время задержки, сек.')
    plt.ylabel('Логарифмическая частота')
    #saveplot('ex_2_3.png')  
    plt.show()    
    
ex_2_3()

def ex_2_4():
    '''Показать сводные статистики набора данных,
       подчиняющегося экспоненциональному распределению'''
    ts = load_data('dwell-times.tsv')['dwell-time']
    print('Среднее:', ts.mean())    
    print('Медиана:', ts.median())
    print('СКО:    ', ts.std())
    
ex_2_4()

Среднее: 93.2014074074
Медиана: 64.0
СКО:     93.9697240252

# Определение служебных функций по обработке дат

def with_parsed_date(df):
    '''Привести поле date к типу date-time'''
    df['date'] = pd.to_datetime(df['date'], errors='ignore')
    return df

def filter_weekdays(df): 
    '''Отфильтровать по будним дням'''
    return df[df['date'].index.dayofweek < 5]  # понедельник..пятница

def mean_dwell_times_by_date(df):
    '''Среднесуточные времена задержки'''
    df.index = with_parsed_date(df)['date']
    return df.resample('D').mean()  # перегруппировать  

def daily_mean_dwell_times(df):
    '''Средние времена задержки с фильтрацией - только по будним дням'''
    df.index = with_parsed_date(df)['date']
    df = filter_weekdays(df)
    return df.resample('D').mean() 

def ex_2_5():
    '''Распределение значений в будние дни'''
    df  = load_data('dwell-times.tsv')    
    means = daily_mean_dwell_times(df)
    print('Среднее: %f' % (means.mean()))    
    print('Медиана: %f' % (means.median()))
    print('СКО:     %f' % (means.std()))
    
ex_2_5()

Среднее: 90.210429
Медиана: 90.136612
СКО:     3.722343

def ex_2_6():
    '''Построить гистограмму значений 
       времени задержки в будние дни'''
    df = load_data('dwell-times.tsv')
    daily_mean_dwell_times(df)['dwell-time'].hist(bins=20)
    plt.xlabel('Время задержки по будним дням, сек.')
    plt.ylabel('Частота')
    #saveplot('ex_2_6.png')  
    plt.show()    
    
ex_2_6()

def ex_2_7():
    '''Подгонка нормальной кривой поверх гистограммы'''
    df = load_data('dwell-times.tsv')
    means = daily_mean_dwell_times(df)['dwell-time'].dropna() 
    ax = means.hist(bins=20, normed=True)
    xs = sorted(means)    # корзины
    df = pd.DataFrame()
    df[0] = xs
    df[1] = stats.norm.pdf(xs, means.mean(), means.std())
    df.plot(0, 1, linewidth=2, color='r', legend=None, ax=ax)
    plt.xlabel('Время задержки по будним дням, сек.')
    plt.ylabel('Плотность')
    plt.grid(True)
    #saveplot('ex_2_7.png')  
    plt.show()
    
ex_2_7()

# Определение функции стандартной ошибки среднего

def variance(xs):
    '''Вычисление дисперсии,
       несмещенная дисперсия при n <= 30'''
    x_hat = xs.mean() 
    n = len(xs)
    n = n-1 if n in range(1, 30) else n  
    square_deviation = lambda x : (x - x_hat) ** 2 
    return sum( map(square_deviation, xs) ) / n

def standard_deviation(xs):
    '''Вычисление стандартного отклонения'''
    return sp.sqrt(variance(xs))

def standard_error(xs):
    '''Вычисление стандартной ошибки'''
    return standard_deviation(xs) / sp.sqrt(len(xs))

def ex_2_8(): 
    '''Вычислить стандартную ошибку 
       средних значений за определенный день'''
    may_1 = '2015-05-01'
    df = with_parsed_date( load_data('dwell-times.tsv') )  
    filtered = df.set_index( ['date'] )[may_1]
    se = standard_error( filtered['dwell-time'] )
    print('Стандартная ошибка (СОШ):', se)

ex_2_8()

Стандартная ошибка (СОШ): 3.62734027309

def confidence_interval(p, xs):
    '''Доверительный интервал'''
    x_hat = xs.mean()
    se = standard_error(xs)
    '''критическое значение z
    Критическое значение z - это число стандартных отклонений, на которые 
    нужно отойти от среднего значения нормального распределения, чтобы захватить 
    долю данных, связанную с нужным доверительным интервалом.'''
    z_crit = stats.norm.ppf(1 - (1-p) / 2)  #q=0.975 -> 1.96
    return [x_hat - z_crit * se, x_hat + z_crit * se]

def ex_2_9():
    '''Вычислить доверительный интервал 
       для данных за определенный день'''
    may_1 = '2015-05-01'
    df    = with_parsed_date( load_data('dwell-times.tsv') )  
    filtered = df.set_index( ['date'] )[may_1]
    ci = confidence_interval(0.95, filtered['dwell-time'])
    print('Доверительный интервал (ДИ): ', ci)
    
ex_2_9()

Доверительный интервал (ДИ):  [83.53415272762004, 97.753065317492741]

def ex_2_10():
    '''Сводные статистики данных, полученных 
       в результате вирусной кампании'''
    ts = load_data('campaign-sample.tsv')['dwell-time']     
    print('n:      ', ts.count())
    print('Среднее:', ts.mean())
    print('Медиана:', ts.median())
    print('СКО:    ', ts.std())
    print('СОШ:    ', standard_error(ts))    

ex_2_10()

n:       300
Среднее: 130.22
Медиана: 84.0
СКО:     136.133707144
СОШ:     7.84657283999

def ex_2_11():
    '''Доверительный интервал данных,
       полученных в результате вирусной кампании'''
    ts = load_data('campaign-sample.tsv')['dwell-time']     
    print('ДИ:', confidence_interval(0.95, ts))

ex_2_11()    

ДИ: [114.84099983154137, 145.59900016845864]

'''Проверка даты''' 
# будние дни: 0..4, выходные дни: 5-6
d = pd.to_datetime('2015 6 6') 
d.weekday() in [5,6]

def ex_2_12(): 
    '''Построить график времени ожидания 
       по всем дням без фильтра'''
    df = load_data('dwell-times.tsv')
    means = mean_dwell_times_by_date(df)['dwell-time']
    means.hist(bins=20)
    plt.xlabel('Ежедневное время ожидания без фильтра, сек.')
    plt.ylabel('Частота')
    #saveplot('ex_2_12.png') 
    plt.show()     

ex_2_12()

def ex_2_13():
    '''Сводные статистики данных,
       отфильтрованных только по выходным дням'''
    df = with_parsed_date( load_data('dwell-times.tsv') )
    df.index = df['date']
    df = df[df['date'].index.dayofweek > 4]   # суббота-воскресенье
    weekend_times = df['dwell-time']
  
    print('n:      ', weekend_times.count())
    print('Среднее:', weekend_times.mean())
    print('Медиана:', weekend_times.median())
    print('СКО:    ', weekend_times.std())
    print('СОШ:    ', standard_error(weekend_times))        
    
ex_2_13()

n:       5860
Среднее: 117.786860068
Медиана: 81.0
СКО:     120.652340772
СОШ:     1.57597703625

def pooled_standard_error(a, b, unbias=False):
    '''Объединенная стандартная ошибка'''
    std1 = a.std(ddof=0) if unbias==False else a.std()  # ddof=0 = смещенное значение
    std2 = b.std(ddof=0) if unbias==False else b.std()
    x = std1 ** 2 / a.count()
    y = std2 ** 2 / b.count()
    return sp.sqrt(x + y)    

def z_stat(a, b, unbias=False):
    '''z-статистика'''
    return (a.mean() - b.mean()) / pooled_standard_error(a, b, unbias)    

def z_test(a, b): 
    '''z-тест'''
    # ИФР нормального распределения
    return stats.norm.cdf([ z_stat(a, b) ])

def ex_2_14():
    '''Сравнение работоспособности двух вариантов
       дизайна веб-сайта на основе z-теста'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    a = groups.get_group(0)
    b = groups.get_group(1) 
    
    print('a n: ', a.count())
    print('b n: ', b.count())
    print('z-статистика:', z_stat(a, b))
    print('p-значение:', z_test(a, b))
    
ex_2_14()

a n:  284
b n:  16
z-статистика: -1.64674381801
p-значение: [ 0.04980536]

def pooled_standard_error_t(a, b):
    '''Объединенная стандартная ошибка для t-теста'''
    return sp.sqrt(standard_error(a) ** 2 + 
                   standard_error(b) ** 2) 

t_stat = z_stat

def ex_2_15():
    '''Вычисление t-статистики 
       двух вариантов дизайна веб-сайта'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    a = groups.get_group(0)
    b = groups.get_group(1)    
    return t_stat(a, b)
    
ex_2_15()

def t_test(a, b):
    '''Выполнение проверки на основе t-теста'''
    df = len(a) + len(b) - 2
    # http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.stats.t.html
    return stats.t.sf([ abs(t_stat(a, b)) ], df)      # функция выживания (1-cdf иногда точнее)

def ex_2_16():
    '''Сравнение работоспособности двух вариантов 
       дизайна веб-сайта на основе t-теста'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    a = groups.get_group(0)
    b = groups.get_group(1)   
    return t_test(a, b)

ex_2_16()

def t_test_verbose(a, sample2=None, mean=None, fn=None):
    '''Служебная функция с подробной информацией 
       результата t-теста'''
    abar = a.mean()
    avar = a.var(ddof=1)  # несмещенное значение
    na = a.size
    adof = na - 1
    conf_int = stats.t.interval(0.95, len(a)-1, 
                                loc=sp.mean(a), scale=stats.sem(a))

    if type(a) == type(sample2):
        bbar = sample2.mean()
        bvar = sample2.var(ddof=1)
        nb = sample2.size
        bdof = nb - 1

        dof = (avar/na + bvar/nb)**2 /   
              (avar**2/(na**2*adof) + bvar**2/(nb**2*bdof))
        return {'p-значение'            : 
                   fn(a, sample2, equal_var=False).pvalue,  #  выполняет t-тест Уэлша   
                'степени свободы df'    : dof,  #t_test( a, b ),   
                'доверительный интервал': conf_int,         # conf-int [-78.9894629402365 10.890871390940724]
                'n1'           : a.count(),
                'n2'           : sample2.count(),
                'среднее x'    : a.mean(),
                'среднее y'    : sample2.mean(),
                'дисперсия x'  : a.var(),
                'дисперсия y'  : sample2.var(),
                't-статистика' : fn( a, sample2, equal_var=False ).statistic} 
    else:
        dof = (avar/na) / (avar/(na*adof))
        return {'p-значение'            : fn(a, mean).pvalue,    
                'степени свободы df'    : dof,    
                'доверительный интервал': conf_int, 
                'n1'                    : a.count(),
                'среднее x'             : a.mean(),
                'дисперсия x'           : a.var(),
                't-статистика'          : fn(a, mean).statistic} 

def ex_2_17():
    '''Двухсторонний t-тест'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    a = groups.get_group(0)
    b = groups.get_group(1)    
    return t_test_verbose(a, sample2=b, fn=stats.ttest_ind)  # выполняет t-тест Уэлша

ex_2_17()

{'n1': 284,
 'n2': 16,
 'p-значение': 0.12756432502462475,
 't-статистика': -1.5985205593851322,
 'дисперсия x': 10463.941024237296,
 'дисперсия y': 6669.8666666666668,
 'доверительный интервал': (76.002631987995969, 99.898776462708256),
 'среднее x': 87.950704225352112,
 'среднее y': 122.0,
 'степени свободы df': 17.761382349686098}

def ex_2_18():
    '''Одновыборочный t-тест'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    b = groups.get_group(1) 
    return t_test_verbose(b, mean=90, fn=stats.ttest_1samp)  

ex_2_18()

{'n1': 16,
 'p-значение': 0.13789520958229415,
 't-статистика': 1.5672973291495713,
 'дисперсия x': 6669.8666666666668,
 'доверительный интервал': (78.481527665903897, 165.51847233409609),
 'среднее x': 122.0,
 'степени свободы df': 15.0}

def ex_2_19():
    '''Построение графика синтетических времен задержки 
       путем извлечения бутстрап-выборок'''
    groups = load_data('new-site.tsv').groupby('site')['dwell-time']
    b = groups.get_group(1) 
    xs = [b.sample(len(b), replace=True).mean() for _ in range(1000)] 
    pd.Series(xs).hist(bins=20)
    plt.xlabel('Бутстрапированные средние времени задержки, сек.')
    plt.ylabel('Частота') 
    #saveplot('ex_2_19.png') 
    plt.show() 
    
ex_2_19()

def ex_2_20():
    '''Выборочные средние значения 
       20 разных вариантов дизайна веб-сайта'''
    df = load_data('multiple-sites.tsv')
    return df.groupby('site').aggregate(sp.mean)

ex_2_20()

import itertools

def ex_2_21():
    '''Проверка вариантов дизайна веб-сайта на основе t-тест 
       по принципу "каждый с каждым"'''
    groups = load_data('multiple-sites.tsv').groupby('site')
    alpha = 0.05
    
    pairs = [list(x)   # найти сочетания из n по k
             for x in itertools.combinations(range(len(groups)), 2)]  

    for pair in pairs:
        gr, gr2 = groups.get_group( pair[0] ), groups.get_group( pair[1] )
        site_a, site_b = pair[0], pair[1]
        a, b = gr['dwell-time'], gr2['dwell-time']  
        p_val = stats.ttest_ind(a, b, equal_var = False).pvalue  
        if p_val < alpha: 
            print('Варианты веб-сайта %i и %i значимо различаются: %f' 
                  % (site_a, site_b, p_val))
    
ex_2_21()

Варианты веб-сайта 0 и 6 значимо различаются: 0.005534
Варианты веб-сайта 0 и 7 значимо различаются: 0.035068
Варианты веб-сайта 0 и 10 значимо различаются: 0.006881
Варианты веб-сайта 0 и 15 значимо различаются: 0.046591
Варианты веб-сайта 2 и 6 значимо различаются: 0.018257
Варианты веб-сайта 2 и 10 значимо различаются: 0.029948
Варианты веб-сайта 3 и 6 значимо различаются: 0.038758
Варианты веб-сайта 4 и 6 значимо различаются: 0.027908
Варианты веб-сайта 4 и 10 значимо различаются: 0.046962
Варианты веб-сайта 6 и 8 значимо различаются: 0.039915
Варианты веб-сайта 6 и 9 значимо различаются: 0.026231
Варианты веб-сайта 6 и 11 значимо различаются: 0.035796
Варианты веб-сайта 6 и 12 значимо различаются: 0.004079
Варианты веб-сайта 6 и 13 значимо различаются: 0.024710
Варианты веб-сайта 6 и 14 значимо различаются: 0.002538
Варианты веб-сайта 6 и 16 значимо различаются: 0.033198
Варианты веб-сайта 6 и 17 значимо различаются: 0.004953
Варианты веб-сайта 6 и 18 значимо различаются: 0.030166
Варианты веб-сайта 6 и 19 значимо различаются: 0.019361
Варианты веб-сайта 7 и 12 значимо различаются: 0.027872
Варианты веб-сайта 7 и 14 значимо различаются: 0.017185
Варианты веб-сайта 7 и 17 значимо различаются: 0.031118
Варианты веб-сайта 9 и 10 значимо различаются: 0.046504
Варианты веб-сайта 10 и 12 значимо различаются: 0.003978
Варианты веб-сайта 10 и 13 значимо различаются: 0.042570
Варианты веб-сайта 10 и 14 значимо различаются: 0.002518
Варианты веб-сайта 10 и 17 значимо различаются: 0.006264
Варианты веб-сайта 10 и 19 значимо различаются: 0.031693
Варианты веб-сайта 12 и 15 значимо различаются: 0.037149
Варианты веб-сайта 14 и 15 значимо различаются: 0.022799
Варианты веб-сайта 15 и 17 значимо различаются: 0.041280

def ex_2_22():
    '''Проверка вариантов дизайна веб-сайта на основе t-теста 
       против исходного (0)'''
    groups = load_data('multiple-sites.tsv').groupby('site')
    alpha = 0.05 
    baseline = groups.get_group(0)['dwell-time']
    for site_a in range(1, len(groups)):
        a = groups.get_group( site_a )['dwell-time']
        p_val = stats.ttest_ind(a, baseline, equal_var = False).pvalue 
        if p_val < alpha: 
            print('Вариант %i веб-сайта значимо отличается от исходного: %f' 
                  % (site_a, p_val))
    
ex_2_22()

Вариант 6 веб-сайта значимо отличается от исходного: 0.005534
Вариант 7 веб-сайта значимо отличается от исходного: 0.035068
Вариант 10 веб-сайта значимо отличается от исходного: 0.006881
Вариант 15 веб-сайта значимо отличается от исходного: 0.046591

def ex_2_23():
    '''Проверка вариантов дизайна веб-сайта на основе t-теста 
       против исходного (0) с поправкой Бонферрони'''
    groups = load_data('multiple-sites.tsv').groupby('site')
    alpha = 0.05 / len(groups)
    baseline = groups.get_group(0)['dwell-time']
    for site_a in range(1, len(groups)):
        a = groups.get_group(site_a)['dwell-time']
        p_val = stats.ttest_ind(a, baseline, equal_var = False).pvalue 
        if p_val < alpha: 
            print('Вариант %i веб-сайта отличается от исходного: %f' 
                  % (site_a, p_val))
    
ex_2_23()

def ex_2_Fisher():
    '''Визуализация разных F-распределений на графике'''
    mu = 0
    d1_values, d2_values = [4, 9, 49], [95, 90, 50]
    linestyles = ['-', '--', ':', '-.']
    x = sp.linspace(0, 5, 101)[1:] 
    ax = None
    for (d1, d2, ls) in zip(d1_values, d2_values, linestyles):
        dist = stats.f(d1, d2, mu)
        df  = pd.DataFrame( {0:x, 1:dist.pdf(x)} )   
        ax = df.plot(0, 1, ls=ls, 
                     label=r'$d_1=%i, d_2=%i$' % (d1,d2), ax=ax)

    plt.xlabel('$x$nF-статистика')
    plt.ylabel('Плотность вероятности n$p(x|d_1, d_2)$')
    plt.grid(True)
    #saveplot('ex_2_Fisher.png') 
    plt.show()

ex_2_Fisher()

def ssdev(xs):
    '''Сумма квадратов отклонений между 
       каждым элементом и средним по выборке'''
    x_hat = xs.mean() 
    square_deviation = lambda x : (x - x_hat) ** 2 
    return sum( map(square_deviation, xs) ) 

def f_test(groups):
    '''Проверка на основе F-теста'''
    m, n = len(groups), sum(groups.count())
    df1, df2 = m - 1, n - m 
    ssw = sum( groups.apply(lambda g: ssdev(g)) )  # внутригрупповая сумма квадратов отклонений
    sst = ssdev( df['dwell-time'] )                # полная сумма квадратов по всему набору
    ssb = sst - ssw                                # межгрупповая сумма квадратов отклонений
    msb = ssb / df1                                # усредненная межгрупповая
    msw = ssw / df2                                # усредненная внутригрупповая
    f_stat = msb / msw
    return stats.f.sf(f_stat, df1, df2)            # функция выживания (иногда точнее, чем 1-cdf )

def ex_2_24():
    '''Проверка вариантов дизайна веб-сайта на основе F-теста'''
    df = load_data('multiple-sites.tsv')
    groups = df.groupby('site')['dwell-time']
    return f_test(groups)
   
ex_2_24()

def ex_2_25():
    '''Визуализация распределений всех вариантов 
       дизайна веб-сайта на одной коробчатой диаграмме'''
    df = load_data('multiple-sites.tsv')
    df.boxplot(by='site', showmeans=True)
    plt.xlabel('Номер дизайна веб-сайта')
    plt.ylabel('Время задержки, сек.')
    plt.title('')
    plt.suptitle('')
    plt.gca().xaxis.grid(False)
    #saveplot('ex_2_25.png') 
    plt.show()

ex_2_25()

def ex_2_26():
    '''Проверка вариантов 0 и 10 дизайна веб-сайта на основе F-теста'''
    df = load_data('multiple-sites.tsv')
    groups   = df.groupby('site')['dwell-time']
    site_0   = groups.get_group(0) 
    site_10  = groups.get_group(10)
    _, p_val = stats.ttest_ind(site_0, site_10, equal_var=False)
    return p_val

ex_2_26()

def ex_2_27():
    '''Проверка вариантов 0 и 6 дизайна веб-сайта на основе F-теста'''
    df = load_data('multiple-sites.tsv')
    groups   = df.groupby('site')['dwell-time']
    site_0   = groups.get_group(0) 
    site_6   = groups.get_group(6)
    _, p_val = stats.ttest_ind(site_0, site_6, equal_var=False)
    return p_val
    
ex_2_27()

def pooled_standard_deviation(a, b):
    '''Объединенное стандартное отклонение 
       (не объединенная стандартная ошибка)'''
    return sp.sqrt( standard_deviation(a) ** 2 +
                    standard_deviation(b) ** 2)

def ex_2_28():
    '''Вычисление интервального индекса d Коэна 
       для варианта дизайна веб-сайта под номером 6'''
    df = load_data('multiple-sites.tsv')
    groups = df.groupby('site')['dwell-time']
    a      = groups.get_group(0)
    b      = groups.get_group(6)
    return (b.mean() - a.mean()) / pooled_standard_deviation(a, b)

ex_2_28()

  Перевод

  Ссылка на автора

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

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

Одним из наиболее важных применений стандартного отклонения является сравнение двух наборов данных. Если два набора данных имеют одинаковое среднее значение, это не значит, что они точно одинаковы, верно? Например, наборы данных199, 200, 201а также0, 200, 400оба имеют одинаковое среднее значение (200), однако первый набор данных имеет очень небольшое стандартное отклонение (s= 1) по сравнению со вторым набором данных (s= 200).

(13–53.5)=(-40.5) and square of (-40.5) will be: 1,640.25

import numpy as np
dataset=[13, 22, 26, 38, 36, 42,49, 50, 77, 81, 98, 110]print('Mean:', np.mean(dataset))
print('Standard Deviation:', np.std(dataset))Mean:53.5
Standard Deviation: 29.694275542602483

City A= 36, 37, 36, 34, 39, 33, 30, 30, 32, 31, 31, 32, 32, 33, 35City B= 41, 35, 28, 29, 25, 36, 36, 32, 38, 40, 40, 34, 31, 28, 30

City_A=[36,37,36,34,39,33,30,30,32,31,31,32,32,33,35] 
City_B=[41,35,28,29,25,36,36,32,38,40,40,34,31,28,30] Mean_City_A=np.mean(City_A)
Mean_City_B=np.mean(City_B) STDV_City_A=np.std(City_A)
STDV_City_B=np.std(City_B)

# Create a figure with customized size
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111) # Set the axis lables
ax.set_xlabel('Day', fontsize = 18)
ax.set_ylabel('Temperature', fontsize = 18)

# X axis is day numbers from 1 to 15
xaxis = np.array(range(1,16))

# Line color for error bar
color_City_A = 'red'
color_City_B = 'darkgreen' # Line style for each dataset
lineStyle_City_A={"linestyle":"--", "linewidth":2, "markeredgewidth":2, "elinewidth":2, "capsize":3}lineStyle_City_B={"linestyle":"-", "linewidth":2, "markeredgewidth":2, "elinewidth":2, "capsize":3}

# Create an error bar for each datasetline_City_A=ax.errorbar(xaxis, City_A, yerr=STDV_City_A, **lineStyle_City_A, color=color_City_A, label='City A')line_City_B=ax.errorbar(xaxis, City_B, yerr=STDV_City_B, **lineStyle_City_B, color=color_City_B, label='City B')

# Label each dataset on the graph, xytext is the label's position for i, txt in enumerate(City_A):   
     ax.annotate(txt, xy=(xaxis[i], City_A[i]), 
             xytext=(xaxis[i]+0.03,City_A[i]+0.3),
             color=color_City_A) for i, txt in enumerate(City_B):        
     ax.annotate(txt, xy=(xaxis[i], City_B[i]), 
            xytext=(xaxis[i]+0.03,
            City_B[i]+0.3),color=color_City_B)

# Customize the legend font and handle 
lengthparams = {'legend.fontsize': 13,          'legend.handlelength': 2}
plt.rcParams.update(params) # Customize the font
font = {'family' : 'Arial', 'weight':'bold','size': 12}
matplotlib.rc('font', **font) # Draw a grid for the graph and set face color to the graph
ax.grid(color='lightgrey', linestyle='-')
ax.set_facecolor('w')

Полное руководство по проверке гипотез для специалистов по данным, использующих Python

Ясно объяснен с примерами вопросов для исследования, шагами решения и полными кодами

Проверка гипотез — важная часть статистики и анализа данных. В большинстве случаев практически невозможно получить данные от всего населения. В этом случае мы берем образец и делаем оценки или утверждения относительно всего населения. Эти предположения или утверждения являются гипотезами. Проверка гипотез — это процесс проверки наличия доказательств для отклонения этой гипотезы.

Проверка гипотез обычно проводится по соотношению и среднему значению.

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

Я объясню процесс проверки гипотез шаг за шагом для всех четырех категорий индивидуально с примерами.

Для этого упражнения я использовал среду Jupyter Notebook. Если у вас его нет, вы можете использовать любой ноутбук или IDE по вашему выбору.

Блокнот Google Collab тоже подойдет. Google Collab — это умный блокнот. В нем предустановлены эти общие библиотеки.

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

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

Для этой демонстрации я использовал набор данных Heart из Kaggle. Не стесняйтесь загрузить набор данных для своей практики. Здесь я импортирую пакеты и набор данных:

import pandas as pd
import numpy as np
import statsmodels.api as sm
import scipy.stats.distributions as distdf = pd.read_csv('Heart.csv')
df.head()

Последний столбец набора данных — «AHD». Это если у человека болезнь сердца. Вопрос исследования для этого раздела:

«Доля населения Ирландии, страдающего сердечными заболеваниями, составляет 42%. В США больше людей страдают сердечными заболеваниями »?

Теперь найдите ответ на этот исследовательский вопрос шаг за шагом.

Шаг 1: определите нулевую гипотезу и альтернативную гипотезу.

В этой задаче нулевая гипотеза состоит в том, что доля населения, страдающего сердечными заболеваниями, в США меньше или равна 42%. Но если мы проверим, что меньше, то будет автоматически. Итак, я делаю это только равным.

Альтернативная гипотеза состоит в том, что доля населения США, страдающего сердечными заболеваниями, составляет более 42%.

Ho: p0 = 0.42  #null hypothesis
Ha: p > 0.42  #alternative hypothesis

Посмотрим, сможем ли мы найти доказательства, чтобы отвергнуть нулевую гипотезу.

Шаг 2: Предположим, что приведенный выше набор данных является репрезентативной выборкой из населения США. Итак, рассчитайте долю населения США, страдающую сердечными заболеваниями.

p_us = len(df[df['AHD']=='Yes'])/len(df)

Доля населения выборки, страдающая сердечными заболеваниями, составляет 0,46 или 46%. Этот процент превышает нулевую гипотезу. Это 42%.

Но вопрос в том, значительно ли больше 42%. Если мы возьмем другую простую случайную выборку, наблюдаемая в настоящее время доля населения (46%) может быть другой.

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

Шаг 3. Рассчитайте статистику теста:

Вот формула тестовой статистики:

Мы используем эту формулу для стандартной ошибки:

В этой формуле p0 составляет 0,42 (согласно нулевой гипотезе), а n — размер выборки. Теперь вычислите Стандартную ошибку и статистику теста:

se = np.sqrt(0.42 * (1-0.42) / len(df))

Найдите статистику теста, используя формулу для статистики теста выше:

#Best estimate
be = p_us  #hypothesized estimate
he = 0.42test_stat = (be - he)/se

Статистика теста составила 1,3665.

Шаг 4: Рассчитайте p-значение

Эта статистика теста также называется z-оценкой. Вы можете найти p-значение из z_table или вы можете найти p-значение из этой формулы в python.

pvalue = 2*dist.norm.cdf(-np.abs(test_stat))

Значение p равно 0,1718. Это означает, что доля выборки (46% или 0,46) составляет 0,1718 нулевых стандартных ошибок выше нулевой гипотезы.

Шаг 5. Сделайте вывод по p-значению

Считайте, что уровень значимости альфа составляет 5% или 0,05. Уровень значимости 5% или меньше означает, что существует вероятность 95% или больше, что результаты не являются случайными.

Здесь p-значение больше, чем наш рассмотренный уровень значимости 0,05. Итак, мы не можем отвергнуть нулевую гипотезу. Это означает, что нет значительной разницы в доле населения, страдающего сердечными заболеваниями, в Ирландии и США.

Проверка гипотез о разнице в двух пропорциях

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

Здесь мы собираемся проверить, отличается ли доля женщин с сердечными заболеваниями в популяции мужчин с сердечными заболеваниями.

Шаг 1. Установите нулевую гипотезу, альтернативную гипотезу и уровень значимости.

Здесь мы хотим проверить, есть ли разница между долей мужчин и женщин, страдающих сердечными заболеваниями. Начнем с предположения, что разницы нет.

Ho: p1 -p2 = 0

Это наша нулевая гипотеза. Здесь p1 — доля женщин с сердечными заболеваниями в популяции, а p2 — доля мужчин, страдающих сердечными заболеваниями.

Какая может быть альтернативная гипотеза?

Альтернативная гипотеза может быть, есть разница.

Ha: p1 - p2 != 0

Давайте воспользуемся уровнем значимости 0,1 или 10%.

Шаг 2. Подготовьте диаграмму, показывающую долю мужчин и женщин с сердечными заболеваниями в общей численности мужского и женского населения.

df['Gender'] = df.Sex.replace({1: "Male", 0: "Female"})
p = df.groupby("Gender")['AHD'].agg([lambda z: np.mean(z=='Yes'), "size"])
p.columns = ["HeartDisease", 'Total']
p

Шаг 3. Рассчитайте статистику теста.

Мы будем использовать ту же формулу для статистики теста, что и раньше. Наилучшая оценка — p1 — p2. Здесь p1 — доля женщин с сердечными заболеваниями в популяции, а p2 — доля мужчин с сердечными заболеваниями.

#Best estimate is p1 - p2. Get p1 and p2 from the chart p above
p_fe = p.HeartDisease.Female
p_male = p.HeartDisease.Male

Стандартная ошибка для двух пропорций населения рассчитывается по следующей формуле:

Здесь p — общая доля населения в выборке с сердечными заболеваниями. n1 и n2 — общее количество женского и мужского населения в выборке.

p = p_us  #calculated in the beginning of the previous example
n1 = p.Total.Female
n2 = p.Total.Male
se = np.sqrt(p_us*(1-p_us)*(1/n1 + 1/n2))

Теперь используйте эту стандартную ошибку и вычислите статистику теста.

#calculate the best estimate
be = p_fe - p_male  #Calculate the hypothesized estimate
#Our null hypothesis is p1 - p2 = 0he = 0  #Calculate the test statistic
test_statistic = (be - he)/se

Вычисленное значение test_statistic равно -0,296. Это означает, что наблюдаемая разница в пропорциях выборки на 0,296 расчетной стандартной ошибки ниже предполагаемого значения.

Шаг 4. Рассчитайте p-значение

pvalue = 2*dist.norm.cdf(-np.abs(test_statistic)

Значение p равно 0,7675. Это означает, что более 76% случаев мы увидим, что наблюдаемые нами результаты верны, учитывая, что нулевая гипотеза верна.

С другой стороны, p-значение больше, чем уровень значимости (0,1). Итак, у нас недостаточно доказательств, чтобы отвергнуть нулевую гипотезу.

Доля мужчин с сердечными заболеваниями существенно не отличается от доли женщин с сердечными заболеваниями.

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

Это простой процесс проверки гипотез. Мы можем выполнить этот тест, если у нас есть конкретное фиксированное среднее значение для сравнения. Давайте поработаем на примере, чтобы понять процесс.

Это вопрос исследования:

«Проверьте, не превышает ли среднее значение RestBP 135». Здесь RestBP показывает артериальное давление в покое. У нас есть столбец RestBP в DataFrame. Решим эту задачу пошагово.

Шаг 1. Сформулируйте гипотезу.

Нам нужно выяснить, больше ли среднее значение RestBP 135. Предположим, что среднее значение RestBP меньше или равно 135.

Итак, нулевая гипотеза может заключаться в том, что среднее значение RestBP равно 135. Потому что, если мы можем доказать, что среднее значение RestBP больше 135, оно автоматически больше 134 или 130.

Если мы найдем достаточно доказательств, чтобы отклонить нулевую гипотезу, мы можем принять, что среднее значение RestBP больше 135. Это альтернативная гипотеза для данного примера.

Ho: mu = 135
Ha: mu > 135

Мы проверим, можем ли мы отклонить нулевую гипотезу, используя уровень значимости 0,05.

Шаг 2. Проверьте предположения.

Есть два предположения:

1. Выборка должна быть простой случайной выборкой.
2. Данные должны быть нормально распределены.

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

sns.distplot(df.RestBP)

Распределение не совсем нормальное. Но это близко к норме.

Хорошая новость в том, что нам не нужно беспокоиться о нормальности данных. Потому что у нас достаточно большой размер выборки (более 25 данных).

Шаг 3. Рассчитайте статистику теста.

Вот формула для расчета статистики теста:

Сначала вычислите стандартную ошибку, используя формулу ниже:

Здесь S — стандартное отклонение выборки, а n — количество выборок.

std= df.RestBP.std()
n = len(df)
se = std/np.sqrt(n)

Теперь используйте эту стандартную ошибку, чтобы найти статистику теста:

#Best estimate
be = df.RestBP.mean() #Hypothesized estimatehe = 135
test_statistic = (be - he)/se

Статистика теста составила -3,27. Посмотрите на формулу тестовой статистики. Сверху он измеряет расстояние между исходным средним и предполагаемым средним. А внизу — стандартная ошибка.

Итак, этот test_statistic означает, что выборочное среднее на 3,27 стандартной ошибки ниже предполагаемого среднего.

Шаг 4. Сделайте вывод на основании статистики теста.

Преобразуйте этот test_statistic в значение вероятности, чтобы увидеть, является ли эта разница необычной или нет. Мы можем получить значение, используя эту формулу Python:

pvalue = 2*dist.norm.cdf(-np.abs(test_statistic))

Значение p равно 0,001, что меньше уровня значимости (0,05).

Итак, мы можем отвергнуть нулевую гипотезу.

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

Итак, мы отвергаем нулевую гипотезу и принимаем альтернативную гипотезу, основанную на этих выборочных данных.

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

В этом примере мы будем использовать те же данные, столбец RestBP. Но на этот раз, чтобы проверить, есть ли разница между средним RestBP самок и средним RestBP самцов.

Шаг 1. Сформулируйте гипотезу

В качестве нулевой гипотезы начните с утверждения, что среднее значение RestBP у женщин и среднее значение RestBP у мужчин одинаковы. Таким образом, разница между этими двумя средними будет равна нулю.

Альтернативная гипотеза состоит в том, что эти два средства не совпадают. Давайте проведем тест с 10% уровнем значимости.

Ho: mu_female - mu_male = 0
Ha: mu_female - mu_male != 0

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

Шаг 2. Рассчитайте статистику теста.

Формула для статистики теста такая же, как и раньше. Но формула для стандартной ошибки другая.

Здесь s1 и s2 — стандартное отклонение выборки для женского и мужского населения соответственно. n1 и n2 — размер выборки женского и мужского населения. Теперь вычислите стандартную ошибку:

pop_fe = df[df.Gender=='Female'].dropna()
pop_male = df[df.Gender=='Male'].dropna()std_fe = pop_fe.RestBP.std()
std_male = pop_male.RestBP.std()se = np.sqrt(std_fe**2/len(pop_fe) + std_male**2/len(pop_male))

Используйте стандартную ошибку, чтобы получить статистику теста.

#calculate the best estimate
mu_fe = pop_fe.RestBP.mean()     #Mean RestBP for females
mu_male = pop_male.RestBP.mean()  #Mean RestBP for malesmu_diff = mu_fe - mu_male    #hypothesized estimate 
mu_diff_hyp = 0   #null hypothesis: difference of two mean = zerotest_statistic = (be-he)/se

Test_statistic — 1.086. Для информации: наблюдаемая разница в средних значениях «mu_diff» составляет 2,52.

Поскольку мы проверяем, отличаются ли средние значения друг от друга, это двусторонний тест.

Значение p — это вероятность того, что статистика теста меньше 1,086 или больше 1,086.

Шаг 3. Сделайте выводы на основе статистики теста.

Рассчитайте p-значение из этой тестовой статистики в python:

pvalue = 2*dist.norm.cdf(-np.abs(test_statistic))

P-значения оказались равными 0,277. Поскольку это двусторонний тест,

p(z < -1.086) = 0.277

p(z > 1.086) = 0.277

p-значение = 0,277 + 0,277 = 0,554

Это означает, что существует приблизительно 55,4% вероятности того, что наблюдаемый результат или более экстремальный верен, когда верна нулевая гипотеза.

С другой стороны, p-значение намного больше, чем уровень значимости. Итак, мы не можем отвергнуть нулевую гипотезу.

Окончательный вывод, основанный на наблюдаемой разнице между средним RestBP женщин и средним RestBP мужчин, мы не можем поддержать идею о том, что существует значительная разница между двумя средними значениями.

Заключение

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

Рекомендуемая литература