• Основы теории вероятностей и распределения вероятностей для случайной величины.
• Биномиальное и нормальное распределения, распределение Пуассона.
• Кумулятивная функция нормального распределения и обратная функция.
• Аппроксимация одного распределения другим в определенных случаях

Биномиальный эксперимент

Кратко:

• Биномиальный эксперимент - это эксперимент, в котором подряд проводят несколько одинаковых и независимых друг от друга испытаний Бернулли.
• Примеры биномиальных экспериментов: пользователи кликают по баннеру, завод производит детали, питон в питон питомнике вырастает добрым или злым.
• Расчёт вероятности определённого порядка успехов и неуспехов.
• Вероятность того, что в биномиальном эксперименте успехи и неудачи произойдут в определённом порядке.
• Вероятности исходов биномиального эксперимента с двумя испытаниями по схеме Бернулли.
• Распределение вероятностей биномиального эксперимента с двумя испытаниями по схеме Бернулли

p = 0.9 # вероятность, что питон добрый
prob = p * p * (1 - p)
print(prob)

Вы узнали, как рассчитывать вероятность исхода биномиального эксперимента, в котором задан определённый порядок успехов и неудач. Теперь посмотрим, как рассчитывать вероятности в случае, если порядок неважен. Рассмотрим на простом примере — испытаний по схеме Бернулли всего 2.

Вероятности исходов биномиального эксперимента с двумя испытаниями по схеме Бернулли

Посмотрим на питонов в питонопитомнике. Вероятности оставим те же: питон вырастает добрым с вероятностью 90 %, питон переходит на тёмную сторону с вероятностью 10 %.

 

Нам подходят два исхода: «добрый-злой» и «злой-добрый» — и не подходят тоже два исхода: «оба добрые» и «оба злые». Вычислим вероятность неподходящих исходов. Вероятность исхода «оба добрые» равна 0.92=0.81, «оба злые» равна 0.12=0.01. Вероятность всех неподходящих исходов равна 81 %+1 %=82 %. На подходящие комбинации остаётся 18 %. Тот же результат, что и в задании!

 

Если бы задача звучала так: «Какова вероятность, что первый питон злой, а второй — добрый?» — то решением был бы всего один исход. Учитывать комбинации «злой-добрый» и «добрый-злой» было бы не нужно: достаточно умножить 0.1 на 0.9, чтобы получить ответ.

 

• Вероятность, что оба питона будут с зелёными глазами, равна 0.25^2=0.0625=6.25 %.
• Вероятность, что цвет глаз у обоих питонов голубой, равна 0.75^2=0.5625=56.25 %.
• Вероятность, что один с зелёными, а другой с голубыми: 0.25⋅0.75⋅2=37.5 %.

Биномиальный эксперимент с тремя испытаниями Бернулли

Комбинацию исходов «два питона добрые, а один злой» можно получить тремя способами:

• «злой, добрый, добрый»;
• «добрый, злой, добрый»;
• «добрый, добрый, злой».

Найдём вероятность этой комбинации. Вероятности для каждого питона вырасти добрым или злым оставим прежними: 0.9 и 0.1 соответственно.

 

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

• P(злой,добрый,добрый)=0.1⋅0.9⋅0.9=0.081
• P(добрый,злой,добрый)=0.9⋅0.1⋅0.9=0.081
• P(добрый,добрый,злой)=0.9⋅0.9⋅0.1=0.081

Получается, в событие «два питона добрые, а один злой» входят три исхода, вероятность каждого из которых равна 8.1 %.

 

Чтобы получить вероятность всего события «два питона добрые, а один злой», нужно сложить вероятности исходов, входящих в него. Вероятности исходов равны, так что просто умножим одну на их количество: 0.081⋅3=0.243=24.3 %

 

Значит, вероятность, что из трёх питонов два добрые, а один злой — чуть меньше 25 %.

 

Посмотрим ещё раз на расчёт:

• Эксперименты независимы, поэтому для конкретного числа добрых питонов из трёх вероятности каждой комбинации одинаковы: мы перемножаем одни и те же числа, просто в разном порядке.
• Три разных пути получения двух добрых питонов из трёх на иллюстрации выше — это три способа выбрать одного злого питона из трёх. Это же количество сочетаний из трёх по одному! Вспомните формулу сочетаний из темы «Комбинаторика»:
• 
  Cnk=n!/k!(n−k)!
  C31 как раз и равно трём.
  Кстати, можно выбирать не одного злого питона из трёх, а двух добрых из трёх — как можно предположить, комбинаций будет столько же. В таком случае получится C32​. Оно тоже равно трём

Мы рассмотрели пример с тремя подряд экспериментами Бернулли. Обобщим изученное.

Вероятность произвольного числа успехов из произвольного числа повторений испытания Бернулли

Пусть мы провели n экспериментов по схеме Бернулли, вероятность успеха в которой равна $.$ Нас интересует вероятность получить k успехов $.$ Действуем по той же схеме, что и для двух и трёх экспериментов. Получаем такой алгоритм:

 

Чтобы найти вероятность получить k успехов из n экспериментов по схеме Бернулли, надо:

1. Перемножить вероятности всех успехов и неудач. Каждый исход состоит из $k$ успехов и n−k неудач. Вероятность успеха равна p, всего их k штук, значит, вероятность всех успехов равна pk. Вероятность неудачи равна $,$ всего их n−k штук, значит, вероятность всех неудач (1−p)n−k. Получаем вероятность исхода с $k$ успехами из $n$ испытаний: pk⋅(1−p)n−k.
2. Вычислить число исходов, которые приведут к нужному событию. Это число сочетаний из n по $.$
3. Перемножить количества из п. 1 и п. 2.

Запишем в виде общей формулы:

 

$.$

Графики биномиальных распределений

from matplotlib import pyplot as plt
from math import factorial
n = 5 # количество попыток
p = 0.5 # вероятность успеха
distr = []
for k in range(0, n + 1):
    choose = factorial(n) / (factorial(k) * factorial(n - k))

from matplotlib import pyplot as plt
from scipy.stats import binom
n = 26 # количество попыток
p = 0.9 # вероятность успеха
distr = []
for k in range(0, n + 1):
    current_value = binom.pmf(k, n, p)
    distr.append(current_value)
plt.bar(range(0, n + 1), distr)

from matplotlib import pyplot as plt
from math import factorial
p = 0.2 # напишите ваш код здесь: какова вероятность заключить контракт?
n = 30 # напишите ваш код здесь: со сколькими компаниями начинаем переговоры?
# напишите ваш код здесь: создайте список distr, в котором будете хранить значения распределения
distr = []
for k in range(0, n + 1):
	   # напишите ваш код здесь: постройте распределение вероятностей
    choose = factorial(n) / (factorial(k) * factorial(n - k))
    prob = choose * p**k * (1 - p)**(n - k) 
    distr.append(prob)
# напишите ваш код здесь: постройте гистограмму распределения вероятностей
plt.bar(range(0, n + 1), distr)

from matplotlib import pyplot as plt
from math import factorial
p = 1/5
n = 60
distr = []
for k in range(0, n + 1):
    choose = factorial(n) / (factorial(k) * factorial(n - k))
    prob = choose * p**k * (1-p)**(n-k) 
    distr.append(prob)
plt.bar(range(0, n + 1), distr)

from matplotlib import pyplot as plt
from math import factorial
n_exams = 6 # напишите ваш код здесь: сколько экзаменов надо сдать?
failure_rate = 0.15 # напишите ваш код здесь: какова вероятность завалить один экзамен?
distr = [] # создайте список distr, в котором будете хранить значения распределения
for k in range(0, n_exams + 1):
     # посчитайте кол-во вариантов выбрать k из n_exams
    choose = factorial(n_exams) / (factorial(k) * factorial(n_exams-k))
    # посчитайте вероятность
    prob = choose * failure_rate**k*(1-failure_rate)**(n_exams-k)
    # добавьте очередную вероятность prob в список distr
    distr.append(prob)
# построение гистограммы распределения вероятностей
plt.bar(range(0, n_exams + 1), distr)

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

Кратко:

В статье обсуждаются непрерывные случайные величины и их распределения.
- Непрерывное равномерное распределение: случайная величина принимает любые значения на определённом промежутке или даже на всей числовой прямой.
- Непрерывное равномерное распределение задаётся функцией плотности вероятности.
- Нормальное распределение: среднее и дисперсия определяют параметры нормального распределения.
- Правило трёх сигм: почти 100% значений попадают в пределы трёх стандартных отклонений от математического ожидания.
- Сравнение распределений: непрерывное равномерное и нормальное распределения имеют свои сходства и различия

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

• 10-метровый провод повреждён, но в какой именно точке — неизвестно. Видно, что на первых двух метрах всё в порядке, как и на последних трёх. Оставшаяся средняя часть провода скрыта за плинтусом. У электрика нет более точной информации о местоположении повреждения, то есть оно с равной вероятностью может располагаться в любом месте скрытого промежутка. Поэтому подходит равномерное распределение с a=2 и b=10−3=7. Значит, здесь X∼U[2,7].
• Николай ждёт звонка с 13:30 до 15:00. Звонок происходит в случайный момент, то есть с равной вероятностью от начала периода может пройти любое время. Здесь X∼U[13.5,15].
• Турист ждёт поезд в метро. Турист не знает города, поэтому считает, что поезда приходят в случайное время. По разговорам людей турист понял, что поезд в метро всегда приходит в течение 3 минут. Он с равной вероятностью может прождать любой промежуток времени до 3 минут. Здесь X∼U[0,3].

Значит, искомая вероятность равна площади закрашенного прямоугольника. Вычислим её: (4−2)⋅1/3=2/3​. Значит, вероятность того, что $X$ попадёт в интервал $,$ равна 2/3​. С этой вероятностью Никита управится со следующей задачей за период от 2 до 4 часов.

 

Напомним немного про обозначения промежутков. То, что величина X лежит в интервале от $2$ до 4, можно записать разными способами: 2≤X≤4 или X∈[2,4]. Это равнозначные обозначения. Их используют и для записи вероятностей. Например, вероятность того, что случайная величина X попадёт в промежуток, можно записать как P(X∈[2,4]). Как мы только что выяснили, P(X∈[2,4])=2/3.

 

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

• Средние значения и значения около них встречаются с высокой вероятностью.
• Чем дальше от среднего, тем меньше вероятность встретить такое значение.
• Вероятность очень маленьких или очень больших значений мала, но есть.

Фотография расставленных по интервалам роста шириной в 1 дюйм 175 студентов одного из колледжей (Blakeslee, 1914)

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

Изменение параметров нормального распределения

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

Кратко:

• Функция нормального распределения задается двумя параметрами: центром и стандартным отклонением.
• Функция распределения (кумулятивная функция) накапливает вероятность по мере роста значения x.
• Для непрерывного нормального распределения кумулятивная функция F(x) = P(X ≤ x) = P(X < x) = P(X ∈ (-∞, x)).
• График плотности вероятности нормального распределения имеет колоколообразную форму.
• Методы norm() и cdf() из пакета scipy.stats помогают рассчитать значения функции нормального распределения.
• Метод cdf() возвращает вероятность, что случайная величина примет значение, меньшее заданного аргумента.
• Вероятность попадания в конечный интервал рассчитывается путем вычитания меньшей кумулятивной вероятности из большей.
• Пример задачи: время реакции беспилотного автомобиля на препятствие распределено нормально, нужно найти вероятность попадания в интервал от 900 до 1100 миллисекунд

Рассчитать значения функции нормального распределения в Python помогут методы из пакета scipy.stats: norm() и cdf().

 

Слово norm в именах методов означает «нормальное распределение» (от англ. normal distribution).

 

Метод norm() задаёт нормальное распределение. На вход он получает два аргумента: математическое ожидание и стандартное отклонение.

 

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

 

Зададим в коде нормальное распределение N(1000,100^2) с помощью метода norm(1000, 100), а затем с помощью метода cdf() найдём вероятность получить 1000 или менее при таком распределении:

from scipy import stats as st
# математическое ожидание
mu = 1000
# стандартное отклонение
sigma = 100
# задаём нормальное распределение с математическим ожиданием 1000 
# и стандартным отклонением 100
distr = st.norm(1000, 100) 
# значение, для которого хотим найти вероятность
x = 1000
# считаем вероятность получить значение, равное x или меньше 
result = distr.cdf(x)
print(result)

0.5

from scipy import stats as st
# математическое ожидание
mu = 1000
# стандартное отклонение
sigma = 100
# задаём нормальное распределение с математическим ожиданием 1000 
# и стандартным отклонением 100
distr = st.norm(1000, 100) 
# значение, для которого хотим найти вероятность
x = 1150
# считаем вероятность получить значение больше x
result = 1 - distr.cdf(x)
print(result)

from scipy import stats as st
mu = 1000
sigma = 100
distr = st.norm(mu, sigma) 
x1 = 900
x2 = 1100
result = distr.cdf(x2) - distr.cdf(x1)
print(result)

from scipy import stats as st
# задайте нормальное распределение с нужными параметрами
mu = 720
sigma = 16900**0.5
print(sigma)
# задайте значения случайной величины - пороги длин броска для получения грамот и медалей
distr = st.norm(mu, sigma) 
x1 = 735
x2 = 750
print('Вероятность получить грамоту:',  1 - distr.cdf(x1))# рассчитайте вероятность 
print('Вероятность получить грамоту без медали:', distr.cdf(x2) - distr.cdf(x1)) # рассчитайте вероятность )
print('Вероятность получить медаль:',  1 - distr.cdf(x2))# рассчитайте вероятность )

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

Кратко:

• Стандартное нормальное распределение - это нормальное распределение с математическим ожиданием, равным 0, и стандартным отклонением, равным 1.
• Любое нормальное распределение можно с помощью простых арифметических операций привести к стандартному.
• Чтобы стандартизировать нормальное распределение, нужно: Вычесть из случайной величины её математическое ожидание μ - график распределения при этом сдвинется так, чтобы центр оказался в нуле.
• Разделить полученную случайную величину на стандартное отклонение σ.
• Стандартизация помогает решить задачи алгебраически.
• Z-значение - это количество стандартных отклонений, на которое некоторое число отклонилось от центра нормального распределения.
• Расчёт Z-значений позволяет сравнить, на сколько значения случайных величин с разными нормальными распределениями отклонились от своих математических ожиданий

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

 

Стандартное нормальное распределение — это нормальное распределение с математическим ожиданием, равным 0, и стандартным отклонением, равным 1. Его принято обозначать Z: Z∼N(0,1)

 

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

 

Чтобы стандартизировать нормальное распределение, нужно:

1. Вычесть из случайной величины её математическое ожидание μ — график распределения при этом сдвинется так, чтобы центр оказался в нуле.
2. Разделить полученную случайную величину на стандартное отклонение $.$ Стандартное отклонение получившегося распределения станет равным 1.

 

Запишем формулой:

 

X∼N(μ,σ^2)⇒(X−μ)/σ∼N(0,1)

 

Посмотрите на графики нормального распределения с μ=20,σ=3 и стандартного нормального распределения, которое получается из него стандартизацией:

 

Как привести к стандартному нормальному распределению случайную величину $Y\sim \mathbb{N}(145,144)?$

 

Правильный ответ

Вычесть 145 и разделить на 12

Математическое ожидание Y равно 145 — сначала его нужно отнять. Дисперсия Y равна $,$ значит, стандартное отклонение равно 144**0.5=12

 

Зачем стандартизировать распределения? Отличный вопрос! Есть класс задач, которые не решить без приведения нормального распределения к стандартному виду. Приведём пример.

 

В школе «Инклюзия» проводят соревнования по прыжкам в длину. Известно, что длины прыжков школьников и школьниц распределены нормально, но с разными параметрами: длина прыжка школьников — со средним 3.9 метра и стандартным отклонением 0.5 метра, а школьниц — со средним 3.5 метров и стандартным отклонением 0.7 метра. Организаторы соревнований хотят установить пороговую длину и выдавать сертификат всем, кто прыгнет дальше этого порога. При этом важно, чтобы вероятность получить сертификат была одинакова для всех детей независимо от пола.

 

В чём особенности этой задачи:

• даны два нормальных распределения;
• параметры распределений — разные;
• распределения связаны одним пороговым условием.

Переведём задачу на язык математики.

Распределение длин прыжков школьников: M∼N(3.9,0.5^2)

Распределение длин прыжков школьниц: F∼N(3.5,0.7^2)

Нужно найти число a, такое, что: P(M≥a)=P(F≥a)

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

 

Для решения нам пригодится новое понятие — Z-значение.

Z-значение

Z-значение — это количество стандартных отклонений, на которое некоторое число отклонилось от центра нормального распределения.

 

Например, если случайная величина X распределена так: X∼N(1000,100^2) Рассмотрим число $,$ оно больше среднего на $,$ что равно одному стандартному отклонению. Значит, Z-значение числа 1100 равно $.$

 

Теперь рассмотрим 900. Оно меньше среднего на $,$ это одно стандартное отклонение в отрицательную сторону. Значит, Z-значение числа 900 равно −1.

 

Z-значение любого значения нормально распределённой случайной величины можно получить по всё той же формуле: (X−μ)/σ∼Z

 

Проверим для тех же примеров:

Z-значение числа $:$ (1100−1000)/100=1;

Z-значение числа $:$ (900−1000)/100=−1

 

Для X∼N(1000,1002) чему будет равно Z-значение числа 950?

 

Правильный ответ

−1/2

Вычислим по формуле: (950−1000)/100=−50/100=−1/2. Значит, Z-значение числа 950 равно $.$ Другими словами, число 950 лежит левее математического ожидания, равного $,$ на половину стандартного отклонения, равного $.$

 

Возьмём два числа, симметричных относительно математического ожидания. Например, для X∼N(1000,100:2) это может быть пара 1100 и 900 или пара 1050 и 950. Числа в такой паре находятся на одинаковом расстоянии от центра распределения. Нормальное распределение симметрично, поэтому Z-значения чисел в такой паре будут равны по модулю, но противоположны по знаку.

 

Выше мы уже рассчитали, что для X∼N(1000,100:2) Z-значение числа 1100 равно $,$ а числа 900 равно −1. Z-значения равны по модулю и отличаются знаками.

 

Аналогично Z-значение числа 1050 равно 0.5, а числа 950 равно −0.5

 

Такие пары чисел, симметричных относительно среднего, отрубают одинаковые «хвосты» распределения слева и справа:

 

Если для нормального распределения X∼N(μ,σ:2) есть значения x0 и x1, такие что $,$ то Z-значение x0 будет равно минус Z-значению x1​.

 

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

 

Допустим, случайная величина X∼N(100,10^2) приняла значение x=95, а случайная величина Y∼N(180,4^2) приняла значение y=178. Их Z-значений равны, и по ним сразу можно понять, что и x, и y отклонились на половину стандартного отклонения влево от математического ожидания:

 

$.$

 

Эти значения отрубили одинаковые «хвосты» от своих распределений. То есть вероятность того, что X окажется меньше (или больше) 95 равна вероятности того, что Y окажется меньше (или больше) 178: P(X≤95)=P(Y≤178); P(X≥95)=P(Y≥178).

 

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

 

Нам известно, что F∼N(3.5,0.7^2) и M∼N(3.9,0.5^2).

 

Нужно найти число $a$, такое, что: P(M≥a)=P(F≥a).

 

Число $a$ должно отрубать справа одинаковый «хвост» от обоих распределений. Это значит, что Z-значение числа a должно быть одинаковым для F и для M.

 

На графике это выглядит так:

 

Стандартизируем оба нормальных распределения, затем учтём, что получающиеся Z-значения для числа $a$ должны быть одинаковыми и для распределений F и M. Получим уравнение:

 

(a−3.5)/0.7=(a−3.9)/0.50.

​

Это уравнение с одной неизвестной, и значит его можно решить. Домножим обе части равенства на 0.5 и на 0.7:

 

0.7⋅(a−3.9)=0.5⋅(a−3.5)

 

Раскроем скобки и приведём подобные слагаемые:

 

0.7a−2.73=0.5a−1.75;

0.2a=0.98;

a=4.9

 

Число $a$ найдено! Если установить порог прыжка в 4.9 метра, сертификат получит одинаковый процент и школьников, и школьниц.

 

Найдём этот процент. Используем метод cdf().

rom scipy import stats as st
# зададим нормальное распределение для школьниц:
mu_F = 3.5
sigma_F = 0.7
distr_F = st.norm(mu_F, sigma_F) 
# зададим нормальное распределение для школьников:
mu_M = 3.9
sigma_M = 0.5
distr_M = st.norm(mu_M, sigma_M) 
# найденная пороговая длина прыжка a, за прыжок дальше которой выдаётся сертификат:
a = 4.9
# найдём и выведем на экран вероятности правых хвостов для a,
# округлив их до 5 цифр после запятой:
print('Вероятность того, что случайная школьница получит сертификат:',
				round(1- distr_F.cdf(a),5))
print('Вероятность того, что случайный школьник получит сертификат:',
				round(1- distr_M.cdf(a), 5))

 

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

 

У компании International два офиса: один — в Ванкувере, другой — в Лондоне. В обоих офисах одинаковое количество сотрудников, а рабочий день с 9 до 18. Когда в Ванкувере 9 утра, в Лондоне — 17 вечера.

 

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

 

Менеджеры выяснили, что время прихода на работу в Ванкувере распределено нормально со средним в 9:00 и стандартным отклонением в 10 минут, а в Лондоне время ухода с работы распределено со средним в 18:00 и стандартным отклонением в 15 минут.

 

Какое время по часовому поясу Ванкувера стоит выбрать для начала онлайн-игры?

 

Первая подсказка

Переведём время в минуты и зададим нормальные распределения на единой числовой оси. Примем за ноль время 9 утра по Ванкуверу. Тогда время прихода сотрудников в этот офис распределено так: $.$

Разница между 9 утра в Ванкувере и 18 вечера в Лондоне — час. Тогда время ухода сотрудников в лондонском офисе распределено так: $.$

 

Вторая подсказка

Нужно найти число $a,$ такое что: P(V≤a)=P(L≥a). Слева в этом уравнении стоит вероятность левого «хвоста», а справа — правого. Чтобы корректно составить уравнение, нам нужны вероятности одинаковых хвостов. На наше счастье, Z-значения симметричны, поэтому при составлении уравнения просто заменим знак одной из частей равенства. Получим уравнение:

(a−0)/10=−(a−60)/15​.

 

Правильный ответ

9:24

Похоже, это оптимальное время для начала игры.

 

Решение

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

Время прихода на работу сотрудников в Ванкувере распределено так: V∼N(0,10^2). 

Время ухода сотрудников в Лондоне распределено так: L∼N(60,15^2)

Найдём число $a,$ такое что: P(V≤a)=P(L≥a).

 

Слева в этом уравнении стоит вероятность левого «хвоста», а справа — правого. Чтобы корректно составить уравнение, нам нужны вероятности одинаковых хвостов. На наше счастье, Z-значения симметричны, поэтому при составлении уравнения просто заменим знак одной из частей равенства. Получим уравнение:

(a−0)/10=−(a−60)/15

Решим его:

600−10a=15a;

a=24

 

За ноль мы принимали 9 часов Ванкувера. Значит, наиболее подходящее время начала онлайн-игры — 9:24 по Ванкуверу, в это время в Лондоне будет 17:24.

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

Кратко:

• Статья посвящена распределению Пуассона и его применению в различных задачах.
• Распределение Пуассона описывает вероятности возможного числа событий в процессах с известной интенсивностью.
• Параметры λ (лямбда) и математическое ожидание совпадают, а дисперсия равна λ.
• Вероятность P(Y=k) для k событий в заданном периоде можно рассчитать по формуле P(Y=k) = λ^kk!· e^λ.
• График распределения Пуассона показывает вероятности для разных значений k.
• Функция распределения Пуассона позволяет находить кумулятивные вероятности и "хвосты" распределения.
• Распределения Пуассона и биномиальное распределение могут быть аппроксимированы друг другом

import matplotlib.pyplot as plt
from math import factorial, exp
 # параметр лямбда (используем lmbd, так как lambda зарезервировано в Python):
lmbd = 3
 distr = []
# для всех чисел от 0 до 19
for k in range(20):
    # рассчитаем вероятности по формуле распределения Пуассона:
    p = (lmbd ** k) / (factorial(k) * exp(lmbd))
    distr.append(p)
 # зададим график из полученных значений с о-образными точками синего цвета
plt.plot(range(20), distr, 'bo')
# используем числа от 0 до 19 как подписи для горизонтальной оси:
plt.xticks(range(20))
# отобразим график:
plt.show()

import matplotlib.pyplot as plt
from scipy.stats import poisson
# параметр лямбда:
lmbd = 3
distr = []
# для всех чисел от 0 до 19
for k in range(20):
    # рассчитаем вероятности этого значения с помощью метода pmf():
    p = poisson.pmf(k, lmbd)
    distr.append(p)
 # зададим график из полученных значений с о-образными точками синего цвета
plt.plot(range(20), distr, 'bo')
# используем числа от 0 до 19 как подписи для горизонтальной оси:
plt.xticks(range(20))
# отобразим график:
plt.show()

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson
# параметр лямбда:
lmbd = 3
distr = []
# для значений от 0 до 19
for x in range(20):

Кумулятивная вероятность помогает находить «хвосты» распределения. Для распределения Пуассона это ответы на вопросы типа «какова вероятность, что за промежуток времени произойдёт меньше (или больше) заданного числа событий?»

 

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

Поскольку параметр λ задаёт интенсивность за определённый период времени, его легко масштабировать на другие временные отрезки. Если в среднем за неделю бывает 3 северных сияния, то за две недели — $6$, а за четыре — 12. Значит, рассчитать вероятности количества сияний для промежутка в две недели поможет распределение Пуассона с параметром λ=6. Очень удобно!

from scipy.stats import poisson
# минимальное необходимое количество:
min_num = 50
# количество дней в июне, июле и августе:
days = 30 + 31 + 31 
# рассчитаем параметр лямбда для этого количества дней:
weeks = days / 7 
lmbd = 4 * weeks
# параметр лямбда может получиться нецелым, это допустимо
print(f'Вероятность правого "хвоста" распределения Пуассона от {min_num}',
      f'при lambda={round(lmbd,1)} равна {1 - poisson.cdf(min_num - 1, lmbd)}'
) 

Вероятность правого "хвоста" распределения Пуассона от 50 при lambda=52.6 равна 0.657101690371793 

from scipy.stats import poisson

guests_threshold = 7 # укажите лимит посетителей
lmbd = 20/4 # укажите или рассчитайте параметр лямбда

# выведите на экран получившуюся вероятность
print(1-poisson.cdf(guests_threshold, lmbd))

Аппроксимация биномиального распределения распределением Пуассона

Кратко:

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

Аппроксимация биномиального распределения нормальным

Кратко:

• Биномиальное распределение описывает множество переменных и часто возникает задача найти вероятность "хвоста" распределения.
• Нормальная аппроксимация удобна и быстра для вычислений.
• Аппроксимация выполняется, если условие выполнено: математическое ожидание биномиального распределения в пределах ± три стандартных отклонения.
• Чем ближе p к 0,5, тем меньше должно быть n для выполнения условия аппроксимации.
• Чем ближе p к 0,5, тем меньше ошибки аппроксимации.
• Аппроксимация нормальным распределением подходит для поиска вероятности "хвоста" биномиального распределения

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

 

from math import sqrt
import matplotlib.pyplot as plt
from scipy.stats import binom, norm
 
# зададим функцию, которая выдаёт <num> чисел,
# равномерно распределённых от <start> до <stop>
# (понадобится для визуализации нормального распределения)
def linspace(start, stop, num):
    step = (stop - start) / (num - 1)
    result = []
    for i in range(num):
        result.append(start + step * i)
    return result
 
# задаём параметры биномиального распределения
n = 50
p = 0.22

# границы промежутка, в котором будем строить аппроксимацию -
# от n*p минус 4 ст.отклонения бином.распределения
# до n*p плюс 4 ст.отклонения бином.распределения
# (там, где вероятности значений почти не равны нулю)
left = int(n * p - 4 * sqrt(n * p * (1 - p)))
right = int(n * p + 4 * sqrt(n * p * (1 - p)))
 
# зададим набор значений биномиальной величины на выбранном промежутке:
# все неотрицательные целые числа от left до right
x = list(range(max(0, left), right))
# рисуем график биномиального распределения
ax = plt.subplot()
# используем синие о-образные маркеры и размер маркера 5
ax.plot(x, binom.pmf(x, n, p), 'bo', ms=5,
	      label='вероятности биномиального распределения')

# зададим мат.ожидание и ст.отклонение нормального распределения равными
# мат.ожиданию и ст.отклонению биномиального распределения
mu = n * p
sigma = sqrt(n * p * (1 - p))
 
# рисуем график нормального распределения:
# получим список из 100 чисел от left до right
x = linspace(left, right, 100)
# используем красную линию ширины 3 и прозрачности 0.3
ax.plot(x, norm.pdf(x, mu, sigma), 'r-', lw=3, alpha=0.3, label='плотность вероятности нормального распределения')
 
# зададим подписи по оси x в пределах ± 3 ст.отклонений от мат.ожидания
# (последний аргумент в методе range - интервал между подписями)
x_ticks = list(range(int(mu - 3 * sigma), int(mu + 3 * sigma), 1))
ax.set_xticks(x_ticks)

# настраиваем оптимальное расположение легенды и выводим график
ax.legend(loc='best')
plt.show()

# выводим, выполняется ли наше условие аппроксимации:
# лежит ли mu ± 3 sigma в промежутке от 0 до n:
if (mu - (3 * sigma) >= 0) and ((mu + (3 * sigma) <= n)):
    print('Условие выполняется')
else:
    print('Условие не выполняется')

• n=20, p=0.5:

Условие выполняется 

• n=50, p=0.03:

Условие не выполняется 

• Если оставить ту же вероятность p=0.03, но увеличить $n$ в десять раз, условие уже будет выполняться и аппроксимация будет хорошей.
  n=500, p=0.03:

Условие выполняется 

Условие не выполняется 

Условие выполняется 

Найдём эту площадь с помощью знакомого вам метода cdf():

from math import sqrt
from scipy import stats as st

# задаём параметры биномиального распределения
n = 45
p = 0.4

# зададим мат.ожидание и ст.отклонение нормального распределения равными
# мат.ожиданию и ст.отклонению биномиального распределения
mu = n * p
sigma = sqrt(n * p * (1 - p))

# задаём нормальное распределение с математическим ожиданием 18 
# и стандартным отклонением ~3.29
distr = st.norm(mu, sigma) 

# границы интервала
y_1 = 15.5
y_2 = 16.5

# считаем вероятность получить значение от y1 до y2
result = distr.cdf(y_2) - distr.cdf(y_1)
print(result)

from math import sqrt
from scipy import stats as st

# задаём параметры биномиального распределения
n = 45
p = 0.4

# зададим мат.ожидание и ст.отклонение нормального распределения равными
# мат.ожиданию и ст.отклонению биномиального распределения
mu = n * p
sigma = sqrt(n * p * (1 - p))

# задаём нормальное распределение с математическим ожиданием 18 
# и стандартным отклонением ~3.29
distr = st.norm(mu, sigma) 

# найдём суммарную ошибку аппроксимации:
error_values = []
for i in range(n + 1):
    error = abs(st.binom.pmf(i, n, p) - (distr.cdf(i + 0.5) - distr.cdf(i - 0.5)))
    error_values.append(error)
 
print('Средняя разница вероятностей на значение:', sum(error_values) / len(error_values))
print('Максимальная разница вероятностей:', max(error_values))

from math import sqrt
from scipy import stats as st

# задаём параметры биномиального распределения
n = 1500000
p = 0.029

# зададим мат.ожидание и ст.отклонение нормального распределения равными
# мат.ожиданию и ст.отклонению биномиального распределения
mu = n * p
sigma = sqrt(n * p * (1 - p))

# задаём нормальное распределение
distr = st.norm(mu, sigma) 

# считаем вероятность хвоста слева от 43 тысяч
result = distr.cdf(43000)
print(result)

from math import sqrt
from scipy import stats as st

n = 550000 # задайте параметры 
p = 0.114 # биномиального распределения

mu = n * p # рассчитайте параметры 
sigma = sqrt(n * p * (1-p)) # нормального распределения

distr = st.norm(mu, sigma)  # задайте нормальное распределение

result = 1- distr.cdf(63000) # найдите вероятность
print(result)

Следующая тема:  САД. Проверка гипотез

Вернуться в раздел: Статистический анализ данных

Вернуться в оглавление: Я.Практикум