Как я могу эффективно вычислить функцию биномиального кумулятивного распределения?
Скажем, что я знаю, что вероятность "успеха" - это P. Я выполняю тест N раз, и я вижу успехи S. Тест сродни бросанию неравномерно взвешенной монеты (возможно, голова - это успех, хвосты - провал).
Я хочу знать приблизительную вероятность увидеть либо успехи S, либо несколько успехов, менее вероятных, чем S успехов.
Так, например, если P равно 0,3, N равно 100, и я получаю 20 успехов, я ищу вероятность получить 20 или меньше успехов.
Если, с другой стороны, P равно 0.3, N равно 100, и я получаю 40 успехов, я ищу вероятность получить 40 наших дальнейших успехов.
Я знаю, что эта проблема связана с поиском области под биномиальной кривой:
- Моя математика-фу не связана с задачей перевода этих знаний в эффективный код.
- Хотя я понимаю, что биномиальная кривая даст точный результат, у меня создается впечатление, что она по сути своей неэффективна. Быстрого метода вычисления приблизительного результата будет достаточно.
Я должен подчеркнуть, что это вычисление должно быть быстрым и в идеале должно быть определено при стандартном вычислении с плавающей запятой 64 или 128 бит.
Я ищу функцию, которая принимает P, S и N - и возвращает вероятность. Поскольку я более знаком с кодом, чем с математической нотацией, я бы предпочел, чтобы в любых ответах использовался псевдокод или код.
Ответы
Ответ 1
Точное биномиальное распределение
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 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
Ответ 2
Я думаю, вы хотите оценить неполную бета-функцию .
Там хорошая реализация с использованием представления непрерывной доли в "Численных рецептах в C", глава 6: "Специальные функции".
Ответ 3
Я не могу полностью ручаться за эффективность, но Scipy имеет модуль для этого
from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)
Ответ 4
Из части вашего вопроса "получение по крайней мере S головок" вам нужна кумулятивная функция биномиального распределения. См. http://en.wikipedia.org/wiki/Binomial_distribution для уравнения, которое описывается как "регулируемая незавершенная бета-функция" (как уже было сказано). Если вы просто хотите рассчитать ответ, не выполняя все решения самостоятельно, Научная библиотека GNU предоставляет функцию: gsl_cdf_binomial_P и gsl_cdf_binomial_Q.
Ответ 5
Эффективный и, что более важно, численный стабильный алгоритм существует в области кривых Безье, используемых в компьютерном проектировании. Он называется алгоритмом де Кастеляу, используемым для оценки полиномов Бернштейна, используемых для определения кривых Безье.
Я считаю, что мне разрешено только одну ссылку за ответ, поэтому начните с Википедия - Политимы Бернштейна
Обратите внимание на очень тесную связь между Биномиальным распределением и полиномами Бернштейна. Затем перейдите по ссылке на алгоритм де Кастеляу.
Позволяет сказать, что я знаю, что вероятность бросать головы с определенной монетой равна P. Какова вероятность того, что я брошу монета T раз и получение по крайней мере S?
- Установить n = T
- Установить бета [i] = 0 для я = 0,... S - 1
- Установить бета [i] = 1 для я = S,... T
- Установите t = p
- Оцените B (t), используя de Casteljau
или не более S головок?
- Установить n = T
- Установить бета [i] = 1 для я = 0,... S
- Установить бета [i] = 0 для я = S + 1,... T
- Установите t = p
- Оцените B (t), используя de Casteljau
Открытый исходный код, вероятно, уже существует. Кривые NURBS (неравномерные Rational B-сплайновые кривые) являются обобщением кривых Безье и широко используются в CAD. Попробуйте openNurbs (лицензия очень либеральная) или не получившая Open CASCADE (несколько менее либеральная и непрозрачная лицензия). Оба набора инструментальных средств находятся на С++, хотя существуют привязки IIRC,.NET.
Ответ 6
DCDFLIB Project имеет функции С# (обертки вокруг кода C) для оценки многих функций CDF, включая биномиальное распределение. Здесь вы можете найти исходный код C и FORTRAN здесь. Этот код хорошо протестирован и точен.
Если вы хотите написать свой собственный код, чтобы избежать зависимости от внешней библиотеки, вы можете использовать обычное приближение к биномии, упомянутому в других ответах. Вот несколько примечаний к насколько хороша аппроксимация при различных обстоятельствах. Если вы идете по этому маршруту и нуждаетесь в коде для вычисления нормального CDF, здесь код Python для этого. Это всего лишь около дюжины строк кода и легко переносится на любой другой язык. Но если вам нужна высокая точность и эффективный код, вам лучше использовать сторонний код, например DCDFLIB. Несколько человеко-лет вступили в производство этой библиотеки.
Ответ 7
Если вы используете Python, не нужно сам его кодировать. Scipy заставил вас прикрыться:
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
Ответ 8
Я был в проекте, где нам нужно было вычислить биномиальный 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 = 10000.
Я не буду вдаваться в полный вывод этого алгоритма, потому что stackoverflow не поддерживает MathJax, но его толчок сначала идентифицирует следующую эквивалентность:
- Для всех k > 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_ {х = я-1}. Это означает, что мы можем делать все в одном цикле с использованием аккумуляторов, и нам не нужно вычислять какие-либо факториалы!
Причина, по которой большая часть вычислений выполняется в лог-пространстве, заключается в улучшении численной устойчивости полиномиальных членов, т.е. p^x
и (1-p)^(1-x)
имеют потенциал быть чрезвычайно большим или крайне малым, что может вызвать ошибки вычислений.
Ответ 9
Попробуйте этот, используемый в GMP. Другая ссылка this.