Модуль random в Numpy | Программирование на Питоне

Модуль random в Numpy

Все курсы > Программирование на Питоне > Занятие 11

Библиотека Numpy предоставляет широкие возможности по работе со случайными числами. Для этого в ней предусмотрен модуль random.

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

При этом прежде чем перейти к изучению кода, давайте познакомимся с некоторыми ключевыми понятиями теории вероятностей.

модуль random
Содержание занятия

Случайная величина

Случайная величина (random variable) — это величина, которая в зависимости от случая может принимать определенное значение с заданной вероятностью. Разберем это определение более подробно и параллельно введем несколько новых терминов.

Предположим, что у нас есть некоторое событие (event), например, выпадение двойки или тройки на игральной кости. Это событие случайно (random process), если его невозможно предугадать. Тогда испытанием (trial) мы назовем бросание игральной кости, а исходом (outcome) испытания — выпавшее число.

Вероятность случайной величины

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

Классическая (теоретическая) вероятность

Пусть k — количество благоприятствующих событию A исходов (то есть исходов, которые влекут за собой наступление события A), а n — общее количество равновозможных исходов, тогда теоретической вероятностью события A будет

$$ P(A) = \frac{k}{n} $$

Разберем это определение на том же примере бросания игральной кости.

Выпадение двойки или тройки благоприятствует (приводит к наступлению) случайного события (A). Таких исходов два (либо двойка, либо тройка) и соответственно k = 2. Всего исходов шесть и вероятность их наступления одинакова (они равновозможны). Значит n = 6. Рассчитаем теоретическую вероятность такого события.

$$ P(A) = \frac{2}{6} = \frac{1}{3} $$

Аналогично, вероятность выпадения орла или решки на «честной» монете (fair coin) равна 1/2, а вероятность вытащить червы из колоды, в которой 36 карт, равна 9/36 или 1/4.

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

Сумма вероятностей событий, образующих полную группу, всегда равна единице. Выпадение орла и выпадение решки образуют полную группу (других событий при подбрасывании монеты быть не может). Вероятность каждого из этих событий равна 1/2, сумма этих вероятностей — единице.

Эмпирическая вероятность

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

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

$$ P(A) = \frac{4}{10} = \frac{2}{5} $$

Испытаний (бросков) было проведено очень мало и эмпирическая вероятность довольно сильно отличается от теоретической. Что же будет, если кратно увеличить количество бросков?

Закон больших чисел

Закон больших чисел (Law of large numbers) утверждает, что при проведении множества испытаний эмпирическая вероятность начнет приближаться к теоретической. Другими словами, если бросать кость достаточное количество раз, то эмпирическая вероятность выпадения двойки или тройки приблизится к 1/3.

Откроем ноутбук к этому занятию

Модуль random библиотеки Numpy

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

Подбрасывание монеты

Подбрасывание монеты можно имитировать, обозначив орел и решку через 0 и 1 и случайно генерируя соответствующие числа. Для этого подойдет функция np.random.randint().

На входе функция принимает границы диапазона (верхняя граница в диапазон не входит), на выходе — генерирует целые псевдослучайные числа в его пределах.

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

Бросание игральной кости

Бросание игральной кости можно имиритировать, изменив диапазон на [1, 7).

Метод Монте-Карло

Теперь давайте вернемся к Закону больших чисел. Для его проверки воспользуемся методом Монте-Карло.

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

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

Будем создавать нашу программу поэтапно.

Шаг 1. Создадим список серий испытаний

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

Создать такой список, вернее массив Numpy, можно с помощью функции np.logspace(). Она похожа на уже изученную нами функцию np.linspace() с той лишь разницей, что возвращает заданное количество чисел, которые равномерно распределены по логарифмической шкале в пределах заданного диапазона.

Пусть наша последовательность будет состоять из чисел 10, возведенных в степень от 1 до 6. То есть первым числом будет 101 = 10, а последним 106 = 1 000 000. Всего таких чисел будет 50 (значение функции np.logspace() по умолчанию).

Промежуточные числа Numpy рассчитает самостоятельно с помощью десятичного логарифма. Например, для получения числа 790604 мы возвели 10 в степень 5.897959007291967.

Убедимся в верности вычислений.

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

экспоненциально возрастающая последовательность

Ровно то, что нам нужно.

Шаг 2. Напишем код для проведения испытаний

Используем два цикла:

  • С помощью первого (внешнего) цикла for мы можем имитировать серии бросков (всего 50 испытаний).
  • Второй (внутренний) цикл отвечает за броски каждой отдельной серии. Первый раз мы бросим кость 10 раз, во второй — 12 и так далее.

Пока ограничимся одной серией. Для этого в конце внешнего цикла поставим оператор break.

Шаг 3. Рассчитаем долю успешных испытаний

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

Обратите внимание, эмпирическая вероятность первой серии, в которой было только 10 бросков (1/2), далека от теоретической (1/3).

Шаг 4. Проведем полноценные испытания алгоритма

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

Помимо этого будем замерять время исполнения кода с помощью магической команды %%time.

Общее время выполнения кода составило 26 секунд.

Шаг 5. Оценим результат визуально

метод Монте-Карло

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

Мы убедились в верности Закона больших чисел с помощью метода Монте-Карло.

Пример векторизованного кода

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

Векторизованный код исполнился за 87,1 миллисекунды, то есть почти в 321 раз (26/0,081) быстрее, чем без использования векторизации.

Убедимся, что результат идентичен.

Задача о двух конвертах

Задача о двух конвертах (two envelopes problem) — популярная задача из области теории вероятностей. Сегодня мы рассмотрим один из ее вариантов и смоделируем решение с помощью модуля random.

Представьте, что вам дали конверт A, вы открыли его и увидели сумму X. После этого вам предлагают заменить этот конверт на другой конверт B, сумма в котором с одинаковой вероятностью либо в два раза меньше, либо в два раза больше суммы X. Стоит ли вам заменить конверт с A на B?

задача о двух конвертах

Аналитическое решение

Решим задачу аналитически. Пусть суммы в конверте B могут быть равны X/2 и 2X, тогда в среднем результат замены конвертов составит

$$ \frac{1}{2} \left( \frac{X}{2} \right) + \frac{1}{2} \left( 2X \right) = \frac{5}{4}X $$

Это больше, чем изначальная сумма X в конверте A, а значит нам стоит взять конверт B.

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

$$ \frac{1}{2} \left( \frac{5000}{2} \right) + \frac{1}{2} \left( 2 \cdot 5000 \right) = \frac{5}{4} \cdot 5000 = 6250 $$

Моделирование с помощью модуля random

Проверим это решение с помощью Питона и Закона больших чисел.

Как мы видим, компьютер подтвердил, что в среднем (при достаточном количестве испытаний) будет выгоднее заменить конверт A на конверт B.

Парадокс двух конвертов

Добавлю, что если изменить условие задачи и не открывать конверт A, то задача превратится в парадокс двух конвертов (two envelopes paradox). Ведь если в среднем выгоднее заменить конверт A на B, то будет выгодно и обратное действие (заменить B на A).

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

Рассмотрение этого парадокса выходит за рамки сегодняшнего занятия.

Метод случайных блужданий

Генератор случайных чисел также применяется в так называемом методе случайных блужданий (random walk models, random walkers). С его помощью мы можем конструировать путь, состоящий из последовательности случайных шагов в определенном пространстве и, таким образом, моделировать различные процессы.

Случайное блуждание по числовой прямой

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

Выведем сделанные шаги.

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

Для наглядности давайте посмотрим на сделанные шаги на числовой прямой.

случайное блуждание (random walk) по числовой прямой

Начальное положение (s0) находится в точке ноль. При первом подбрасывании монеты выпадает 0, и на первом шаге (s1) мы остаемся на месте. На втором шаге (s2) выпадает единица, и мы сдвигаемся в точку один. На третьем шаге (s3) мы перемещаемся в точку два. После этого, на шагах четыре (s4) и пять (s5) мы остаемся на месте.

Доедем ли мы из Москвы до Санкт-Петербурга

Симуляция пути из Москвы в Санкт-Петербург

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

  • Нам нужно преодолеть ровно 700 километров
  • Если при бросании кости выпадет единица или двойка, мы сместимся на 5 километров вперед
  • Если тройка или четверка, то на 5 километров назад
  • В случае если выпадет пять или шесть, мы повторно бросим кость и сдвинемся на выпавшее число километров умноженное на пять
  • Кроме того, есть вероятность внезапной поломки. Она составляет менее 0,001, то есть менее 0,1%. В этом случае на эвакуаторе нас отвезут обратно в Москву для ремонта

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

Теперь давайте с помощью Питона построим такой механизм.

  • Как уже было сказано, для бросания игральной кости мы будем использовать функцию np.random.randint() с параметрами 1 и 7.
  • Для того чтобы не уйти в минус, используем функцию max(). Первым параметром передадим 0, вторым — километр, на котором мы окажемся, если сдвинемся назад (то есть если выпадет твойка или четверка). В случае если второе число окажется отрицательным, то функция max() выберет ноль.
  • Поломку мы будем имитировать с помощью функции np.random.rand(). Эта функция генерирует равномерное распределение (об этом ниже) чисел в интервале [0, 1) и, если выпавшее число меньше, чем 0,001, мы вернемся в начало пути.

Посмотрим на код.

Функция случайного блуждания
Создание одного случайного блуждания

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

Выведем результат с помощью графика.

график одного случайного блуждания

Посмотрим на расстояние, которое удалось преодолеть.

Очевидно в этот раз мы не смогли доехать до Санкт-Петербурга.

Моделирование нескольких случайных блужданий

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

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

Для наглядности мы можем использовать DataFrame библиотеки Pandas.

датафрейм из 1000 случайных блужданий

Первый столбец соответствует первому случайному блужданию и на шаге 100 (с индексом 99) мы достигли 430-ого километра.

Выведем результат на графике.

визуализация нескольких случайных блужданий

Обратите внимание на «упавшие» кривые. Это как раз те случаи, когда наше воображаемое транспортное средство ломалось, и мы возвращались в Москву на ремонт.

Вероятность преодоления пути

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

гистограмма распределений случайных блужданий

Выясним, какое расстояние в среднем проезжало наше транспортное средство.

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

Транспортное средство не слишком надежно, лишь в 12,8% случаев мы доезжали до конечной точки.

Продолжим знакомиться с теорией вероятностей и параллельно изучать возможности модуля random.

Дискретные и непрерывные случайные величины

Случайные величины делятся на дискретные и непрерывные.

Дискретная случайная величина

Дискретная случайная величина — это случайная величина, значения которой конечны и счетны.

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

На числовой прямой такая величина может быть отмечана только на определенных точках.

дискретная случайная величина

Непрерывная случайная величина

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

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

непрерывная случайная величина

Вероятностное распределение и модуль random

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

Дискретное вероятностное распределение

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

С помощью функции np.unique() посчитаем, сколько раз выпадет каждый из шести возможных исходов.

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

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

равномерное дискретное распределение

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

Равномерное дискретное распределение

Как вы видите, вероятность наступления каждого исхода примерно одинакова. Такое распределение называется равномерным (discrete uniform distribution). Именно его и использует функция np.random.randint() при генерации чисел.

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

Посмотрим на нотацию равномерного распределения

$$ X \sim U(a, b) $$

В данном случае a и b — минимальное и максимальное значения случайной величины X. При бросании кости a и b равны одному и шести соответственно.

$$ X \sim U(1, 6) $$

Функция вероятности

Помимо графика, дискретное вероятностное распределение удобно описывать с помощью функции вероятности (probability mass function, pmf). Такая функция возвращает вероятность того, что случайная величина примет определенное значение.

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

$$ pmf = \frac {1}{n} $$

Параметр n можно рассчитать по формуле

$$ n = b-a+1 $$

Вероятность каждого исхода, таким образом, будет равна 1/6.

$$ pmf = \frac{1}{6-1+1} = \frac{1}{6} $$

Рассмотрим эту функцию на графике.

функция вероятности равномерного дискретного распределения

Функцию вероятности рассматриваемой нами дискретной величины также можно представить в форме таблицы.

таблица распределения дискретной случайной величины

Еще раз обращу ваше внимание на то, что сумма вероятностей любого распределения равна единице.

Функция распределения

Одновременно, вероятностное распределение бывает удобно описать с помощью (кумулятивной, накопительной) функции распределения (cumulative distribution function, cdf). Эта функция возвращает вероятность того, что случайная величина примет значение меньше заданного. Говоря неформально, она показывает «накопившуюся» вероятность исходов.

Для равномерного дискретного распределения функция распределения выглядит следующим образом

$$ cdf(k; a, b) = \frac{[k]-a+1}{b-a+1} $$

где k — это возможный исход в границах между а и b, то есть k ∈ [a, b].

Приведем пример. При бросании кости рассчитаем вероятность выпадения тройки или меньшего значения (то есть 1, 2 или 3).

$$ cdf(3; 1, 6) = \frac{3-1+1}{6-1+1} = \frac{3}{6} = \frac{1}{2} $$

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

функция распределения равномерного дискретного распределения
Математическое ожидание

Кроме того, вероятностное распределение характеризуется математическим ожиданием (expected value). Матожидание — это среднее значение (mean value) случайной величины X, взвешенное по вероятности каждого из возможных значений.

$$ {\mathbb E}[X] = \sum_{i=1}^{\infty} x_i p_i $$

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

$$ {\mathbb E}[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5 $$

Для равномерного распределения есть и сокращенная формула.

$$ {\mathbb E}[X] = \frac{a+b}{2} = \frac{1+6}{2} = 3.5 $$

На Питоне мы можем сложить результаты поэлементного умножения.

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

Посмотрим на среднее значение сгенерированного распределения.

Дисперсия

Еще одной характеристикой вероятностного распределения является дисперсия (variance). Для равномерного распределение она вычисляется по следующей формуле.

$$ {\mathbb D}[X] = \frac{n^2-1}{12}$$

В нашем случае дисперсия равна

$$ {\mathbb D}[X] = \frac{6^{2}-1}{12} = \frac{35}{12} \approx 2,917 $$

Сравним с нашим распределением

Впрочем, далеко не все дискретные распределения равномерны. Посмотрим на распределение Бернулли.

Распределение Бернулли

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

Давайте рассмотрим этот процесс на практике. На этот раз мы будем подбрасывать монету, но не обычную или «честную монету», а такую, в которой вероятность выпадения орла (обозначим ее p) равна 0,7, а вероятность решки — 0,3. После каждого подбрасывания запишем получившийся результат.

В Numpy нет отдельной функции для распределения Бернулли, поэтому напишем собственную функцию.

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

  • Объявим функцию bernoulli() с двумя параметрами, p и iter. Первый параметр будет отвечать за выпадение орла, второй — за количество подбрасываний.
  • С помощью функции np.random.rand() мы будем генерировать значение непрерывной равномерной величины (об этом опять же чуть ниже) и если получившееся значение меньше p, мы запишем, что выпал орел, в противном случае, что решка.

Посмотрим на реализацию на Питоне.

Вызовем эту функцию с параметрами p и iter, равными 0,7 и 10000 соответственно.

Посмотрим на распределение Бернулли на графике. Воспользуемся столбчатой диаграммой.

распределение Бернулли

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

Пример с одинаковой вероятностью

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

Обозначим орла через H (head), а решку через T (tail).

дерево вероятностей, подбрасывание монеты

Вероятность каждой из комбинаций равна 1/8, так как при одном подбрасывании вероятность выпадения орла или решки одинакова (монета симметрична, fair coin).

$$ P = \frac{1}{2} \times \frac{1}{2} \times \frac{1}{2} = \frac{1}{8} $$

Какова вероятность выпадения двух орлов (вне зависимости от порядка, в котором они выпадали)? Вначале нам нужно посмотреть, в скольких комбинациях оказалось два орла (назовем такие испытания «успехами», successes).

дерево вероятностей, подбрасывание симметричной монеты

Таких комбинаций три (HHT, HTH, THH). Осталось количество «успехов» умножить на вероятность каждого исхода.

$$ P(X = 2) = 3 \times \frac{1}{8} = \frac{3}{8} $$

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

$$ P(X = 0) = \frac{1}{8} $$

Вероятность одного и трех орлов после трех испытаний равна

$$ P(X = 1) = 3 \times \frac{1}{8} = \frac{3}{8} $$

$$ P(X = 3) = \frac{1}{8} $$

Для полноты картины замечу, что, например, P(X = 3) читается как вероятность (P) того, что случайная величина (X) примет значение три.

Теперь давайте выведем эти вероятности на график. По оси x разместим возможные значения случайной величины X, а по оси y соответствующие вероятности P(X = x). У нас получился график биномиального распределения.

график симметричного биномиального распределения

Таким образом, биномиальное распределение (binomial distribution) показывает вероятность количества успехов каждой из возможных комбинаций исходов в серии одинаковых независимых испытаний Бернулли.

Пример с разной вероятностью

Давайте немного изменим условия эксперимента. На этот раз монета будет неправильной или несимметричной (biased coin) и вероятность выпадения орла составит 0,7, а решки — 0,3. Построим дерево вероятностей (tree diagram).

дерево вероятностей, подбрасывание несимметричной монеты

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

$$ P(X = 1) = 3 \times 0,063 = 0,189 $$

$$ P(X = 2) = 3 \times 0,063 = 0,441 $$

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

$$ P(X = 0) = 0,027 $$

$$ P(X = 3) = 0,343 $$

Построим такое биномиальное распределение на графике.

график несимметричного биномиального распределения

Обратите внимание, распределение несимметрично с более пологой частью слева. Еще говорят, что распределение скошено влево (skewed left). Если бы вероятность выпадения орла была 0,3, а решки — 0,7, то распределение имело бы обратную форму и было скошено вправо (skewed right).

Формула биномиального распределения

Пока мы бросали монету по три раза, то вероятность любой комбинации исходов несложно посчитать с помощью дерева вероятностей. При этом если подбрасывать, например, 10 раз, количество вариантов достигнет 210 = 1024. Для такого распределения понадобится формула.

Возможно вы обратили внимание, что для выведения формулы нам нужно два компонента:

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

Количество каждой из возможных комбинаций можно найти с помощью треугольника Паскаля. Например, при трех подбрасываниях (n = 3), мы видим, что количество возможных комбинаций будет равно {1, 3, 3, 1}, что соответствует нашему дереву вероятностей.

треугольник Паскаля для расчета количества комбинаций

Можно также воспользоваться формулой биномиальных коэффициентов, в которой через n мы обозначим количество подбрасываний, а через k — количество «успехов» (выпадений орлов).

$$ \binom{n}{k} = \frac{n!}{k!(n-k)!} $$

Например, убедимся, что выпадение двух орлов (k = 2) при трех подбрасываниях (n = 3) возможно в трех комбинациях.

$$ \binom{3}{2} = \frac{3!}{2!(3-2)!} = 3 $$

Одновременно, нам нужно понимать вероятность каждой комбинации. Ее можно рассчитать по формуле

$$ p^k (1-p)^{n-k} $$

В данном случае мы возводим вероятность успеха (p) в степень количества успехов (k) и вероятность неудачи (1−p) в степень количества неудач (n−k). Например, если вероятность выпадения орла (p) равна 0,7, то вероятность выпадения двух орлов (k = 2) в трех бросках (n = 3) равна

$$ 0,7^2 (1-0,7)^{3-2} = 0,49 \times 0,3 = 0,147 $$

Остается перемножить количество комбинаций и вероятность каждой из них

$$ P(X = 2) = 3 \times 0,147 = 0,441 $$

Таким образом, вся формула целиком выглядит так

$$ P(X = k) = \binom{n}{k} p^k (1-p)^{n-k} $$

Эта же формула описывает функцию вероятности (probability mass function, pmf) биномиального распределения.

Матожидание и дисперсия

Матожидание рассчитывается по формуле

$$ {\mathbb E}[X] = np $$

Например, после трех бросков мы можем ожидать, что в среднем орел выпадет 3 x 0,7 = 2,1 раза. Дисперсию можно рассчитать через

$$ {\mathbb D}[X] = np(1-p) \rightarrow {\mathbb D}[X] = 3 \times 0,7 \times (1-0,7) = 0,63 $$

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

Выполним эти же расчеты с помощью Питона. Функция np.random.binomial() принимает три параметра:

  • n — количество испытаний в одном эксперименте (например, подбрасываний монеты)
  • p — вероятность успеха (например, выпадения орла)
  • size — количество экспериментов (серий по n подбрасываний)

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

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

Посмотрим на распределение на графике.

создание биномиального распределения через модуль random

Рассчитаем среднее значение и дисперсию.

Дополнительно, как и обещал, покажу как можно повторить эксперимент Бернулли с помощью функции np.random.binomial(). Для этого достаточно указать параметр n = 1 (то есть только одно подбрасывание).

Выведем распределение на графике.

создание распределения Бернулли через модуль random

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

Непрерывное вероятностное распределение

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

Равномерное непрерывное распределение

Непрерывное равномерное распределение (continuous uniform distribution) описывает случайную величину, вероятность значений которой одинакова на заданном интервале от a до b.

$$ X \sim U(a, b) $$

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

$$ X \sim U(0, 12) $$

Плотность вероятности

Непрерывное распределение (в отличие от дискретного) задается плотностью вероятности (probability density function, pdf). Для равномерного непрерывного распределения плотность вероятности задается вот такой несложной функцией.

$$ pdf(x) = \begin{cases} \frac{1}{b-a}, x \in [a, b] \\ 0, x \notin [a, b] \end{cases} $$

В примере с ожиданием автобуса вероятность его приезда в любой момент в пределах заданного интервала равна

$$ pdf(x) = \begin{cases} \frac{1}{12-0} = \frac{1}{12}, x \in [0, 12] \\ 0, x \notin [0, 12] \end{cases} $$

На графике равномерное распределение представляет собой прямоугольник, площадь которого всегда равна единице.

график плотности вероятности равномерного непрерывного распределения

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

Например, вероятность приезда автобуса при ожидании до 12 минут включительно составляет 1.00 или 100%, потому что такой промежуток включает всю площадь прямоугольника.

Теперь давайте рассчитаем вероятность ожидания автобуса до 7 минут включительно. Нас будет интересовать интервал от 0 до 7 минут и соответствующий участок площади прямоугольника.

расчет вероятности по графику равномерного непрерывного распределения

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

$$ P(7) = \frac{1}{12} \times 7 \approx 0,583 $$

Матожидание и дисперсия

Остается рассчитать матожидание (среднее время ожидания автобуса) и дисперсию.

$$ {\mathbb E}[X] = \frac{a + b}{2} = \frac{0 + 12}{2} = 6 $$

$$ {\mathbb D}[X] = \frac{(b-a)^2}{12} = \frac{(12-0)^2}{12} = 12 $$

Реализация на Питоне

Воспользуемся функцией np.random.uniform() для того, чтобы создать равномерное распределение с параметрами U(0, 12).

Посмотрим на результат с помощью гистограммы.

создание равномерного непрерывного распределения с помощью np.random.uniform()

Посмотрим на среднее значение и дисперсию.

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

Разница между np.random.random(), np.random.rand() и np.random.uniform()

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

Функция np.random.random(size = None) создает равномерное распределение в полуоткрытом интервале [0, 1). Параметр size задает размер этого распределения (количество экспериментов).

Функция np.random.rand() практически идентична.

Для функции np.random.uniform(low = 0.0, high = 1.0, size = None) интервал [0, 1) является интервалом по умолчанию, при этом можно задать любой другой промежуток (как мы и сделали выше).

Приведем несколько примеров.

Напоследок замечу, что равномерное непрерывное распределение является частным случаем бета-распределения с параметрами Beta(1, 1). О том что это такое мы поговорим с вами на курсе статистики вывода.

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

На занятии по описательной статистике мы начали изучать количественные данные с примера роста мужчин в России. При этом интересно, что рост людей, как и многие другие величины (например, вес человека, артериальное давление, некоторые природные явления и многое другое) имеют так называемое нормальное распределение (normal distribution).

Функция плотности нормального распределения

Функция плотности (pdf) нормального распределения случайной величины X определяется функцией Гаусса и поэтому нормальное распределение также называется распределением Гаусса (Gauss distribution).

$$ pdf(x, \mu, \sigma) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2} (\frac{x-\mu}{\sigma})^2} $$

Как видно из формулы, единственными неизвестными параметрами являются μ («мю», матожидание) и σ («сигма», среднее квадратическое отклонение, СКО). Именно они и определяют нормальное распределение.

$$ X \sim {\mathcal N}(\mu, \sigma) $$

Здесь важно напомнить, что среднее квадратическое отклонение (standard deviation) равно квадратному корню из дисперсии (variance).

$$ \sigma = \sqrt{\sigma^2} $$

Функция np.random.normal()

На Питоне нормальное распределение создается с помощью функции np.random.normal(). Мы уже использовали ее, в частности, для создания данных о росте мужчин и женщин в ноутбуке⧉ к десятому занятию вводного курса. Повторим этот код, увеличив размер массива до 100 000.

Посмотрим на результат.

Для визуализации нормального распределения можно использовать несколько типов графиков. Например, можно использовать уже знакомую нам гистограмму.

две гистрограммы нормального распределения на одном графике

Также можно использовать график плотности (density plot) распределения случайной величины.

два графика плотности распределения

Еще одной полезной визуализацией является так называемый boxplot (ящик с усами или диаграмма размаха).

датафрейм с данными о росте
два boxplot (диаграммы размаха) на одном графике

Boxplot позволяет увидеть медиану (median), первый и третий квартили (Quartile 1, Q1 и Quartile 3, Q3), межквартильный размах (Interquartile Range, IQR), а также, что очень важно, так называемые выбросы (outliers), то есть значения, сильно отличающиеся от среднего (на графике выше они обозначены точками). Ни гистограмма, ни график плотности вероятности этой информации не выводят.

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

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

С σ («сигма»), обозначающей среднее квадратическое отклонение мы уже знакомы, а вот что такое квартиль, межквартильный размах и о том, что измеряют Q1 − 1.5 x IQR и Q3 + 1.5 x IQR мы поговорим на следующем курсе по анализу и обработке данных.

Дополнительно приведу пример того, как можно совместить boxplot с гистограммой на Питоне.

boxplot и гистрограмма на одном графике
Расчет вероятности

Теперь давайте рассчитаем вероятность того, что рост случайно встретившегося нам человека на улице составляет менее 190 см.

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

Вначале рассчитаем теоретическую вероятность. Для создания «идеального» теоретического распределения воспользуемся библиотекой scipy.

кривая нормального распределения

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

площадь под кривой нормального распределения до отметки 190 см

Для нахождения площади можно воспользоваться одним из многочисленных онлайн-калькуляторов⧉.

калькулятор площади под кривой нормального распределения

На Питоне это вычисление можно выполнить с помощью функции распределения (cumulative density function, cdf).

Обратное вычисление, то есть нахождение значения (роста) по площади выполняется с помощью квантиль-функции (percent point function, ppf).

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

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

площадь под кривой нормального распределения, начиная с отметки 190 см

Если нужно вычислить площадь между двумя значениями a и b, например, между 170 и 190 см, из большего, находящегося правее значения функции распределения, можно вычесть из меньшее, находящееся левее.

Посмотрим на этот участок на графике.

площадь под кривой нормального распределения между отметками 170 см и 190 см

Теперь вернемся к нашим данным и рассчитаем эмпирическую вероятность встретить человека с ростом не более 190 см.

В целом результат близок к найденному аналитическому решению.

Функция плотности и функция распределения

Рассмотрим связь функции плотности (pdf) и функции распределения (cdf). Напомню, что функция плотности нормального распределения определяется по формуле

$$ P(x) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2} (\frac{x-\mu}{\sigma})^2} $$

При этом вот что мы успели про нее узнать:

  1. Вероятность того, что случайная величина примет значение не более заданного в интервале (−∞; x] равна площади под кривой функции плотности на этом промежутке
  2. Эта площадь вычисляется с помощью функции распределения

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

Значит функция распределения (cdf) есть интеграл функции плотности (pdf).

Математически это выражается так.

$$ D(x) = \int_{-\infty}^{x} P(x) dx $$

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

связь функции плотности (pdf) и функции распределения (cdf) через интегрирование

Вначале обратимся к графику слева. Как мы видим, вероятность встретить человека не более 180 см составляет 0,5 (закрашенный синим участок). Одновременно, проинтегрировав функцию плотности на интервале (−∞; 180], на графике справа мы видим, что «накопленная» вероятность составляет 0,5, и именно эту вероятность нам показывает график функции распределения на отметке x = 180.

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

Если функция распределения есть интеграл функции плотности, то плотность вероятности (pdf) является производной функции распределения (cdf).

Приведем формулу.

$$ P(x) = D'(x) $$

Хотя опять же тему дифференцирования (то есть нахождения производной) мы будем разбирать на курсе по оптимизации, ничто не мешает нам начать знакомиться с ней прямо сейчас.

Совместим обе функции на одном графике.

связь функции распределения (cdf) и функции плотности (pdf) через производную

Здесь вначале посмотрим на функцию распределения (cdf, оранжевый график). В промежутке от 150 до 210 см эта функция непрерывно возрастает, при этом она возрастает с разной скоростью. До точки x = 180 (она называется точкой перегиба, inflection point) скорость возрастания функции увеличивается, после нее убывает.

Именно это изменение описывает производная от нее функция плотности (pdf, синий график). На участке от 150 до 180 она возрастает, а потом в интервале от 180 до 210 постоянно убывает.

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

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

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

Вероятность конкретного значения

Важный и немного контринтуитивный момент. Вероятность того, что случайная величина непрерывного распределения примет конкретное значение (то есть нам встретится человек определенного роста) равна нулю. Это легко продемонстрировать. Ранее мы находили вероятность встретить человека ростом от 170 до 190 сантиметров, вычитая меньшую площадь под кривой из большей.

$$ P(170 \leq x \leq 190) = P(x \leq 190)-P(x \leq 170) \approx 0,68 $$

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

$$ P(x = 190) = P(x \leq 190)-P(x \leq 190) = 0 $$

Формирование выборки

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

модуль random: формирование выборки

Пусть таким набором элементов будет мешок с разноцветными шарами.

При этом формировать выборку можно двумя способами. В первом случае мы случайным образом достаем элемент, но, прежде чем взять следующий, кладем этот элемент обратно. Такой процесс называют выборкой с возвращением (sampling with replacement), и имитировать его можно с помощью функции np.random.choice().

Обратите внимание, белый шар повторяется дважды.

Кроме того, можно сформировать выборку без возвращения (sampling without replacement). В этом случае мы не кладем элемент обратно, а откладываем в сторону и только потом достаем следующий элемент. Для этого функции np.random.choice() нужно задать параметр replace = False.

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

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

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

Определения и нотация

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

Во-первых, напомню, что данные могут представлять собой генеральную совокупность (population) или выборку из нее (sample). Если мы берем несколько выборок из одной генеральной совокупности и измеряем для каждой из них определенный параметр (например, среднее арифметическое, sample mean), то совокупность этих средних формирует выборочное распределение (sampling distribution).

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

центральная предельная теорема: среднее средних нескольких выборок

Теперь перейдем к сути.

ЦПД и нормальное распределение

Мы уже знаем, что среднее средних нескольких выборок (mean of sampling distribution of sample means) из одной генеральной совокупности будет стремиться к истинному среднему этой генеральной совокупности.

$$ \mu_\bar{x} = \mu $$

Однако это не все. При соблюдении двух условий, а именно, (1) выборки сформированы случайным образом, (2) размер каждой выборки составляет не менее 30 элементов, выборочные средние будут следовать нормальному распределению.

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

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

$$ \sigma_\bar{x} = \frac{\sigma}{\sqrt{n}} $$

Математически эти выводы можно записать так

$$ X \sim dist(\mu, \sigma) \rightarrow \bar{X} \sim {\mathcal N}(\mu, \frac{\sigma}{\sqrt{n}}) $$

Проверим на Питоне

А теперь давайте проверим истинность ЦПД с помощью Питона. Для начала создадим скошенное вправо распределение (right skewed distribution). Такое распределение может характеризовать, например, зарплаты людей.

Посмотрим на график этого распределения.

скошенное вправо распределение (right skewed distribution)

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

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

Рассчитаем СКО.

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

Посмотрим на результат.

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

Обратите внимание, распределение выборочных средних гораздо ближе к нормальному распределению, нежели исходное распределение salaries. Остается рассчитать значения, к которым, согласно ЦПД, должны стремиться среднее значение и СКО выборочных средних.

Посмотрим так ли это.

Как мы видим, ЦПД выполняется.

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

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

Важно, что при формировании выборки без возвращения ЦПД будет выполняться, если размер одной выборки состаляет не более 5% от размера генеральной совокупности.

Для формирования распределения выборок без возвращения напишем собственную функцию sample_means(). На входе она будет принимать следующие параметры:

  • data — набор данных (генеральную совокупность);
  • n_samples — количество выборок
  • sample_size — размер одной выборки
  • replace = True — с возвращением делать выборки или без
  • random_state = None — воспроизводимость результата

Протестируем эту функцию.

распределение выборочных средних без возвращения

Теперь сгенерируем несколько распределений выборочных средних с 20, 100 и 500 выборками в распределении и размером выборки в 2, 10 и 30 значений.

график: распределения выборочных средних с 20, 100 и 500 выборками в распределении и размером выборки в 2, 10 и 30 значений

Кроме того, давайте посмотрим на параметры этих распределений в табличной форме и сравним с «целевыми» показателями, основанными на ЦПД.

датафрейм: распределения выборочных средних с 20, 100 и 500 выборками в распределении и размером выборки в 2, 10 и 30 значений

Как мы видим, на выполняемость ЦПД влияет не только размер выборки, но и количество этих выборок. Очевидно, что наибольшую близость к расчетным показателям демонстрируют распределения из 500 выборок.

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

Любое нормальное распределение со средним значением μ и СКО σ можно привести к стандартному нормальному распределению (standard normal distribution) со средним значением ноль и СКО равным единице.

$$ Z \sim {\mathcal N}(0, 1) $$

Для этого воспользуемся следующей формулой.

$$ z = \frac{x-\mu}{\sigma} $$

Таким образом мы приводим каждое значение x к соответствующей z-оценке (z-score), вычитая среднее μ и деля результат на СКО σ. Например, приведем данные о росте мужчин к стандартному виду (для этого воспользуемся векторизацией и трансляцией кода).

Посмотрим на результат на графике.

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

Стоит сказать, что ровно такого же результата можно добиться, применив метод .fit_transform() класса ScandardScaler модуля preprocessing библиотеки sklearn.

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

Добавлю, что стандартное нормальное распределение можно также создать с помощью функции np.random.standard_normal().

Посмотрим на результат.

стандартное нормальное распределение через функцию np.random.standard_normal()

Создаваемый массив также может быть многомерным.

Как и в случае обычного нормального распределения, мы можем найти значение (в данном случае z-score) по площади с помощью квантиль-функции.

Убедиться в верности результата можно с помощью функции распределения.

Критерии нормальности распределения

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

Одновременно пока что мы определяли нормальность распределения «на глаз» (с помощью гистограммы, графика плотности или boxplot), хотя в некоторых случаях полезно иметь более надежный критерий нормальности. Это важно, например, при оценке:

  • Распределения остатков модели линейной регрессии
  • Распределения наблюдений в классах модели линейного дискриминантного анализа (Linear Discriminant Analysis, LDA)

Рассмотрим два способа оценки нормальности распределения.

Способ 1. График нормальной вероятности

График нормальной вероятности (normal probability plot) показывает соотношение упорядоченных по возрастанию данных и соответствующих им квантилей нормального распределения.

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

Алгоритм создания графика нормальной вероятности следующий:

  1. Сортируем исходные данные по возрастанию
  2. Находим накопленную вероятность (cumulative probability) каждого значения
  3. С помощью квантиль-функции выясняем, какому квантилю соответствовала бы эта вероятность, если бы распределение было нормальным
  4. По оси x отмечаем квантили, по оси y — отсортированные данные

Накопленную вероятность будем вычислять по формуле

$$ P = \frac{i-0.375}{n+0.25} $$

где i — это индекс (начиная с единицы) значения в перечне отсортированных данных, а n — количество наблюдений.

С помощью Питона построим график нормальной вероятности для данных о росте.

график нормальной вероятности, построенный вручную

Как мы видим, в целом данные распределены нормально (что разумеется ожидаемо, потому что мы генерировали их с помощью функции np.random.normal()).

Можно также воспользоваться функцией probplot() модуля stats библиотеки scipy.

Построим график для тех же данных о росте.

график нормальной вероятности с помощью функции probplot()

Эта функция использует другую формулу для вычисления накопленной вероятности (а именно Filliben’s estimate), поэтому квантили незначительно, но будут отличаться от написанного нами алгоритма.

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

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

Здесь мы видим, что график нормальной вероятности довольно сильно отклоняется от «идеальных» значений, лежащих на диагонали.

Остается построить график для распределения выборочных средних.

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

Как и должно быть, это распределение гораздо ближе к нормальному.

Способ 2. Тест Шапиро-Уилка

Тест Шапиро-Уилка (Shapiro-Wilk test) позволяет сделать статистически значимый вывод о нормальности распределения.

  • Нулевая гипотеза предполагает, что распределение нормально
  • Альтернативная гипотеза утверждает обратное

Тест Шапиро-Уилка чувствителен к количеству элементов (N) в наборе данных и теряет точность при N > 5000.

Проведем тест для распределений роста и распределения зарплат при пороговом значении 0,05, однако вначале создадим распределения с меньшим количеством элементов.

Начнем с роста.

Теперь посмотрим на зарплаты.

Как и ожидалось, в первом случае мы не смогли отвергнуть нулевую гипотезу о нормальном распределении (p-value > 0,05), во втором случае мы можем это сделать (p-value < 0,05).

Нормальное приближение биномиального распределения
Теорема Муавра-Лапласа

Теорема Муавра-Лапласа (de Moivre–Laplace theorem), частный случай Центральной предельной теоремы, утверждает, что при определенных условиях нормальное распределение может быть использовано в качестве приближения биномиального распределения (Normal Approximation to Binomial Distribution).

Напомню формулу биномиального распределения

$$ P(X = k) = \binom{n}{k} p^k (1-p)^{n-k} $$

где n — количество испытаний, k — количество успехов, а p — вероятность успеха. Так вот если np ≥ 5 и n(1−p) ≥ 5, то выполняется следующее

$$ B(n, p) \sim {\mathcal N}(np, \sqrt{np(1-p)}) $$

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

Проиллюстрируем теорему с помощью Питона. Будем подбрасывать несимметричную монету (p = 0,8) по 3, 5, 10, 15, 25 и 50 раз и сравнивать получившиеся распределения с нормальным.

нормальное приближение биномиального распределения для n = 3 и n = 5
нормальное приближение биномиального распределения для n = 10 и n = 15
нормальное приближение биномиального распределения для n = 25 и n = 50

Как мы видим, вначале (например, при n = 3), распределение было ожидаемо скошенным, при этом по мере увеличения количества бросков оно стало все больше «вписываться» в график плотности нормального распределения.

Поправка на непрерывность распределения

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

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

поправка на непрерывность распределения 1

Получаем

$$ P_B(X ≤ 1) = P_B(X = 0) + P_B(X = 1) = \frac{1}{8} + \frac{3}{8} = \frac{4}{8} = \frac{1}{2} $$

При этом, если мы хотим рассчитать эту же площадь с помощью кривой нормального распределения, нам нужно сместить границу таким образом, чтобы захватить все интересующие нас столбцы. В данном случае прибавить 0,5. То есть PN (X < 1,5).

поправка на непрерывность распределения 2

Воспользуемся теоремой Муавра-Лапласа для расчета среднего и СКО (напомню n = 3, p = 0,5).

$$ \mu = np = 3 \times 0,5 = 1,5 $$

$$ \sigma = \sqrt{np(1-p)} = \sqrt{3 \times 0,5 \times 0,5} = 0,75 $$

Остается воспользоваться Питоном для нахождения площади под кривой нормального распределения.

Посмотрим на результат на графике.

поправка на непрерывность распределения: нахождение площади под кривой

Если же нас попросят найти вероятность конкретного значения, например, выпадения двух орлов PB(X = 2) (то есть площадь одного столбца гистограммы), то при расчете площади под кривой нормального распределения нужен интервал PN(1,5 < X < 2,5). Другими словами мы прибавили по 0,5 с обеих сторон.

поправка на непрерывность распределения 3

Рассчитаем площадь с помощью Питона.

Напомню, площадь столбца составляет PB(X = 2) = 3/8 или 0,375, что довольно существенно отличается от площади под кривой. Это связано с тем, что мы взяли слишком маленькое n и условия np ≥ 5 и n(1−p) ≥ 5 в данном случае не выполняются.

По мере увеличения n поправка на непрерывность распределения становится все менее значимой.

Выведем результат на графике.

поправка на непрерывность распределения: нахождение площади под кривой 2
Пример приближения

Для закрепления пройденного материала рассмотрим более жизненный пример. Предположим, вам поставили партию из 500 единиц оборудования. При этом вам известно, что в среднем 2% (то есть 0,02) оборудования имеют различные дефекты. Какова вероятность того, что в партии не менее 15 бракованных единиц оборудования?

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

$$ B(n, p) = B(500; 0,02) \rightarrow P_B(X \geq 15) $$

Теперь решим эту задачу с помощью нормального приближения. По теореме Муавра-Лапласа выводим следующее.

$$ B(500; 0,02) \sim {\mathcal N}(10; \sqrt{9.8}) \rightarrow P_N(X > 14,5) $$

Таким образом, задача сводится к нахождению площади под кривой с учетом внесенной поправки на непрерывность распределения (15 − 0,5 = 14,5).

Посмотрим, как это выглядит на графике.

пример нормального приближения биномиального распределения

Перейдем к комбинаторике.

Комбинаторика в Питоне

np.random.shuffle() и np.random.permutation()

Обе эти функции перемешивают (shuffle) или переставляют (permute) элементы списка или массива. Основное отличие заключается в том, что np.random.shuffle() изменяет исходный массив, а np.random.permutation() создает «перемешанную» копию исходного массива:

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

Подгрузим данные о пациентах с диабетом из модуля datasets библиотеки sklearn. Как это делается, мы подробно разобрали на занятии по линейной регрессии на вводном курсе.

Вызовем созданную функцию и проверим ее работу.

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

модуль random, пример применения функции permutations()
модуль random, второй пример применения функции permutations()

Похоже, что все работает как надо.

Модуль itertools

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

Перестановки

Напомню, что перестановка (permutation) в комбинаторике позволяет узнать, сколькими способами можно разместить элементы множества в определенном порядке.

Перестановки без замены

Вначале возьмем те случаи, когда один элемент может повторяться только один раз. В этом случае мы говорим о перестановках без замены (permutations without replacement). Они бывают двух видов: без повторений и с повторениями.

Перестановки без повторений

Рассмотрим простой пример. У нас есть три элемента A, B и C. Сколькими способами можно переставить эти объекты так, чтобы порядок элементов не повторялся?

Так как элементы в исходном множестве не повторяются, то такая перестановка называется перестановкой без повторения (permutation without repetition).

Очевидно, на первое место мы можем поставить любой из трех элементов, на второе только оставшиеся два, на третьем месте окажется только один элемент. Получается 3 x 2 x 1 = 6. Такой расчет не что иное, как факториал числа три (3!).

$$ P(n) = n! \rightarrow P(3) = 3! = 6 $$

На Питоне мы можем вычислить факториал нужного нам числа через функцию factorial() модуля math.

Кроме того, мы можем воспользоваться функцией permutations() модуля itertools. В этом случае мы можем вывести все возможные перестановки элементов множества.

Если все же нас интересуют не сами перестановки, а их количество, можно воспользоваться функцией len().

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

Приведем формулу.

$$ P(n, r) = \frac{n!}{(n-r)!} \rightarrow P(6, 3) = \frac{6!}{(6-3)!} = 120 $$

На Питоне мы можем применить ту же функцию permutations(), передав ей через параметр r размер выборки (количество мест).

Посмотрим на общее количество перестановок.

Дополнительно замечу, что в случае если выборка r также велика, как и множество n (первый пример), то формулу перестановок можно записать как

$$ P(n, n) = \frac{n!}{(n-n)!} $$

В знаменателе получается 0!, и для того чтобы избежать деления на ноль, принято считать, что факториал нуля равен единице (0! = 1).

Перестановки с повторениями

Давайте дополнительно усложним задачу и рассчитаем количество перестановок букв в слове «молоко». В данном случае мы не можем воспользоваться факториалом, потому что буква «о» повторяется трижды, и соответственно мы можем разместить ее 3! способами. Как следствие, для получения правильного ответа нам нужно разделить количество перестановок на шесть (3! = 6).

В таком случае мы говорим про перестановки с повторениями (permutations with repetitions).

$$ \frac{6!}{3!} = 120 $$

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

$$ \frac{P(n, r)}{n_1!, n_2!, \dots, n_r!}, \text{где } n_1 + n_2! + \dots + n_r! = n $$

Для расчета перестановок с повторениями напишем собственную функцию.

Рассчитаем количество перестановок с повторениями букв в слове «молоко».

Перестановки с заменой

Перестановки с заменой (permutations with replacement), когда элементы могут использоваться более одного раза, есть не что иное, как прямое или декартово произведение множества само на себя.

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

Декартово произведение (Cartesian product) двух множеств предполагает, что мы создадаем все возможные упорядоченные пары исходных множеств.

В частности, если у нас есть одно множество A с элементами {1, 2, 3} и множество B с элементами {x, y}, то их произведением A x B будут все возможные упорядоченные комбинации этих элементов.

произведение двух множеств

При умножении множества само на себя один раз мы говорим про квадрат Декарта (Cartesian square). В частности, если умножить множество действительных чисел R само на себя (R x R), то получится координатная плоскость или декартова система координат (Cartesian plane).

координатная плоскость

Замечу, что множество можно умножить само на себя много раз, в этом случае говорят про декартову степень (Cartesian power).

Так вот, перестановка с заменой и есть декартова степень исходного множества.

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

декартово произведение двух множеств

При этом так как это перестановки и порядок имеет значение, (клубника, ваниль) и (ваниль, клубника) это разные блюда.

В модуле itertools для такой операции есть