Глава 5. Моделирование распределений

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

Вы на примерах убедитесь в следующем:

• Количество попаданий и промахов в соревнованиях по стендовой стрельбе хорошо моделируется с помощью биномиального распределения.

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

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

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

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

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

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

Для моделирования воспользуемся следующей функцией, которая принимает в качестве аргументов количество целей (n) и вероятность попадания в каждую из них (p), а возвращает последовательность из единиц и нулей, обозначающих попадания и промахи соответственно:

def flip(n, p):

    choices = [1, 0]

    probs = [p, 1 - p]

    return np.random.choice(choices, n, p=probs)

Вот пример моделирования раунда из 25 мишеней, в котором вероятность попадания в каждую из них равна 90 %:

flip(25, 0.9)

array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,

      1, 1, 1])

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

from empiricaldist import Pmf

 

seq = flip(1000, 0.9)

pmf = Pmf.from_seq(seq)

pmf

	вероятность
0	0.101
1	0.899

Теперь для моделирования (simulate) раунда стендовой стрельбы реализуем новую функцию, вызывающую функцию flip и возвращающую количество попаданий:

def simulate_round(n, p):

    seq = flip(n, p)

    return seq.sum()

Предположим, что в крупном соревновании 200 участников стреляют по 5 раундов каждый с одинаковой вероятностью попадания в цель, p=0.9. Можем смоделировать подобное соревнование, вызвав simulate_round 1000 раз:

n = 25

р = 0.9

results_sim = [simulate_round(n, p) for i in range(1000)]

Средний балл близок к 22.5 — произведению n на p:

np.mean(results_sim), n * p

(22.522, 22.5)

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

from empiricaldist import Pmf

 

pmf_sim = Pmf.from_seq(results_sim, name="результаты моделирования")

 

pmf_sim.bar()

decorate(xlabel="Попадания", ylabel="Вероятность")

Максимум распределения близок к среднему значению, а само распределение скошено влево.

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

from scipy.special import comb

 

def binomial_pmf(k, n, p):

    return comb(n, k) * (p**k) * ((1 - p) ** (n - k))

В библиотеке SciPy реализована функция comb, которая вычисляет количество сочетаний (combinations) из n элементов, взятых по k за раз, что часто произносится как «из n по k».

Функция binomial_pmf вычисляет вероятность при заданном p, попасть k раз, сделав n попыток. Если вызвать эту функцию, подав на вход массив из k значений, то можно создать объект Pmf, представляющий собой распределение вероятностей исходов:

ks = np.arange(16, n + 1)

ps = binomial_pmf(ks, n, p)

pmf_binom = Pmf(ps, ks, name="биномиальная модель")

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

from thinkstats import two_bar_plots

 

two_bar_plots(pmf_sim, pmf_binom)

decorate(xlabel="Попадания", ylabel="Вероятность")

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

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

filename = "Shooting_at_the_2020_Summer_Olympics_Mens_skeet"

tables = pd.read_html(filename)

table = tables[6]

table.head()

	Rank (Место)	Athlete (Спортсмен)	Country (Страна)	1	2	3	4	5	Total[3] (Итого)	Shoot-off (Перестрелка)	Notes (Примечания)
0	1	Éric Delaunay	France	25	25	25	24	25	124	+6	Q, OR
1	2	Tammaro Cassandro	Italy	24	25	25	25	25	124	+5	Q, OR
2	3	Eetu Kallioinen	Finland	25	25	24	25	24	123	NaN	Q
3	4	Vincent Hancock	United States	25	25	25	25	22	122	+8	Q
4	5	Abdullah Al-Rashidi	Kuwait	25	25	24	25	23	122	+7	Q

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

columns = ["1", "2", "3", "4", "5"]

results = table[columns].values.flatten()

Для 30 участников мы получили результаты 150 раундов, по 25 выстрелов в каждом — всего 3574 попаданий (hits) при 3750 попытках (shots):

total_shots = 25 * len(results)

total_hits = results.sum()

n, total_shots, total_hits

(25, 3750, 3575)

Таким образом, доля успешных попыток составляет 95,3 %:

p = total_hits / total_shots

p

0.9533333333333334

Теперь создадим объект Pmf для биномиального распределения с n = 25 и только что вычисленным значением p:

ps = binomial_pmf(ks, n, p)

pmf_binom = Pmf(ps, ks, name="биномиальная модель")

Полученное распределение можно сравнить с Pmf фактических результатов:

pmf_results = Pmf.from_seq(results, name="фактические результаты")

 

two_bar_plots(pmf_results, pmf_binom)

decorate(xlabel="Попадания", ylabel="Вероятность")

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

Распределение Пуассона

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

Начнем с моделирования игры продолжительностью 60 минут (3600 секунд), исходя из предположения, что команды забивают в среднем по 6 шайб за игру и что вероятность заброшенной шайбы в секунду, p, постоянна во времени:

n = 3600

m = 6

p = m / 3600

p

0.0016666666666666668

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

def simulate_goals(n, p):

    return flip(n, p).sum()

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

goals = [simulate_goals(n, p) for i in range(1001)]

np.mean(goals)

6.021978021978022

Для описания этих результатов можно было бы использовать биномиальное распределение, но, когда n велико, а p мало, результаты также хорошо аппроксимируются распределением Пуассона. Оно определяется параметром, обычно обозначаемым греческой буквой λ («лямбда»). В коде ниже этому параметру соответствует переменная lam (lambda нельзя использовать в качестве имени переменной в Python, поскольку это зарезервированное слово). lam в данном примере представляет собой показатель результативности, составляющий 6 голов за игру.

ФВ распределения Пуассона легко вычислить. Имея значение lam, для вычисления вероятности того, что в игре будет забито k голов, можно использовать следующую функцию:

from scipy.special import factorial

 

def poisson_pmf(k, lam):

    return (lam**k) * np.exp(-lam) / factorial(k)

В библиотеке SciPy есть функция factorial, которая вычисляет произведение целых чисел от 1 до k.

Если вызвать poisson_pmf, указав массив k значений, то можно построить объект Pmf, представляющий распределение исходов:

lam = 6

ks = np.arange(20)

ps = poisson_pmf(ks, lam)

pmf_poisson = Pmf(ps, ks, name="пуассоновская модель")

А также легко подтвердить, что среднее значение распределения близко к 6:

pmf_poisson.normalize()

pmf_poisson.mean()

5.999925498375129

Теперь сравним результаты моделирования с распределением Пуассона с тем же средним:

pmf_sim = Pmf.from_seq(goals, name="моделирование")

 

two_bar_plots(pmf_sim, pmf_poisson)

decorate(xlabel="Шайбы", ylabel="Вероятность")

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

Сначала я скачал результаты всех игр регулярного чемпионата Национальной хоккейной лиги (НХЛ) 2023–2024 годов (без игр плей-офф) с сайта HockeyReference. Далее я извлек информацию о шайбах, забитых за 60 минут основного времени матча, без учета дополнительного времени и серии буллитов. Результаты записаны в файле формата HDF, в котором один ключ соответствует одной игре. По ключу можно получить список моментов времени, в секундах с начала игры, в которые были заброшены шайбы. Инструкции по загрузке данных приведены в Jupyter-блокноте для этой главы.

Вот как считывать ключи из файла:

filename = "nhl_2023_2024.hdf"

 

with pd.HDFStore(filename, "r") as store:

    keys = store.keys()

 

len(keys), keys[0]

(1312, '/202310100PIT')

В регулярном сезоне было проведено 1312 игр. Каждый ключ содержит дату игры и аббревиатуру из трех букв, обозначающую принимающую команду. Функция read_hdf по ключу выдает список моментов времени для забитых шайб:

times = pd.read_hdf(filename, key=keys[0])

times

0     424

1    1916

2    2137

3    3005

4    3329

5    3513

dtype: int64

В первой игре сезона было забито шесть шайб: первая — после 424 секунд игры, последняя — после 3513 секунд, всего за 87 секунд до конца игры.

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

goals = []

 

for key in keys:

    times = pd.read_hdf(filename, key=key)

    n = len(times)

    goals.append(n)

Среднее количество шайб за игру немного превышает 6:

lam = np.mean(goals)

lam

6.0182926829268295

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

ps = poisson_pmf(ks, lam)

pmf_poisson = Pmf(ps, ks, name="пуассоновская модель")

И вот как это выглядит в сравнении с ФВ данных:

pmf_goals = Pmf.from_seq(goals, name="забитые шайбы")

 

two_bar_plots(pmf_goals, pmf_poisson)

decorate(xlabel="Шайбы", ylabel="Вероятность")

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

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

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

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

n = 3600

m = 6

p = m / 3600

p

0.0016666666666666668

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

def simulate_first_goal(n, p):

    return flip(n, p).argmax()

Мы получаем верный результат, поскольку функция flip возвращает последовательность из единиц и нулей, а значит, максимальное значение почти всегда равно 1. Если в последовательности имеется одна единица (забитая шайба) или больше, то argmax возвращает индекс первой из них. В противном случае возвращается 0, но такое бывает редко, так что мы проигнорируем этот вариант.

При помощи функции simulate_first_goal смоделируем 1001 игру и составим список моментов времени первых забитых шайб:

first_goal_times = [simulate_first_goal(n, p) for i in range(1001)]

mean = np.mean(first_goal_times)

mean

597.7902097902098

Среднее время до первой шайбы составляет около 600 секунд, или 10 минут. И это логично: если мы ожидаем 6 шайб за 60 минут игры, то в среднем можно ожидать одну забитую шайбу каждые 10 минут.

Когда n велико, а p мало, можно математически доказать, что время ожидания первого гола следует экспоненциальному распределению.

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

def exponential_cdf(x, lam):

    return 1 - np.exp(-lam * x)

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

lam = 1 / mean

lam

0.0016728276636563566

Если подать на вход этой функции массив значений времени, можно получить аппроксимацию распределения времени первой забитой шайбы. NumPy-функция linspace создает массив равноотстоящих значений; в данном примере она вычисляет 201 значение от 0 до 3600, включая границы диапазона:

from empiricaldist import Cdf

 

ts = np.linspace(0, 3600, 201)

ps = exponential_cdf(ts, lam)

cdf_expo = Cdf(ps, ts, name="экспоненциальная модель")

Теперь можно сравнить результаты моделирования с только что вычисленным экспоненциальным распределением:

cdf_sim = Cdf.from_seq(first_goal_times, name="моделирование")

 

cdf_expo.plot(ls=":", color="gray")

cdf_sim.plot()

 

decorate(xlabel="Время первой шайбы (с)", ylabel="ИФР")

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

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

firsts = []

 

for key in keys:

    times = pd.read_hdf(filename, key=key)

    if len(times) > 0:

        firsts.append(times[0])

    else:

        firsts.append(np.nan)

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

lam = 1 / np.nanmean(firsts)

lam

0.0015121567467720825

Теперь можно построить ИФР экспоненциального распределения с тем же голевым темпом, что и в данных:

ps = exponential_cdf(ts, lam)

cdf_expo = Cdf(ps, ts, name="экспоненциальная модель")

Чтобы посчитать ИФР данных, применим аргумент dropna=False, что приводит к добавлению nan-значений на конце временного диапазона:

cdf_firsts = Cdf.from_seq(firsts, name="данные", dropna=False)

cdf_firsts.tail()

	вероятность
3286.0	0.996951
3581.0	0.997713
NaN	1.000000

На следующем рисунке можно сравнить экспоненциальное распределение с распределением реальных данных:

cdf_expo.plot(ls=":", color="gray")

cdf_firsts.plot()

 

decorate(xlabel="Время первой шайбы (с)", ylabel="ИФР")

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

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

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

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

Для имитации такой модели в течение n дней можно задействовать следующую функцию, возвращающую итоговое значение прироста (gains) веса:

def simulate_growth(n):

    choices = [1, 2, 3]

    gains = np.random.choice(choices, n)

    return gains.sum()

В модуле random библиотеки NumPy имеется функция choice, которая генерирует массив n результатов случайного выбора из последовательности значений — в данном примере из списка choices.

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

sim_weights = [simulate_growth(100) for i in range(1001)]

m, s = np.mean(sim_weights), np.std(sim_weights)

m, s

(199.37062937062936, 8.388630840376777)

Среднее значение близко к 200 фунтам, а стандартное отклонение составляет около 8 фунтов. Чтобы увидеть, соответствуют ли веса нормальному распределению, воспользуемся следующей функцией, которая принимает на вход выборку и создает объект Cdf, представляющий нормальное распределение с теми же значениями среднего и стандартного отклонения, что и выборка. Распределение рассчитывается в диапазоне от четырех стандартных отклонений ниже среднего (low) до четырех стандартных отклонений выше среднего (high):

from scipy.stats import norm

 

def make_normal_model(data):

    m, s = np.mean(data), np.std(data)

    low, high = m - 4 * s, m + 4 * s

    qs = np.linspace(low, high, 201)

    ps = norm.cdf(qs, m, s)

    return Cdf(ps, qs, name="нормальная модель")

Вот пример ее использования:

cdf_model = make_normal_model(sim_weights)

Теперь создадим объект Cdf, содержащий распределение результатов моделирования:

cdf_sim_weights = Cdf.from_seq(sim_weights, name="моделирование")

Напишем функцию для сравнения распределений. Аргументы cdf_model и cdf_data являются объектами Cdf. xlabel — это строка, а options — это словарь параметров отображения cdf_data:

def two_cdf_plots(cdf_model, cdf_data, xlabel="", **options):

    cdf_model.plot(ls=":", color="gray")

    cdf_data.plot(**options)

    decorate(xlabel=xlabel, ylabel="ИФР")

А вот результат ее работы:

two_cdf_plots(cdf_model, cdf_sim_weights, xlabel="Вес (фунты)")

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

Но сначала посмотрим, насколько хорошо нормальное распределение описывает реальные данные. Для примера рассмотрим распределение веса новорожденных по данным Национального исследования роста семьи (NSFG). Воспользуемся функцией read_fem_preg для чтения данных, а затем выберем столбец totalwgt_lb, в котором записан вес при рождении в фунтах:

import nsfg

 

preg = nsfg.read_fem_preg()

birth_weights = preg["totalwgt_lb"].dropna()

Средний вес новорожденного составляет около 7.27 фунта, а стандартное отклонение — 1.4 фунта, но, как вы уже знаете, в этом датасете есть выбросы, которые, скорее всего, являются ошибочными значениями:

m, s = np.mean(birth_weights), np.std(birth_weights)

m, s

(7.265628457623368, 1.40821553384062)

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

from scipy.stats import trimboth

 

trimmed = trimboth(birth_weights, 0.01)

m, s = np.mean(trimmed), np.std(trimmed)

m, s

(7.280883100022579, 1.2430657948614345)

Среднее значение таких усеченных (trimmed) данных немного меньше, а стандартное отклонение существенно меньше. Воспользуемся этими данными для построения нормальной модели:

cdf_model = make_normal_model(trimmed)

И сравним ее с Cdf реальных данных:

cdf_birth_weight = Cdf.from_seq(birth_weights, name='данные')

 

two_cdf_plots(cdf_model, cdf_birth_weight, xlabel="Вес при рождении (фунты)")

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

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

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

Следующая функция моделирует процесс такого пропорционального роста, при котором тыква прибавляет в весе 3 % в плохую погоду, 5 % — в удовлетворительную и 7 % — в хорошую. Снова будем исходить из предположения, что погода в любой день с равной вероятностью будет плохой, удовлетворительной или хорошей:

def simulate_proportionate_growth(n):

    choices = [1.03, 1.05, 1.07]

    gains = np.random.choice(choices, n)

    return gains.prod()

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

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

sim_weights = [simulate_proportionate_growth(100) for i in range(1001)]

np.mean(sim_weights), np.std(sim_weights)

(130.80183363824722, 20.956047434921466)

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

Можно доказать математически, что вес в данном случае соответствует логарифмически нормальному (логнормальному) распределению. Это означает, что логарифмы веса соответствуют нормальному распределению. Для проверки вычислим логарифм веса, его среднее значение и стандартное отклонение. Можно было бы применить логарифм по любому основанию, но лучше использовать основание 10, поскольку в этом случае будет проще интерпретировать результаты:

log_sim_weights = np.log10(sim_weights)

m, s = np.mean(log_sim_weights), np.std(log_sim_weights)

m, s

(2.1111299372609933, 0.06898607064749827)

Теперь сравним распределение логарифма веса с нормальным распределением с теми же средним значением и стандартным отклонением:

cdf_model = make_normal_model(log_sim_weights)

cdf_log_sim_weights = Cdf.from_seq(log_sim_weights, name="моделирование")

 

two_cdf_plots(

    cdf_model, cdf_log_sim_weights, xlabel="Вес тыквы (log10 фунтов)"

)

Модель, как и ожидалось, прекрасно согласуется с результатами моделирования.

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

Национальный центр профилактики хронических заболеваний и укрепления здоровья (The National Center for Chronic Disease Prevention and Health Promotion) США проводит ежегодное исследование в рамках Системы наблюдения за поведенческими факторами риска (Behavioral Risk Factor Surveillance System, BRFSS). В 2008 году было опрошено 414 509 респондентов, которым были заданы вопросы об их демографических характеристиках, состоянии здоровья и факторах риска. Среди собранных ими данных — данные о весе в килограммах 398 484 респондентов. Инструкции по загрузке этих данных приведены в Jupyter-блокноте к этой главе.

В модуле thinkstats реализована функция, считывающая данные BRFSS и возвращающая Pandas-объект DataFrame:

from thinkstats import read_brfss

 

brfss = read_brfss()

Вес взрослых в килограммах записан в столбце wtkg2:

adult_weights = brfss["wtkg2"].dropna()

m, s = np.mean(adult_weights), np.std(adult_weights)

m, s

(78.9924529968581, 19.546132387397257)

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

cdf_model = make_normal_model(adult_weights)

cdf_adult_weights = Cdf.from_seq(adult_weights, name="вес взрослого")

 

two_cdf_plots(cdf_model, cdf_adult_weights, xlabel="Вес взрослого (кг)")

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

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

log_adult_weights = np.log10(adult_weights)

cdf_model = make_normal_model(log_adult_weights)

 

cdf_log_adult_weights = Cdf.from_seq(log_adult_weights, name="логарифм веса

взрослого")

two_cdf_plots(cdf_model, cdf_log_adult_weights, xlabel="Вес взрослого (log10 кг)")

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

Зачем нужны модели

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

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

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

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

Кроме того, теоретические распределения, как вы увидите в главе 14, поддаются математическому анализу.

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

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

Глоссарий

Биномиальное распределение (binomial distribution)

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

Распределение Пуассона (Poisson distribution)

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

Экспоненциальное распределение (exponential distribution)

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

Нормальное распределение (normal distribution)

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

Логнормальное распределение (lognormal distribution)

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

Упражнения

Упражнение 5.1

В файле респондента датасета NSFG в столбце numfmhh указано «количество членов семьи» в домохозяйстве каждой респондентки. Вот как можно задействовать функцию read_fem_resp для чтения файла и метод query для выбора респонденток, которым на момент опроса было 25 лет и старше:

from nsfg import read_fem_resp

 

resp = read_fem_resp()

older = resp.query("age >= 25")

num_family = older["numfmhh"]

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

Упражнение 5.2

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

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

intervals = []

 

for key in keys:

    times = pd.read_hdf(filename, key=key)

    if len(times) > 1:

        intervals.extend(times.diff().dropna())

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

Упражнение 5.3

Распределение человеческого роста больше похоже на нормальное или логнормальное? Чтобы получить ответ на этот вопрос, выберем из BRFSS данные о росте следующим образом:

adult_heights = brfss["htm3"].dropna()

m, s = np.mean(adult_heights), np.std(adult_heights)

m, s

(168.82518961012298, 10.35264015645592)

Постройте ИФР этих значений и сравните ее с нормальным распределением с теми же средним значением и стандартным отклонением. Затем вычислите логарифмы роста и постройте их распределение в сравнении с нормальным распределением. По результатам визуального сравнения, какая модель лучше согласуется с данными?