Comment puis-je calculer efficacement la fonction de distribution cumulative binomiale?

https://stackoverflow.com/questions/1095650

11-09-2019
|

Question

Disons que je sais que la probabilité d'un « succès » est P. je lance le test N fois, et je vois des succès S. Le test est semblable à une pièce de monnaie de façon inégale lancer pondérée (peut-être la tête est un succès, la queue est un échec).

Je veux savoir la probabilité approximative de voir soit succès S, ou un certain nombre de succès moins susceptibles que les réussites.

Ainsi, par exemple, si P est de 0,3, N est 100, et je reçois 20 succès, je cherche la probabilité d'obtenir 20 ou moins succès.

Si, sur l'autre hadn, P est de 0,3, N est 100, et je reçois 40 succès, je cherche la probabilité d'obtenir 40 plus nos succès.

Je suis conscient du fait que ce problème est lié à trouver l'aire sous une courbe binomiale, cependant:

Mon maths-fu n'est pas à la tâche de traduire ces connaissances en code efficace
Bien que je comprenne une courbe binomiale donnerait un résultat exact, je l'impression que ce serait intrinsèquement inefficace. Une méthode rapide pour calculer un résultat approximatif suffit.

Je tiens à souligner que ce calcul doit être rapide, et devrait idéalement être déterminable à la norme 64 ou 128 bits calcul en virgule flottante.

Je cherche une fonction qui prend P, S et N - et renvoie une probabilité. Comme je suis plus familier avec le code de la notation mathématique, je préfère que les réponses emploient pseudo-code ou code.

La solution

Distribution binomiale exacte

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

Estimation normale, bonne pour les grandes 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

Poisson Estimation: Bon pour n grand et petit 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

Autres conseils

J'étais sur un projet où nous devions être en mesure de calculer le CDF binomiale dans un environnement qui n'a pas eu une fonction factoriel ou gamma définie. Il m'a fallu quelques semaines, mais je fini par venir avec l'algorithme suivant qui calcule la CDF exactement (à savoir aucune approximation nécessaire). Python est fondamentalement aussi bon que pseudocode, droit?

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

échelles de performance avec x. Pour les petites valeurs de x, cette solution est d'un ordre de grandeur plus rapide que scipy.stats.binom.cdf, avec des performances similaires à environ x = 10.000.

Je ne vais pas aller dans une dérivation complète de cet algorithme, car stackoverflow ne supporte pas MathJax, mais l'idée maîtresse de c'est d'abord d'identifier l'équivalence suivante:

Pour tout k> 0, sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])

Ce que nous pouvons réécrire comme:

sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k

ou dans l'espace journal:

np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)

Parce que la CDF est une sommation de PMF, nous pouvons utiliser cette formulation pour calculer le coefficient binomial (le journal dont est b dans la fonction ci-dessus) pour PMF_ {x = i} du coefficient nous avons calculé pour PMF_ {x = i-1}. Cela signifie que nous pouvons faire tout à l'intérieur d'une seule boucle en utilisant les accumulateurs, et on n'a pas besoin de calculer les factorielles!

La raison pour la plupart des calculs sont effectués dans l'espace de journal est d'améliorer la stabilité numérique des termes polynomiaux, à savoir p^x et (1-p)^(1-x) ont le potentiel d'être très grande ou très petite, ce qui peut provoquer des erreurs de calcul.

EDIT: Est-ce un nouvel algorithme? J'ai farfouillé de façon intermittente depuis avant posté, et je me demande de plus en plus si je devais écrire cette place plus formelle et le soumettre à un journal.

Je pense que vous voulez évaluer les fonction bêta incomplète .

Il y a une belle mise en œuvre en utilisant une représentation en fraction continue « Recettes numériques En C », chapitre 6:. « Fonctions spéciales »

Je ne peux pas se porter garant totalement pour l'efficacité, mais Scipy a une Module pour cette

from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)

Un efficace et, plus important encore, l'algorithme stable numérique existe dans le domaine de Courbes Bézier utilisés dans la conception assistée par ordinateur. Il est appelé de l'algorithme de Casteljau utilisé pour évaluer les Bernstein Polynomials utilisé pour définir des courbes de Bézier.

Je crois que je ne suis autorisé un lien par réponse pour commencer par Wikipedia - Bernstein Polynomials

Remarquez la relation très étroite entre la distribution binomiale et le Polynôme de Bernstein. Cliquez ensuite par le lien sur l'algorithme de de Casteljau.

Disons que je sais que la probabilité de jeter une tête avec une pièce de monnaie particulière est P. Quelle est la probabilité de me jeter les temps T de la pièce et obtenir au moins têtes de S?

Set n = T
Set beta [i] = 0 pour i = 0, ... S - 1
Définir bêta [i] = 1 pour i = S, T ...
Définir t = p
Évaluer B (t) à l'aide de Casteljau

ou la plupart des chefs S?

Set n = T
Définir bêta [i] = 1 pour i = 0, ... S
Set beta [i] = 0 pour i = S + 1, ... T
Définir t = p
Évaluer B (t) à l'aide de Casteljau

Le code source ouvert existe probablement déjà. Les courbes NURBS (non uniformes rationnelles Courbes B-splines) sont une généralisation des courbes de Bézier et sont largement utilisés dans CAD. Essayez openNURBS (la licence est très libérale) ou à défaut Open CASCADE (une licence un peu moins libérale et opaque). Les deux boîtes à outils sont en C ++, cependant, IIRC, les liaisons .NET existent.

Si vous utilisez Python, pas besoin de ce code vous-même. Scipy pensé à vous:

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

A partir de la partie de votre question « obtenir au moins la tête de S » vous voulez que la fonction de distribution binomiale cummulative. Voir http://en.wikipedia.org/wiki/Binomial_distribution pour l'équation, qui est décrit comme étant en termes de « régularisé fonction bêta incomplète » (comme cela a déjà répondu). Si vous voulez juste pour calculer la réponse sans avoir à mettre en œuvre la solution complète vous, la GNU Scientific Library fournit la fonction:. Gsl_cdf_binomial_P et gsl_cdf_binomial_Q

Le DCDFLIB Projet a des fonctions C # (emballages autour du code C) pour évaluer un grand nombre fonctions CDF, y compris la distribution binomiale. Vous pouvez trouver l'original du code C et FORTRAN . Ce code est bien testé et précis.

Si vous voulez écrire votre propre code pour ne pas dépendre d'une bibliothèque externe, vous pouvez utiliser l'approximation normale au binomiale mentionné dans d'autres réponses. Voici quelques notes sur à quel point le rapprochement est dans diverses circonstances. Si vous choisissez cette voie et ont besoin d'un code pour calculer la CDF normale, voici code Python pour le faire . Il est seulement une douzaine de lignes de code et pourrait facilement être porté à une autre langue. Mais si vous voulez une grande précision et un code efficace, vous êtes mieux en utilisant le code tiers comme DCDFLIB. Plusieurs hommes-années se sont écoulées dans la production de cette bibliothèque.

celui-ci , utilisé dans les BPF. Une autre référence est cette .

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.

Licencié sous: CC-BY-SA avec attribution

Non affilié à StackOverflow