Глава 4. Интегральная функция распределения

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

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

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

Если вы когда-либо проходили стандартизированный тест, то, вероятно, получали результаты в виде первичного балла (raw score) и процентильного ранга. В данном контексте процентильный ранг — это процент тестируемых, которые получили тот же балл, что и вы, или ниже. То есть нахождение «в 90-м про­центиле» означает, что вы сдали экзамен так же или лучше, чем 90 % испы­туемых.

Чтобы разобраться в процентилях и процентильных рангах, рассмотрим пример со скоростью бега. Несколько лет назад я участвовал в забеге James Joyce Ramble — состязании в беге на 10 километров, проходящем в Массачусетсе. После окончания забега я скачал его результаты, чтобы посмотреть на свое время в сравнении с временем других бегунов.

Инструкции по загрузке этих данных можно найти в Jupyter-блокноте для главы 4. В модуле relay.py реализована функция, которая считывает результаты и возвращает Pandas-объект DataFrame:

from relay import read_results

 

results = read_results()

results.head()

	Place	Div/Tot	Division	Guntime	Nettime	Min/Mile	MPH
0	1	1/362	М2039	30:43	30:42	4:57	12.121212
1	2	2/362	М2039	31:36	31:36	5:06	11.764706
2	3	3/362	М2039	31:42	31:42	5:07	11.726384
3	4	4/362	М2039	32:28	32:27	5:14	11.464968
4	5	5/362	М2039	32:52	32:52	5:18	11.320755

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

speeds = results["MPH"].values

Я финишировал с чистым временем (net time) 42:44, поэтому мою строку можно найти так:

my_result = results.query("Nettime == '42:44'")

my_result

	Place	Div/Tot	Division	Guntime	Nettime	Min/Mile	MPH
96	97	26/256	М4049	42:48	42:44	6:53	8.716707

Индекс моей строки равен 96, так что мою скорость можно получить так:

my_speed = speeds[96]

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

(speeds <= my_speed).sum()

1537

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

(speeds <= my_speed).mean() * 100

94.12124923453766

Результат, составивший около 94 %, — это мой процентильный ранг в забеге.

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

def percentile_rank(x, seq):

    """Процентильный ранг x.

 

    x: значение

    seq: последовательность значений

 

    returns: процентильный ранг 0–100

    """

    return (seq <= x).mean() * 100

В столбце Division датафрейма results указана категория каждого участника, определяемая по полу и возрасту. Например, моя категория была M4049, составленная из бегунов-мужчин в возрасте от 40 до 49 лет. При помощи метода query выберем строки бегунов моей категории и извлечем из них значения скорости:

my_division = results.query("Division == 'M4049'")

my_division_speeds = my_division["MPH"].values

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

percentile_rank(my_speed, my_division_speeds)

90.234375

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

def percentile(p, seq):

    n = len(seq)

    i = (1 - p / 100) * (n + 1)

    return seq[round(i)]

Здесь n — количество элементов в последовательности, а i — индекс элемента с заданным процентильным рангом. Когда мы указываем процентильный ранг, то соответствующее значение в последовательности называется процентилем:

percentile(90, my_division_speeds)

8.591885441527447

В моей категории 90-й процентиль составлял около 8.6 мили в час.

Прошло несколько лет с момента моего участия в том забеге, и я нахожусь в категории М5059. Поэтому посмотрим, с какой скоростью я должен был бы бежать, чтобы иметь тот же процентильный ранг в моей новой категории. Мы можем ответить на этот вопрос, переведя мой процентильный ранг в категории M4049, равный примерно 90.2 %, в значение скорости в категории M5059:

next_division = results.query("Division == 'M5059'")

next_division_speeds = next_division["MPH"].values

 

percentile(90.2, next_division_speeds)

8.017817371937639

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

next_division.query("MPH > 8.01").tail(1)

	Place	Div/Tot	Division	Guntime	Nettime	Min/Mile	MPH
222	223	18/171	М5059	46:30	46:25	7:29	8.017817

Этот бегун финишировал со временем 46:25 и занял 18-е место из 171 участника своей категории.

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

Интегральная функция распределения

Интегральная функция распределения, ИФР (cumulative distribution function, CDF), — это еще один способ, описывающий распределение набора значений, помимо таблицы частот и функции вероятности. При заданном значении x ИФР показывает долю значений, меньших или равных x. Для начала посмотрим на пример с короткой последовательностью:

t = [1, 2, 2, 3, 5]

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

from empiricaldist import Pmf

 

pmf = Pmf.from_seq(t)

pmf

	вероятность
1	0.2
2	0.4
3	0.2
5	0.2

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

pmf[2]

0.4

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

У Pmf есть метод make_cdf, который вычисляет накопленную сумму вероятностей:

cdf = pmf.make_cdf()

cdf

	вероятность
1	0.2
2	0.6
3	0.8
5	1.0

Результатом является объект Cdf, производный от Series библиотеки Pandas. Поэтому для поиска значения можно использовать оператор квадратных скобок:

cdf[2]

0.6000000000000001

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

Объект Cdf можно вызывать как функцию, используя круглые скобки:

cdf(3)

array(0.8)

Интегральная функция распределения определена для всех чисел, а не только для тех, которые присутствуют в последовательности:

cdf(4)

array(0.8)

Чтобы визуализировать объект Cdf, можно использовать метод step, который отображает Cdf в виде ступенчатой функции:

cdf.step()

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

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

from empiricaldist import Cdf

 

cdf_speeds = Cdf.from_seq(speeds)

А вот как выглядит распределение (вертикальная линия обозначает мою скорость):

cdf_speeds.step()

plt.axvline(my_speed, ls=":", color="gray")

decorate(xlabel="Скорость (миль/ч)", ylabel="ИФР")

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

cdf_speeds(my_speed) * 100

94.12124923453766

Таким образом, Cdf можно считать объектом, который на основе заданного значения вычисляет процентильный ранг, только возвращая долю значений в диапазоне от 0 до 1, а не процент от 0 до 100.

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

Например, если кто-то утверждает, что он пробежал так же быстро или быстрее, чем 50 % участников, можно определить его скорость следующим образом:

cdf_speeds.inverse(0.5)

array(6.70391061)

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

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

Сравнение ИФР

ИФР особенно эффективны для сравнения распределений. В качестве примера сравним распределение веса при рождении для первенцев и всех последующих детей. Мы снова загрузим датасет NSFG и создадим на его базе три объекта DataFrame: все дети, родившиеся живыми (live), первые (firsts) и последующие (others) дети:

from nsfg import get_nsfg_groups

 

live, firsts, others = get_nsfg_groups()

Из датафреймов firsts и others выберем полный вес при рождении в фунтах и воспользуемся методом dropna для удаления nan-значений:

first_weights = firsts["totalwgt_lb"].dropna()

first_weights.mean()

7.201094430437772

other_weights = others["totalwgt_lb"].dropna()

other_weights.mean()

7.325855614973262

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

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

from empiricaldist import Pmf

 

first_pmf = Pmf.from_seq(first_weights, name="первые")

other_pmf = Pmf.from_seq(other_weights, name="последующие")

Но, как видим, здесь они не очень наглядны:

from thinkstats import two_bar_plots

 

two_bar_plots(first_pmf, other_pmf, width=0.06)

decorate(xlabel="Вес (фунты)", ylabel="Вероятность")

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

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

Хорошей альтернативой такому подходу станет график ИФР:

first_cdf = first_pmf.make_cdf()

other_cdf = other_pmf.make_cdf()

Вот как он выглядит:

first_cdf.plot(ls="--")

other_cdf.plot(alpha=0.5)

decorate(xlabel="Вес (фунты)", ylabel="ИФР")

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

Процентильная статистика

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

Альтернатива — использование статистических показателей, основанных на процентилях распределения, которые, как правило, более робастны. Это означает, что они менее чувствительны к выбросам. Чтобы убедиться в этом, снова загрузим данные NSFG и не будем проводить их очистку:

from nsfg import read_stata

 

dct_file = "2002FemPreg.dct"

dat_file = "2002FemPreg.dat.gz"

 

preg = read_stata(dct_file, dat_file)

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

birthwgt_lb = preg["birthwgt_lb"]

birthwgt_oz = preg["birthwgt_oz"]

Если мы создадим объект Hist со значениями из birthwgt_oz, то увидим, что они содержат специальные значения 97, 98 и 99, указывающие на недостающие данные:

from empiricaldist import Hist

 

Hist.from_seq(birthwgt_oz).tail(5)

birthwgt_oz	freqs
14.0	475
15.0	378
97.0	1
98.0	1
99.0	46

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

Hist.from_seq(birthwgt_lb).tail(5)

birthwgt_lb	freqs
15.0	1
51.0	1
97.0	1
98.0	1
99.0	57

Теперь представьте два сценария. В первом мы очищаем эти переменные, заменяя пропущенные и недопустимые значения на nan, а затем вычисляем общий вес в фунтах. Делением объекта birthwgt_oz_clean на 16 мы переводим вес из унций в фунты в десятичной системе счисления:

birthwgt_lb_clean = birthwgt_lb.replace([51, 97, 98, 99], np.nan)

birthwgt_oz_clean = birthwgt_oz.replace([97, 98, 99], np.nan)

 

total_weight_clean = birthwgt_lb_clean + birthwgt_oz_clean / 16

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

total_weight_bogus = birthwgt_lb + birthwgt_oz / 16

Такой «фиктивный» датасет содержит только 49 некорректных значений, что составляет около 0.5 % от общего объема данных:

count1, count2 = total_weight_bogus.count(), total_weight_clean.count()

diff = count1 - count2

 

diff, diff / count2 * 100

(49, 0.5421553441026776)

Теперь вычислим среднее значение в обоих сценариях:

mean1, mean2 = total_weight_bogus.mean(), total_weight_clean.mean()

mean1, mean2

(7.319680587652691, 7.265628457623368)

Фиктивные значения оказывают умеренное влияние на среднее. Если мы примем среднюю величину очищенных данных за верную, то среднее фиктивных данных будет отличаться менее чем на 1 %:

(mean1 - mean2) / mean2 * 100

0.74394294099376

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

std1, std2 = total_weight_bogus.std(), total_weight_clean.std()

std1, std2

(2.096001779161835, 1.4082934455690173)

(std1 - std2) / std2 * 100

48.83274403900607

Стандартное отклонение фиктивных данных отличается почти на 50 %, что намного заметнее. Наконец, вот асимметрия двух датасетов:

def skewness(seq):

    deviations = seq - seq.mean()

    return np.mean(deviations**3) / seq.std(ddof=0) ** 3

 

skew1, skew2 = skewness(total_weight_bogus), skewness(total_weight_clean)

skew1, skew2

(22.251846195422484, -0.5895062687577697)

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

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

Другой способ — использование показателей, основанных на процентилях. К ним относятся:

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

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

• Квартильная асимметрия использует квартили распределения (25-й, 50-й и 75-й процентили) для количественной оценки асимметрии.

Объект Cdf удобен для вычисления таких процентильных показателей. В демонстрационных целях создадим объекты Cdf из фиктивного (bogus) и очищенного (clean) датасетов:

cdf_total_weight_bogus = Cdf.from_seq(total_weight_bogus)

cdf_total_weight_clean = Cdf.from_seq(total_weight_clean)

Следующая функция принимает на вход объект Cdf и использует его метод inverse для вычисления 50-го процентиля, который является медианой (во всяком случае, это один из способов определения медианы датасета):

def median(cdf):

    m = cdf.inverse(0.5)

    return m

Теперь вычислим медиану обоих датасетов:

median(cdf_total_weight_bogus), median(cdf_total_weight_clean)

(array(7.375), array(7.375))

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

Межквартильный размах, МКР (interquartile range, IQR), — это разница между 75-м и 25-м процентилями. Следующая функция принимает объект Cdf и возвращает значение МКР:

def iqr(cdf):

    low, high = cdf.inverse([0.25, 0.75])

    return high - low

А вот значения межквартильного размаха для двух датасетов:

iqr(cdf_total_weight_bogus), iqr(cdf_total_weight_clean)

(1.625, 1.625)

Как правило, выбросы оказывают меньшее влияние на МКР, чем на стандартное отклонение (в данном случае вообще не оказывают).

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

• медианы (median);

• средней точки (midpoint) между 25-м и 75-м процентилями;

• половины МКР (semi-IQR).

def quartile_skewness(cdf):

    low, median, high = cdf.inverse([0.25, 0.5, 0.75])

    midpoint = (high + low) / 2

    semi_iqr = (high - low) / 2

    return (midpoint - median) / semi_iqr

Вот чему равна квартильная асимметрия для двух датасетов:

qskew1 = quartile_skewness(cdf_total_weight_bogus)

qskew2 = quartile_skewness(cdf_total_weight_clean)

qskew1, qskew2

(-0.07692307692307693, -0.07692307692307693)

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

Случайные числа

Объекты Cdf обеспечивают эффективный способ генерации случайных чисел из распределений. Сначала мы генерируем случайные числа из равномерного распределения от 0 до 1. Затем в полученных точках вычисляем обратную ИФР. Следующая функция реализует этот алгоритм:

def sample_from_cdf(cdf, n):

    ps = np.random.random(size=n)

    return cdf.inverse(ps)

Для примера сгенерируем случайную выборку (sample) скоростей бега:

sample = sample_from_cdf(cdf_speeds, 1001)

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

cdf_sample = Cdf.from_seq(sample)

 

cdf_speeds.plot(label="оригинал", ls="--")

cdf_sample.plot(label="выборка", alpha=0.5)

 

decorate(xlabel="Скорость (миль/ч)", ylabel="ИФР")

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

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

percentile_ranks = cdf_speeds(sample) * 100

А вот ИФР процентильных рангов:

cdf_percentile_rank = Cdf.from_seq(percentile_ranks)

cdf_percentile_rank.plot()

 

decorate(xlabel="Процентильный ранг", ylabel="ИФР")

ИФР процентильных рангов близка к прямой линии, проведенной между 0 и 1. И это логично: ведь в любом распределении доля значений с процентильным рангом менее 50 % равна 0.5, с рангом менее 90 % — 0.9 и т.д.

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

sample = cdf_speeds.sample(1001)

Глоссарий

Процентильный ранг (percentile rank)

Процент значений в распределении, которые меньше или равны заданной величине.

Процентиль (percentile)

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

Интегральная функция распределения, ИФР (cumulative distribution function, CDF)

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

Квантиль (quantile)

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

Робастный (robust)

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

Межквартильный разброс, МКР (interquartile range, IQR)

Разница между 75-м и 25-м процентилями, используемая для оценки разброса распределения.

Упражнения

Упражнение 4.1

Сколько вы весили при рождении? Если вы не знаете, позвоните своей маме или кому-то еще, кто знает. А если такого человека нет, то для этого упражнения вы можете использовать мой вес при рождении — 8.5 фунта.

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

from nsfg import get_nsfg_groups

 

live, firsts, others = get_nsfg_groups()

Упражнение 4.2

В датасете NSFG столбец babysex указывает пол рожденных живыми детей — мужской или женский. Воспользуемся методом query для выбора строк для младенцев мужского (male) и женского (female) пола:

male = live.query("babysex == 1")

female = live.query("babysex == 2")

len(male), len(female)

(4641, 4500)

Создайте объекты Cdf, представляющие распределение веса при рождении для младенцев мужского и женского пола. Постройте график этих двух ИФР. Каковы различия в форме и положении распределений?

Если ребенок мужского пола весит 8.5 фунта, каков его процентильный ранг? Каков вес ребенка женского пола, имеющего такой же процентильный ранг?

Упражнение 4.3

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

from nsfg import read_fem_preg

 

preg = read_fem_preg()

Упражнение 4.4

Вот данные о скорости бега людей, которые финишировали в забеге James Joyce Ramble, упомянутом ранее в этой главе:

speeds = results["MPH"].values

Создайте объект Cdf, представляющий распределение этих скоростей, и используйте его для вычисления медианы, МКР и квартильной асимметрии. Является распределение скошенным влево или вправо?

Упражнение 4.5

Предполагается, что числа, сгенерированные при помощи функции np.random.random, равномерно распределены в диапазоне от 0 до 1, это означает, что ИФР выборки должна быть прямой линией. Посмотрим, так ли это. Вот выборка из 1001 числа. Постройте для нее график ИФР. Представляет ли он прямую линию?

t = np.random.random(1001)

---