Привет, Хаброжители!
Мы открыли предзаказ на книгу «Думай как аналитик. Статистика и данные с примерами на 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
вероятность | |
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="Вероятность")
Распределение Пуассона хорошо согласуется с данными, что означает, что эта модель точно описывает процесс забивания шайб в хоккее.
Оформить предзаказ на книгу: «Думай как аналитик. Статистика и данные с примерами на Python. 3-е изд.» можно на нашем сайте.
Скидка 35% по промокоду - Предзаказ
