Как я могу эффективно рассчитать биномиальную кумулятивную функцию распределения?

StackOverflow https://stackoverflow.com/questions/1095650

Вопрос

Допустим, я знаю, что вероятность «успеха» равна P.Я запускаю тест N раз и вижу S успехов.Тест подобен подбрасыванию монеты разного веса (возможно, орел – успех, решка – провал).

Я хочу знать приблизительную вероятность увидеть либо S успехов, либо ряд успехов, менее вероятных, чем S успехов.

Например, если P равно 0,3, N равно 100 и я получил 20 успехов, я ищу вероятность получения 20. или меньше успехи.

Если, с другой стороны, P равно 0,3, N равно 100 и я получу 40 успехов, я ищу вероятность получить 40 наших дополнительных успехов.

Однако я знаю, что эта проблема связана с поиском площади под биномиальной кривой:

  1. Моя математика не справляется с задачей перевода этих знаний в эффективный код.
  2. Хотя я понимаю, что биномиальная кривая даст точный результат, у меня складывается впечатление, что она по своей сути неэффективна.Достаточно быстрого метода расчета приблизительного результата.

Я должен подчеркнуть, что эти вычисления должны быть быстрыми и в идеале должны определяться с помощью стандартных 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.На создание этой библиотеки ушло несколько человеко-лет.

Пытаться Вот этот, используемый в GMP.Еще одна ссылка этот.

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.
Лицензировано под: CC-BY-SA с атрибуция
Не связан с StackOverflow
scroll top