Глава 14. Аналитические методы

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

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

График нормальной вероятности

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

Чтобы продемонстрировать первое обстоятельство, посмотрим на некоторые переменные из датасета о пингвинах. Затем мы исследуем математические свойства нормального распределения. Инструкции по загрузке данных приведены в Jupyter-блокноте этой главы.

Загрузим данные:

penguins = pd.read_csv("penguins_raw.csv")

penguins.shape

(344, 17)

Датасет содержит результаты измерений трех видов пингвинов. В этом примере выберем пингвинов вида Адели:

adelie = penguins.query('Species.str.startswith("Adelie")').copy()

len(adelie)

152

Чтобы проверить, соответствует ли вес пингвинов нормальному распределению, вычислим эмпирическую интегральную функцию распределения (ИФР) результатов измерений:

from empiricaldist import Cdf

weights = adelie["Body Mass (g)"].dropna()

cdf_weights = Cdf.from_seq(weights)

Также вычислим аналитическую ИФР нормального распределения с теми же средним значением и стандартным отклонением:

m, s = weights.mean(), weights.std()

m, s

(3700.662251655629, 458.5661259101348)

from scipy.stats import norm

dist = norm(m, s)

qs = np.linspace(m - 3.5 * s, m + 3.5 * s)

ps = dist.cdf(qs)

Вот как выглядит ИФР фактических данных в сравнении с нормальной моделью:

model_options = dict(color="gray", alpha=0.5, label="модель")

plt.plot(qs, ps, **model_options)

cdf_weights.plot(label="данные")

decorate(ylabel="ИФР")

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

В целом построение ИФР данных и ИФР модели — хороший способ проверить, насколько модель соответствует данным. Однако этот метод зависит от того, хорошо ли мы оценили параметры модели — в данном примере среднее значение и стандартное отклонение.

Альтернативой ему является график нормальной вероятности (normal probability plot), который не зависит от нашего умения оценивать параметры. На графике нормальной вероятности значения y — отсортированные выборочные данные:

ys = np.sort(weights)

А значения x — это соответствующие квантили нормального распределения, вычисленные с использованием метода ppf объекта norm. Этот метод считает квантильную функцию (percent point function, PPF), которая является обратной к ИФР:

n = len(weights)

ps = (np.arange(n) + 0.5) / n

xs = norm.ppf(ps)

Если результаты измерений на самом деле подчиняются нормальному распределению, то координаты y и x должны лежать на прямой линии. Чтобы посмотреть, так ли это, можно с помощью функции linregress подобрать линию регрессии:

from scipy.stats import linregress

results = linregress(xs, ys)

intercept, slope = results.intercept, results.slope

fit_xs = np.linspace(-3, 3)

fit_ys = intercept + slope * fit_xs

Следующий рисунок демонстрирует величины x и y вместе с подогнанной линией:

plt.plot(fit_xs, fit_ys, **model_options)

plt.plot(xs, ys, label="данные")

decorate(xlabel="Стандартное нормальное распределение", ylabel="Масса тела (г)")

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

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

def normal_probability_plot(sample, **options):

"""Построить график нормальной вероятности с подогнанной линией."""

n = len(sample)

ps = (np.arange(n) + 0.5) / n

xs = norm.ppf(ps)

ys = np.sort(sample)

results = linregress(xs, ys)

intercept, slope = results.intercept, results.slope

fit_xs = np.linspace(-3, 3)

fit_ys = intercept + slope * fit_xs

plt.plot(fit_xs, fit_ys, color="gray", alpha=0.5)

plt.plot(xs, ys, **options)

decorate(xlabel="Стандартное нормальное распределение")

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

grouped = adelie.groupby("Sex")

weights_male = grouped.get_group("MALE")["Body Mass (g)"]

normal_probability_plot(weights_male, ls="--", label="Самцы")

weights_female = grouped.get_group("FEMALE")["Body Mass (g)"]

normal_probability_plot(weights_female, label="Самки")

decorate(ylabel="Вес (г)")

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

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

Нормальное распределение

Следующий класс описывает объект, представляющий нормальное распределение. В качестве атрибутов используются параметры mu и sigma2 — среднее значение и дисперсия распределения. Название sigma2 напоминает, что дисперсия — это квадрат стандартного отклонения, которое обычно обозначается как sigma:

class Normal:

"""Представляет нормальное распределение"""

def __init__(self, mu, sigma2):

self.mu = mu

self.sigma2 = sigma2

def __repr__(self):

return f"Normal({self.mu}, {self.sigma2})"

__str__ = __repr__

Для примера создадим объект Normal, представляющий нормальное распределение с тем же средним значением и дисперсией, что и вес самцов пингвинов:

m, s = weights_male.mean(), weights_male.std()

dist_male = Normal(m, s**2)

dist_male

Normal(4043.4931506849316, 120278.25342465754)

И еще один объект Normal со средним значением и дисперсией как у веса самок пингвинов:

m, s = weights_female.mean(), weights_female.std()

dist_female = Normal(m, s**2)

dist_female

Normal(3368.8356164383563, 72565.63926940637)

Далее в класс Normal добавим метод, который генерирует случайную выборку из нормального распределения. Новые методы к существующему классу мы будем добавлять с помощью «магической» команды Jupyter add_method_to, которая определена в модуле thinkstats. Эта команда не является функционалом языка Python — ее можно использовать только в Jupyter-блокнотах:

%%add_method_to Normal

def sample(self, n):

sigma = np.sqrt(self.sigma2)

return np.random.normal(self.mu, sigma, n)

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

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

sample_sum = dist_male.sample(1000) + dist_female.sample(1000)

normal_probability_plot(sample_sum)

decorate(ylabel="Суммарный вес (г)")

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

%%add_method_to Normal

def __add__(self, other):

"""Распределение суммы двух нормальных распределений."""

return Normal(self.mu + other.mu, self.sigma2 + other.sigma2)

У распределения суммы среднее значение представляет собой сумму средних значений, а дисперсия — сумму дисперсий. Теперь, после того как определен специальный метод __add__, мы можем использовать оператор + для «сложения» двух распределений, то есть для вычисления распределения их суммы:

dist_sum = dist_male + dist_female

dist_sum

Normal(7412.328767123288, 192843.8926940639)

Чтобы подтвердить правильность полученного результата, воспользуемся следующим методом, который строит график аналитической ИФР нормального распределения:

%%add_method_to Normal

def plot_cdf(self, n_sigmas=3.5, **options):

mu, sigma = self.mu, np.sqrt(self.sigma2)

low, high = mu - n_sigmas * sigma, mu + n_sigmas * sigma

xs = np.linspace(low, high, 101)

ys = norm.cdf(xs, mu, sigma)

plt.plot(xs, ys, **options)

Вот графики аналитической и эмпирической ИФР суммы случайных выборок:

dist_sum.plot_cdf(**model_options)

Cdf.from_seq(sample_sum).plot(label="выборка")

decorate(xlabel="Суммарный вес (г)", ylabel="ИФР")

Выглядит так, что параметры, которые мы вычислили, верны. Это подтверждает, что можно сложить два нормальных распределения, суммировав их средние значения и дисперсии.

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

n = len(weights_male)

sample_sums_male = [dist_male.sample(n).sum() for i in range(1001)]

Следующий метод создает объект Normal, представляющий распределение сумм. Чтобы вычислить его параметры, умножим среднее значение и дисперсию на n:

%%add_method_to Normal

def sum(self, n):

"""Распределение суммы n значений."""

return Normal(n * self.mu, n * self.sigma2)

Вот распределение суммы n значений веса:

dist_sums_male = dist_male.sum(n)

А вот как оно соотносится с эмпирическим распределением случайной выборки:

dist_sums_male.plot_cdf(**model_options)

Cdf.from_seq(sample_sums_male).plot(label="выборка")

decorate(xlabel="Суммарный вес (г)", ylabel="ИФР")

Аналитическое распределение соответствует распределению выборки, что подтверждает верность метода sum. Таким образом, если составить выборку из n измерений, можно вычислить распределение их суммы.

Распределение выборочных средних

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

%%add_method_to Normal

def __mul__(self, factor):

"""Умножение на константу."""

return Normal(factor * self.mu, factor**2 * self.sigma2)

%%add_method_to Normal

def __truediv__(self, factor):

"""Деление на константу."""

return self * (1/factor)

Чтобы вычислить распределение произведения, мы умножаем среднее значение на значение factor, а дисперсию — на квадрат factor. Это свойство можно использовать для вычисления распределения выборочных средних:

dist_mean_male = dist_sums_male / n

Чтобы убедиться в правильности полученного результата, вычислим также средние значения случайных выборок:

sample_means_male = np.array(sample_sums_male) / n

И сравним нормальную модель с эмпирической ИФР выборочных средних:

dist_mean_male.plot_cdf(**model_options)

Cdf.from_seq(sample_means_male).plot(label="выборка")

decorate(xlabel="Средний вес (г)", ylabel="ИФР")

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

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

standard_error = np.sqrt(dist_mean_male.sigma2)

standard_error

40.591222045992765

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

Таким образом, рассчитаем стандартную ошибку выборочного среднего:

standard_error = weights_male.std() / np.sqrt(n) standard_error

40.59122204599277

Теперь рассмотрим еще один результат, который можно получить при помощи нормальных распределений, — распределение разностей.

Распределение разностей

Объединив все описанные в предыдущем разделе шаги, можно рассчитать распределение выборочных средних для веса самок пингвинов:

n = len(weights_female)

dist_mean_female = dist_female.sum(n) / n

dist_mean_female

Normal(3368.835616438356, 994.0498530055667)

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

%%add_method_to Normal

def __sub__(self, other):

"""Распределение разности."""

return Normal(self.mu - other.mu, self.sigma2 + other.sigma2)

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

• Если изменить знак второго распределения, среднее значение также сменит знак, но дисперсия останется такой же.

• Далее, если добавить первое распределение, дисперсия суммы будет равна сумме дисперсий.

Если у вас еще остались сомнения, давайте их развеем. Вот аналитическое распределение разностей:

dist_diff_means = dist_mean_male - dist_mean_female

dist_diff_means

Normal(674.6575342465753, 2641.697160192656)

А вот случайная выборка разностей:

sample_sums_female = [dist_female.sample(n).sum() for i in range(1001)]

sample_means_female = np.array(sample_sums_female) / n

sample_diff_means = sample_means_male - sample_means_female

Следующий рисунок демонстрирует эмпирическую ИФР случайной выборки и аналитическую ИФР нормального распределения:

dist_diff_means.plot_cdf(**model_options)

Cdf.from_seq(sample_diff_means).plot(label="выборка")

decorate(xlabel="Разность средних весов (г)", ylabel="ИФР")

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

%%add_method_to Normal

def ppf(self, xs):

sigma = np.sqrt(self.sigma2)

return norm.ppf(xs, self.mu, sigma)

5-й и 95-й процентили образуют 90%-ный доверительный интервал:

ci90 = dist_diff_means.ppf([0.05, 0.95])

ci90

array([590.1162635, 759.19880499])

С помощью случайной выборки приходим к похожим результатам:

np.percentile(sample_diff_means, [5, 95])

array([589.01470284, 760.1276391 ])

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

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

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

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

Но многие распределения обладают следующим свойством. Если сгенерировать n значений, а затем сложить их, распределение суммы будет сходиться к нормальному по мере увеличения n. Точнее, если исходное распределение имеет среднее значение m и дисперсию s2, то распределение суммы сходится к нормальному распределению со средним значением n * m и дисперсией n * s2.

Такая закономерность называется центральной предельной теоремой, ЦПТ (Central Limit Theorem, CLT). ЦПТ — один из наиболее полезных инструментов статистического анализа, но с некоторыми оговорками.

• Значения должны быть получены из одного и того же распределения (хотя это требование может быть ослаблено).

• Значения должны быть получены независимо друг от друга. Если они коррелируют, ЦПТ неприменима (хотя и может работать, если корреляция не слишком сильна).

• Значения должны быть получены из распределения с конечными средним значением и дисперсией. Таким образом, ЦПТ неприменима к некоторым распределениям с «длинным хвостом».

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

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

lam = 1

df_sample_expo = pd.DataFrame()

for n in [1, 10, 100]:

df_sample_expo[n] = [np.sum(np.random.exponential(lam, n))

for _ in range(1001)]

Вот средние значения для каждого списка сумм:

df_sample_expo.mean()

1 1.010255

10 9.993695

100 100.416396

dtype: float64

У этого распределения среднее значение равно 1. Поэтому если сложить 10 значений, среднее значение суммы будет близко к 10, а если сложить 100 значений, то среднее значение суммы будет близко к 100.

На следующем рисунке показаны графики нормальной вероятности для трех списков сумм (определение функции normal_plot_samples приведено в Jupyter-блокноте этой главы):

normal_plot_samples(df_sample_expo, ylabel="Сумма экспоненциальных значений")

При n=1 распределение суммы является экспоненциальным, поэтому график нормальной вероятности не выглядит как прямая линия. Но при n=10 распределение суммы становится близким к нормальному, а при n=100 оно уже почти неотличимо от нормального.

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

mu, sigma = 3.0, 1.0

df_sample_lognormal = pd.DataFrame()

for n in [1, 10, 100]:

df_sample_lognormal[n] = [

np.sum(np.random.lognormal(mu, sigma, n)) for _ in range(1001)

]

Вот графики нормальной вероятности для тех же размеров выборки:

normal_plot_samples(df_sample_lognormal, ylabel="Сумма логнормальных значений")

При n=1 нормальная модель плохо описывает распределение, и при n=10 ситуация ненамного лучше. Даже при n=100 хвосты распределения явно отклоняются от модельной прямой.

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

Ограничения центральной предельной теоремы

Распределение Парето еще более скошено, чем логнормальное. При некоторых значениях параметра распределение Парето не имеет конечного среднего значения или дисперсии. В подобных случаях центральная предельная теорема неприменима.

Чтобы это показать, сгенерируем выборку из распределения Парето с параметром alpha=1, имеющего бесконечные среднее значение и дисперсию:

alpha = 1.0

df_sample = pd.DataFrame()

for n in [1, 10, 100]:

df_sample[n] = [np.sum(np.random.pareto(alpha, n)) for _ in range(1001)]

Вот как выглядят графики нормальной вероятности при разных размерах выборки:

normal_plot_samples(df_sample, ylabel="Выборочная сумма распределения Парето")

Даже при n=100 распределение суммы совсем не похоже на нормальное.

Я уже упоминал, что центральная предельная теорема неприменима, если значения коррелируют. Чтобы в этом убедиться, воспользуемся функцией generate_expo_correlated. Она генерирует значения согласно экспоненциальному распределению, при этом последовательная корреляция, то есть корреляция между следующими один за другим элементами в полученной выборке, равна заданному значению rho. Функция определена в Jupyter-блокноте данной главы.

В следующем цикле создается объект DataFrame со столбцом для каждого размера выборки и 1001 суммой в каждом столбце:

rho = 0.8

df_sample = pd.DataFrame()

for n in [1, 10, 100]:

df_sample[n] = [np.sum(generate_expo_correlated(n, rho))

for _ in range(1001)]

Вот графики нормальной вероятности для распределения этих сумм:

normal_plot_samples(df_sample, ylabel="Сумма коррелированных значений")

При rho=0.8 существует сильная корреляция между последовательными элементами, и распределение суммы сходится медленно. Если также присутствует сильная корреляция между отдаленными элементами последовательности, оно может вообще не сойтись.

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

Применение центральной предельной теоремы

Чтобы понять, почему центральная предельная теорема очень полезна, вернемся к примеру из раздела «Проверка разницы средних» главы 9 на с. 166 — проверке наблюдаемой разницы в средней продолжительности беременности для первых и последующих детей. Снова воспользуемся данными NSFG — инструкции по их загрузке приведены в Jupyter-блокноте этой главы.

Воспользуемся функцией get_nsfg_groups, чтобы считать данные и разделить их на соответствующие первым (firsts) и последующим (others) детям:

from nsfg import get_nsfg_groups

live, firsts, others = get_nsfg_groups()

Как вы уже видели, первые дети в среднем рождаются немного позже: наблюдаемая разница составляет около 0.078 недели:

delta = firsts["prglngth"].mean() - others["prglngth"].mean()

delta

0.07803726677754952

Чтобы понять, можно ли появление такого различия объяснить случайностью, примем в качестве нулевой гипотезы, что среднее значение и дисперсия продолжительности беременности на самом деле одинаковы в обеих группах. Тогда их оценку можно получить, используя все случаи рождения живого ребенка (live):

all_lengths = live["prglngth"]

m, s2 = all_lengths.mean(), all_lengths.var()

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

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

def sampling_dist_mean(data, n):

mean, var = data.mean(), data.var()

dist = Normal(mean, var)

return dist.sum(n) / n

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

n1 = firsts["totalwgt_lb"].count()

dist_firsts = sampling_dist_mean(all_lengths, n1)

А вот выборочное распределение для второго и последующих детей:

n2 = others["totalwgt_lb"].count()

dist_others = sampling_dist_mean(all_lengths, n2)

Вычислим выборочное распределение их разности:

dist_diff = dist_firsts - dist_others

dist_diff

Normal(0.0, 0.003235837567930557)

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

Чтобы убедиться, что это распределение аппроксимирует выборочное распределение, оценим его при помощи ресемплинга:

sample_firsts = [np.random.choice(all_lengths, n1).mean() for i in range(1001)]

sample_others = [np.random.choice(all_lengths, n2).mean() for i in range(1001)]

sample_diffs = np.subtract(sample_firsts, sample_others)

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

dist_diff.plot_cdf(**model_options)

Cdf.from_seq(sample_diffs).plot(label="выборка")

plt.axvline(delta, ls=":")

plt.axvline(-delta, ls=":")

decorate(xlabel="Разность продолжительностей беременности", ylabel="ИФР")

В этом примере размеры выборок велики, а асимметрия данных умеренная, поэтому выборочное распределение хорошо аппроксимируется нормальным распределением. Следовательно, можно использовать нормальную ИФР для вычисления p-значения. Следующий метод высчитывает ИФР нормального распределения:

%%add_method_to Normal

def cdf(self, xs):

sigma = np.sqrt(self.sigma2)

return norm.cdf(xs, self.mu, sigma)

Вот вероятность разности величиной с delta, то есть площадь под правым хвостом выборочного распределения, в условиях нулевой гипотезы:

right = 1 - dist_diff.cdf(delta)

right

0.08505405315526993

А вот вероятность отрицательной разности размером с -delta, то есть площадь под левым хвостом:

left = dist_diff.cdf(-delta)

left

0.08505405315526993

Значения переменных left и right одинаковы, поскольку нормальное распределение симметрично. Сумма этих двух значений — вероятность, что разность по модулю равна величине delta:

left + right

0.17010810631053985

Итоговое p-значение, равное 0.170, согласуется с оценкой, вычисленной при помощи ресемплинга в разделе «Проверка разницы средних» главы 9 на с. 166.

Способ, которым мы вычислили это p-значение, аналогичен t-критерию Стьюдента для независимых выборок (independent sample t-test). В библиотеке SciPy определена функция ttest_ind, которая принимает на вход две выборки и вычисляет p-значение для различия в их средних:

from scipy.stats import ttest_ind

result = ttest_ind(firsts["prglngth"], others["prglngth"])

result.pvalue

0.16755412639415004

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

Критерий корреляции

В разделе «Проверка значимости корреляции» главы 9 на с. 170 для проверки значимости корреляции между весом новорожденного и возрастом матери мы использовали критерий перестановок (permutation test) и обнаружили, что она статистически значима, а p-значение меньше 0.001.

Теперь то же самое можно проделать аналитически. Для этого нужно выполнить следующие шаги: сначала сгенерировать две выборки размером n из нормальных распределений, далее вычислить коэффициент корреляции Пирсона r, а затем преобразовать (transform) его с помощью функции ниже:

def transform_correlation(r, n):

return r * np.sqrt((n - 2) / (1 - r**2))

Преобразованный коэффициент корреляции следует t-распределению с параметром n-2. Чтобы посмотреть, как оно выглядит, нам потребуется следующая функция, которая генерирует некоррелированные выборки из стандартного нормального распределения:

def generate_data(n):

"""Некоррелированные последовательности из стандартного нормального

распределения."""

xs = np.random.normal(0, 1, n)

ys = np.random.normal(0, 1, n)

return xs, ys

А следующая функция вычисляет их корреляцию:

def correlation(data):

xs, ys = data

return np.corrcoef(xs, ys)[0, 1]

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

n = 100

rs = [correlation(generate_data(n)) for i in range(1001)]

Далее вычислим преобразованные коэффициенты корреляции:

ts = transform_correlation(np.array(rs), n)

Чтобы проверить, соответствуют ли значения в массиве ts t-распределению, воспользуемся следующей функцией. В ней создается объект, представляющий ИФР t-распределения:

from scipy.stats import t as student_t

def make_student_cdf(df):

"""Расчет ИФР для t-распределения Стьюдента."""

ts = np.linspace(-3, 3, 101)

ps = student_t.cdf(ts, df=df)

return Cdf(ps, ts)

Параметр t-распределения обозначается как df, то есть «степени свободы» (degrees of freedom). На следующем рисунке показана ИФР t-распределения с параметром n-2, а также эмпирическая ИФР преобразованных коэффициентов корреляции:

make_student_cdf(df=n — 2).plot(**model_options)

cdf_ts = Cdf.from_seq(ts)

cdf_ts.plot(метка="случайные выборки")

decorate(xlabel="Преобразованный коэффициент корреляции", ylabel="ИФР")

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

Если извлечь выборки из других распределений, то их преобразованные коэффициенты корреляции не будут в точности соответствовать t-распределению, но будут сходиться к t-распределению по мере увеличения размера выборки. Посмотрим, применимо ли это утверждение к корреляции возраста матери и веса новорожденного. Из объекта DataFrame данных о живых новорожденных выберем строки с действительными (valid) значениями:

valid = live.dropna(subset=["agepreg", "totalwgt_lb"])

n = len(valid)

9038

Фактический коэффициент корреляции составляет около 0,07:

data = valid["agepreg"].values, valid["totalwgt_lb"].values

r_actual = correlation(data)

r_actual

0.0688339703541091

Как вы узнали из раздела «Проверка значимости корреляции», нулевую гипотезу можно смоделировать при помощи перестановки выборочных значений:

def permute(data):

"""Перетасовка значений x."""

xs, ys = data

new_xs = xs.copy()

np.random.shuffle(new_xs)

return new_xs, ys

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

permuted_corrs = [correlation(permute(data)) for i in range(1001)]

Далее можно вычислить преобразованные корреляции:

ts = transform_correlation(np.array(permuted_corrs), n)

На следующем рисунке показана эмпирическая ИФР, полученная для последовательности ts, вместе с ИФР t-распределения с параметром n-2:

make_student_cdf(n - 2).plot(**model_options)

Cdf.from_seq(ts).plot(label="перемешанные данные")

decorate(xlabel="Преобразованный коэффициент корреляции", ylabel="ИФР")

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

t_actual = transform_correlation(r_actual, n)

Теперь используем ИФР t-распределения для вычисления вероятности того, что при верности нулевой гипотезы p-значение будет равно величине t_actual:

right = 1 - student_t.cdf(t_actual, df=n - 2)

right

2.861466619208386е-11

Также рассчитаем вероятность получить отрицательное значение, равное -t_actual:

left = student_t.cdf(-t_actual, df=n — 2)

left

2.8614735536574016е-11

Сумма этих двух вероятностей равна вероятности того, что коэффициент корреляции равен величине r_actual, независимо от знака:

left + right

5.722940172865787е-11

Библиотека SciPy предоставляет функцию, которая выполняет те же вычисления и возвращает p-значение наблюдаемой корреляции:

from scipy.stats import pearsonr

corr, p_value = pearsonr(*data)

p_value

5.722947107314431е-11

Результаты практически одинаковые.

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

Однако на практике это может и не требоваться. Как правило, p-значение меньше 0.001 говорит о том, что наблюдаемый эффект вряд ли является случайным. Но знать, насколько именно эта случайность маловероятна, вовсе не обязательно.

Критерий хи-квадрат

В разделе «Проверка пропорций» главы 9 на с. 173 мы проверяли, является ли игральная кость поддельной, отталкиваясь от множества наблюдаемых (observed) исходов:

from empiricaldist import Hist

qs = np.arange(1, 7)

freqs = [8, 9, 19, 5, 8, 11]

observed = Hist(freqs, qs)

observed.index.name = "исход"

observed

исход	частота
1	8
2	9
3	19
4	5
5	8
6	11

Сначала мы вычисляли ожидаемую (expected) частоту всех исходов (outcomes):

num_rolls = observed.sum()

outcomes = observed.qs

expected = Hist(num_rolls / 6, outcomes)

Затем использовали следующую функцию для вычисления статистики хи-квадрат:

def chi_squared_stat(observed, expected):

diffs = (observed - expected) ** 2

ratios = diffs / expected

return np.sum(ratios.values.flatten())

observed_chi2 = chi_squared_stat(observed, expected)

Статистика хи-квадрат широко используется для такого рода данных, поскольку их выборочное распределение при нулевой гипотезе сходится к распределению, которое легко вычислить, — неслучайно оно называется распределением хи-квадрат. Чтобы понять, как оно выглядит, воспользуемся следующей функцией, которая моделирует (симулирует) бросание неподдельной игральной кости:

def simulate_dice(observed):

n = np.sum(observed)

rolls = np.random.choice(observed.qs, n, replace=True)

hist = Hist.from_seq(rolls)

return hist

В цикле ниже моделирование и последующее вычисление статистики хи-квадрат результатов повторяются существенное количество раз:

simulated_chi_squared = [

chi_squared_stat(simulate_dice(observed), expected) for i in range(1001)

]

cdf_simulated = Cdf.from_seq(simulated_chi_squared)

Чтобы проверить, следуют ли полученные результаты распределению хи-квадрат, применим следующую функцию, вычисляющую ИФР распределения хи-квадрат с параметром df:

from scipy.stats import chi2 as chi2_dist

def chi_squared_cdf(df):

"""Дискретная аппроксимация ИФР распределения хи-квадрат."""

xs = np.linspace(0, 21, 101)

ps = chi2_dist.cdf(xs, df=df)

return Cdf(ps, xs)

При n возможных исходах статистика хи-квадрат результатов симуляции должна соответствовать распределению хи-квадрат с параметром n-1:

n = len(observed)

cdf_model = chi_squared_cdf(df=n - 1)

Вот эмпирическая ИФР статистики хи-квадрат результатов симуляции, а также ИФР распределения хи-квадрат:

cdf_model.plot(**model_options)

cdf_simulated.plot(label="симуляция")

decorate(xlabel="Статистика хи-квадрат", ylabel="ИФР")

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

p_value = 1 - chi2_dist.cdf(observed_chi2, df=n - 1)

p_value

0.04069938850404997

В библиотеке SciPy есть функция, которая выполняет эти же вычисления:

from scipy.stats import chisquare

chi2_stat, p_value = chisquare(f_obs=observed, f_exp=expected)

Результат очень близок к вычисленному нами p-значению:

p_value

0.040699388504049985

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

Вычислительные и аналитические методы

В этой книге основное внимание уделяется численным методам, таким как ресемплинг и перестановка. Они имеют ряд преимуществ перед аналитическими:

Простота объяснения и понимания

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

Надежность и универсальность

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

Возможность отладки

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

Но у них есть и свои недостатки:

• Они могут требовать больше вычислительных ресурсов.

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

Принимая во внимание все эти плюсы и минусы, я рекомендую придерживаться следующей процедуры:

1. Используйте вычислительные методы для разведочного анализа данных. Если ответ вас удовлетворяет, а время выполнения приемлемо, то на этом можно остановиться.

2. Если время выполнения оказывается неприемлемым, поищите способы оптимизации. Использование аналитических методов — один из таких возможных способов.

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

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

Глоссарий

График нормальной вероятности (normal probability plot)

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

t-критерий Стьюдента для независимых выборок (independent sample t-test)

Метод вычисления p-значения наблюдаемой разницы между средними значениями двух независимых групп.

t-распределение (t-distribution)

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

Распределение хи-квадрат (chi-squared distribution)

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

Статистика хи-квадрат (chi-squared statistic)

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

Упражнения

Упражнение 14.1

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

grouped = adelie.groupby("Sex")

lengths_male = grouped.get_group("MALE")["Flipper Length (mm)"]

lengths_female = grouped.get_group("FEMALE")["Flipper Length (mm)"]

observed_diff = lengths_male.mean() - lengths_female.mean()

observed_diff

4.616438356164366

С помощью функции sampling_dist_mean создайте объекты Normal, представляющие выборочные распределения средней длины крыльев в двух группах, обращая при этом внимание на несовпадение размеров групп. Затем вычислите выборочное распределение разности и 90%-ный доверительный интервал.

Упражнение 14.2

Пользуясь данными NSFG, мы вычисляли корреляцию между весом новорожденного и возрастом матери, а также с помощью t-распределения определяли p-значение. Теперь проделаем то же самое с весом новорожденного и возрастом отца, который указан в столбце hpagelb:

valid = live.dropna(subset=["hpagelb", "totalwgt_lb"])

n = len(valid)

8933

Наблюдаемая корреляция составляет около 0.065:

data = valid["hpagelb"].values, valid["totalwgt_lb"].values

r_actual = correlation(data)

r_actual

0.06468629895432174

Вычислите преобразованное значение коэффициента корреляции t_actual. Используйте ИФР t-распределения для вычисления p-значения — является ли эта корреляция статистически значимой? Проверьте полученные результаты с помощью функции pearsonr библиотеки SciPy.

Упражнение 14.3

В одном из упражнений главы 11 мы рассматривали гипотезу Триверса — Уилларда, которая предполагает, что у многих млекопитающих соотношение полов зависит от «состояния матери», то есть таких факторов, как возраст, физические параметры, состояние здоровья и социальный статус. Ряд исследований показал наличие этого эффекта у людей, но результаты неоднозначны.

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

Упражнение 14.4

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

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

До вмешательства студенты-мужчины выше оценили свое участие в задачах по программированию в рамках проектов, чем студентки-женщины: средний балл для мужчин составил 3.57 при стандартной ошибке (СО) 0.28, а средний балл для женщин — 1.91 при стандартной ошибке 0.32.

После вмешательства гендерный разрыв сократился: средний балл для мужчин составил 3.44 (СО 0.16), средний балл для женщин — 3.18 (СО 0.16).

1. Чтобы представить выборочные распределения оценок средних значений до и после вмешательства для студентов мужского и женского пола, создайте четыре объекта Normal. Поскольку у нас есть стандартные ошибки для оценок средних значений, то для получения параметров выборочных распределений нам не нужно знать размер выборки.

2. Рассчитайте выборочное распределение гендерного разрыва — разности средних баллов — до и после вмешательства.

3. Затем вычислите выборочное распределение разности разностей, то есть изменение размера разрыва. Вычислите 95%-ный доверительный интервал и p-значение.

Есть ли доказательства того, что после проведенного вмешательства гендерный разрыв сократился?