Оценка вероятности по выборочным данным

границы доверительного интервала,  вычисление на Python

Получение оценки вероятности по результатам наблюдений одна из распространенных задач в статистическом анализе. Рассмотрим следующий пример*:

С 1871 по 1900 г. в Швейцарии родились 1359671 мальчик и 1285086 девочек (см. Polya,Handbuch der biol. Arbeitsmethoden, S.742). Что можно сказать о величине вероятности рождения мальчика?

Формулировка задачи неявно предполагает, что вероятностные условия для событий рождения мальчика или девочки сохраняются неизменными для всех (1359671+1285086) рождений. Это означает, что представленные в задаче данные можно рассматривать как результат осуществления (1359671+1285086) экспериментов по схеме Бернулли с некоторыми вероятностями выпадения "орла" \(p\) и "решки" \(q\). Пусть "орлом" будет появление мальчика, а "решкой"  — появление девочки.

Если \(k\) число рождений мальчиков из \(n\) рождений, то на основании теоремы Бернулли частота \(\nu=k/n\) по мере увеличения числа экспериментов \(n\) стремится к искомой вероятности \(p\).

Чтобы построить доверительный интервал для оценки вероятности воспользуемся интегральной теоремой Муавра-Лапласа: $$
P(\left|\nu - p\right|<\Delta)\approx\dfrac{1}{\sqrt{2\pi}}\int\limits_{-\Delta\sqrt{\dfrac{n}{pq}}}^{\Delta\sqrt{\dfrac{n}{pq}}}e^{-\dfrac{x^2}{2}}dx=2\Phi_0(\Delta\sqrt{\dfrac{n}{pq}}),
$$
где \(\Phi_0(z)=\dfrac{1}{\sqrt{2\pi}}\int\limits_0^ze^{-\dfrac{x^2}{2}}dx\).

Приведенное выражение, однако, обладает тем практическим недостатком, что содержит неизвестные параметры \(p\), \(q\). Однако этот недостаток можно обойти...

Оценим произведение \(pq\), учитывая, что \(q=1-p\). Рассмотрим \(f(p)=p(1-p)\). Данная функция имеет максимум на сегменте \([0,1]\) равный \(1/4\). Таким образом, \(pq\leq1/4\).

Следовательно, принимая во внимание монотонность функции \(\Phi_0\) и продолжая начатое выше неравенство, получим $$
P(\left|\nu - p\right|<\Delta)\approx2\Phi_0(\Delta\sqrt{\dfrac{n}{pq}})\geq
2\Phi_0(2\Delta\sqrt n)
$$

Таким образом, задавшись доверительной вероятностью \(\alpha\) можно найти соответствующий интервал (с центром \(\nu\)), которому с вероятностью \(\alpha\) принадлежит искомая вероятность рождения мальчика. Радиус (т.е. половина его длины) доверительного интервала  может быть найден как решение уравнения: \(2\Phi_0(2\Delta\sqrt n)=\alpha\).

Оценка вероятности в схеме Бернулли по результатам эксперимента: Python-реализация

Воспользуемся статистической библиотекой, входящей в пакет SciPy. К сожалению, в этой библиотеке отсутствует необходимая нам функция \(\Phi_0(x)\) (что вообще говоря не удивительно, поскольку  \(\Phi_0(x)\) функция весьма специфичного вида), однако
определена функция \(F(x)=\dfrac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^xe^{-\dfrac{x^2}{2}}dx\), являющаяся функцией стандартного нормального распределения.

В силу ее определения \(\Phi_0(x)\) может быть переписано в виде  \(\Phi_0(x)=F(x)-F(0)\). Таким образом, исходное уравнение
\(2\Phi_0(2\Delta\sqrt n)=\alpha\) примет вид \(2(F(2\Delta\sqrt n)-F(0))=\alpha\). Для \(F(x)\) в библиотеке SciPy определена обратная функция, что позволяет решить последнее уравнение и найти \(\Delta\).  В представленном ниже коде z.cdf(x) — это \(F(x)\), а z.ppf(x) — \(F^{-1}(x)\) — обратная функция к \(F(x)\), т.е. \(F^{-1}(F(x))=x\).
Следовательно,
$$
\Delta=\dfrac{F^{-1}(\alpha/2+F(0))}{2\sqrt n},\qquad (\mbox{к сведению,}\, F(0)=1/2),
$$
а соответствующий код вычислений выглядит следующим образом:

# -*- coding: utf-8 -*-
from math import sqrt     #Доступ к функции вычисления квадратного корня
import scipy.stats as st  #Подключение статистической библиотеки
n1 = 1359671 		  #Количество рожденных мальчиков
n2 = 1285086		  #Количество рожденных девочек
n = n1+n2			  #Общее количество детей
alpha = 0.95    		  #Доверительная вероятность
z = st.norm    	#Объект, содержащий функции, связанные с нормальным распред.
Delta = z.ppf(alpha/2+z.cdf(0))/(2.0*sqrt(n)) #радиус доверительного интервала
LowBoundP = float(n1)/n-Delta #Нижняя граница для вероятности
UpperBoundP = float(n1)/n+Delta #Верхняя граница для вероятности
print 'Искомый интервал:[',LowBoundP,',',UpperBoundP,']' #Вывод результата

В итоге, получим следующий результат:

Искомый интервал:[ 0.513497944732 , 0.514703133855 ]

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


*Задача взята из [Ван дер Варден Б.Л. Математическая статистика. ИЛ., 1960. 435 с.]

blog comments powered by Disqus