Как я могу эффективно рассчитать биномиальную кумулятивную функцию распределения?
-
11-09-2019 - |
Вопрос
Допустим, я знаю, что вероятность «успеха» равна P.Я запускаю тест N раз и вижу S успехов.Тест подобен подбрасыванию монеты разного веса (возможно, орел – успех, решка – провал).
Я хочу знать приблизительную вероятность увидеть либо S успехов, либо ряд успехов, менее вероятных, чем S успехов.
Например, если P равно 0,3, N равно 100 и я получил 20 успехов, я ищу вероятность получения 20. или меньше успехи.
Если, с другой стороны, P равно 0,3, N равно 100 и я получу 40 успехов, я ищу вероятность получить 40 наших дополнительных успехов.
Однако я знаю, что эта проблема связана с поиском площади под биномиальной кривой:
- Моя математика не справляется с задачей перевода этих знаний в эффективный код.
- Хотя я понимаю, что биномиальная кривая даст точный результат, у меня складывается впечатление, что она по своей сути неэффективна.Достаточно быстрого метода расчета приблизительного результата.
Я должен подчеркнуть, что эти вычисления должны быть быстрыми и в идеале должны определяться с помощью стандартных 64- или 128-битных вычислений с плавающей запятой.
Я ищу функцию, которая принимает P, S и N и возвращает вероятность.Поскольку я больше знаком с кодом, чем с математическими обозначениями, я бы предпочел, чтобы в любых ответах использовался псевдокод или код.
Решение
Точное биномиальное распределение
def factorial(n):
if n < 2: return 1
return reduce(lambda x, y: x*y, xrange(2, int(n)+1))
def prob(s, p, n):
x = 1.0 - p
a = n - s
b = s + 1
c = a + b - 1
prob = 0.0
for j in xrange(a, c + 1):
prob += factorial(c) / (factorial(j)*factorial(c-j)) \
* x**j * (1 - x)**(c-j)
return prob
>>> prob(20, 0.3, 100)
0.016462853241869437
>>> 1-prob(40-1, 0.3, 100)
0.020988576003924564
Нормальная оценка, подходит для больших n
import math
def erf(z):
t = 1.0 / (1.0 + 0.5 * abs(z))
# use Horner's method
ans = 1 - t * math.exp( -z*z - 1.26551223 +
t * ( 1.00002368 +
t * ( 0.37409196 +
t * ( 0.09678418 +
t * (-0.18628806 +
t * ( 0.27886807 +
t * (-1.13520398 +
t * ( 1.48851587 +
t * (-0.82215223 +
t * ( 0.17087277))))))))))
if z >= 0.0:
return ans
else:
return -ans
def normal_estimate(s, p, n):
u = n * p
o = (u * (1-p)) ** 0.5
return 0.5 * (1 + erf((s-u)/(o*2**0.5)))
>>> normal_estimate(20, 0.3, 100)
0.014548164531920815
>>> 1-normal_estimate(40-1, 0.3, 100)
0.024767304545069813
Оценка Пуассона:Хорошо для больших n и маленьких p
import math
def poisson(s,p,n):
L = n*p
sum = 0
for i in xrange(0, s+1):
sum += L**i/factorial(i)
return sum*math.e**(-L)
>>> poisson(20, 0.3, 100)
0.013411150012837811
>>> 1-poisson(40-1, 0.3, 100)
0.046253037645840323
Другие советы
Я работал над проектом, в котором нам нужно было иметь возможность вычислить биномиальный CDF в среде, в которой не были определены факториал или гамма-функция.Мне потребовалось несколько недель, но в итоге я придумал следующий алгоритм, который точно рассчитывает CDF (т.аппроксимация не требуется).Python, по сути, так же хорош, как и псевдокод, верно?
import numpy as np
def binomial_cdf(x,n,p):
cdf = 0
b = 0
for k in range(x+1):
if k > 0:
b += + np.log(n-k+1) - np.log(k)
log_pmf_k = b + k * np.log(p) + (n-k) * np.log(1-p)
cdf += np.exp(log_pmf_k)
return cdf
Производительность масштабируется с помощью x.Для малых значений x это решение примерно на порядок быстрее, чем scipy.stats.binom.cdf
, с аналогичной производительностью около x=10 000.
Я не буду вдаваться в полное описание этого алгоритма, поскольку stackoverflow не поддерживает MathJax, но суть его заключается в том, чтобы сначала определить следующую эквивалентность:
- Для всех к > 0
sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])
Что мы можем переписать как:
sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k
или в пространстве журнала:
np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)
Поскольку CDF представляет собой сумму PMF, мы можем использовать эту формулировку для расчета биномиального коэффициента (логарифм которого равен b
в приведенной выше функции) для PMF_{x=i} из коэффициента, который мы рассчитали для PMF_{x=i-1}.Это означает, что мы можем сделать все в одном цикле, используя аккумуляторы, и нам не нужно вычислять факториалы!
Причина, по которой большинство вычислений выполняется в логарифмическом пространстве, заключается в улучшении численной стабильности полиномиальных членов, т.е. p^x
и (1-p)^(1-x)
потенциально могут быть чрезвычайно большими или чрезвычайно маленькими, что может вызвать вычислительные ошибки.
РЕДАКТИРОВАТЬ:Это новый алгоритм?Я время от времени ковырялся в этом с тех пор, как опубликовал это, и все чаще задаюсь вопросом, стоит ли мне написать это более формально и отправить в журнал.
Я думаю, вы хотите оценить неполная бета-функция.
В главе 6 книги «Числовые рецепты на языке C» есть хорошая реализация с использованием представления непрерывной дроби:«Специальные функции».
Я не могу полностью ручаться за эффективность, но у Сципи есть модуль для этого
from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)
Эффективный и, что более важно, численный устойчивый алгоритм существует в области Кривые Безье используется в компьютерном проектировании.Это называется алгоритм де Кастельжо используется для оценки Полиномы Бернштейна используется для определения кривых Безье.
Я считаю, что мне разрешена только одна ссылка на ответ, поэтому начните с Википедия - Полиномы Бернштейна
Обратите внимание на очень тесную связь между биномиальным распределением и полиномами Бернштейна.Затем перейдите по ссылке на алгоритм де Кастельжо.
Допустим, я знаю, что вероятность выпадения орла на конкретной монете равна P.Какова вероятность того, что я бросаю время монеты и получает хотя бы головы?
- Установите n = Т
- Установите beta[i] = 0 для i = 0,...С - 1
- Установите beta[i] = 1 для i = S,...Т
- Установить т = р
- Оцените B(t) с помощью де Кастельжо.
или максимум S-головки?
- Установите n = Т
- Установите beta[i] = 1 для i = 0,...С
- Установите beta[i] = 0 для i = S + 1,...Т
- Установить т = р
- Оцените B(t) с помощью де Кастельжо.
Открытый исходный код, вероятно, уже существует. NURBS-кривые (Неоднородные рациональные B-сплайновые кривые) являются обобщением кривых Безье и широко используются в САПР.Попробуйте openNurbs (очень либеральная лицензия) или провалите Open CASCADE (несколько менее либеральная и непрозрачная лицензия).Оба инструментария написаны на C++, однако существуют привязки IIRC и .NET.
Если вы используете Python, вам не нужно писать код самостоятельно.Сципи тебя прикрыл:
from scipy.stats import binom
# probability that you get 20 or less successes out of 100, when p=0.3
binom.cdf(20, 100, 0.3)
>>> 0.016462853241869434
# probability that you get exactly 20 successes out of 100, when p=0.3
binom.pmf(20, 100, 0.3)
>>> 0.0075756449257260777
Из части вашего вопроса «получить как минимум S голов» вам нужна кумулятивная биномиальная функция распределения.Видеть http://en.wikipedia.org/wiki/Binomial_distribution для уравнения, которое описывается как «регуляризованная неполная бета-функция» (как уже было сказано).Если вы просто хотите вычислить ответ, не реализуя все решение самостоятельно, Научная библиотека GNU предоставляет функцию:gsl_cdf_binomial_P и gsl_cdf_binomial_Q.
А Проект DCDFLIB имеет функции C# (обертки вокруг кода C) для оценки многих функций CDF, включая биномиальное распределение.Вы можете найти исходный код C и FORTRAN. здесь.Этот код хорошо протестирован и точен.
Если вы хотите написать свой собственный код, чтобы избежать зависимости от внешней библиотеки, вы можете использовать обычное приближение к биному, упомянутому в других ответах.Вот некоторые заметки по насколько хорошее приближение при различных обстоятельствах.Если вы пойдете по этому пути и вам понадобится код для вычисления обычного CDF, вот Код Python за это.Это всего лишь около дюжины строк кода, и его можно легко перенести на любой другой язык.Но если вам нужна высокая точность и эффективность кода, лучше использовать сторонний код, например DCDFLIB.На создание этой библиотеки ушло несколько человеко-лет.
import numpy as np
np.random.seed(1)
x=np.random.binomial(20,0.6,10000) #20 flips of coin,probability of
heads percentage and 10000 times
done.
sum(x>12)/len(x)
The output is 41% of times we got 12 heads.