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