Привет, Хаброжители!

Мы открыли предзаказ на книгу «Думай как аналитик. Статистика и данные с примерами на Python. 3-е изд.», хотим немного рассказать вам о ней и поделиться интересным отрывком.

Многие бросают Data Science из-за академической статистики с её громоздкими формулами. Но если вы уже знаете Python, у вас есть огромное преимущество. Код — это всё, что нужно для понимания работы с данными.

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

Что внутри?

• Реальный опыт: Весь цикл разведочного анализа (EDA) от очистки датасетов до проверки гипотез разбирается на практических примерах.

• Актуальный стек: Быстрое погружение в Pandas, NumPy, SciPy и визуализацию данных.

• Интерактивность: Все главы доступны в виде Jupyter-блокнотов — можно читать, сразу запускать код и делать упражнения в одном окне.

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

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

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

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

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

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

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

В качестве первого примера рассмотрим такой вид спорта, как стендовая стрельба, в котором участники используют гладкоствольные ружья для стрельбы по глиняным тарелочкам, подбрасываемым в воздух. На международных соревнованиях, в том числе Олимпийских играх, проводится пять раундов по 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

Теперь для моделирования (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()

В таблице содержится по одной строке для каждого участника и по одному столбцу для каждого из пяти раундов. Выберем столбцы с результатами раундов и с помощью 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="Вероятность")

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

Оформить предзаказ на книгу: «Думай как аналитик. Статистика и данные с примерами на Python. 3-е изд.» можно на нашем сайте.

Скидка 35% по промокоду - Предзаказ