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

Disons que je sais que la probabilité d'un "succès" de P. je lance le test de N fois, et je vois réussites. Le test est semblable à jeter un inégalement pondérée de la pièce de monnaie (peut-être les chefs est un succès, la queue est un échec).

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

Ainsi, par exemple, si P est de 0,3, ce qui N est de 100, et j'ai 20 réussites, je suis à la recherche de la probabilité d'obtenir 20 ou moins succès.

Si, sur l'autre n'avait, P est de 0,3, ce qui N est de 100, et j'obtiens 40 succès, je suis à la recherche de la probabilité d'obtenir 40 nos plus de succès.

Je suis conscient que ce problème est lié à la découverte de la région en vertu d'une loi binomiale courbe, cependant:

Mes maths-fu n'est pas à la tâche de traduire cette connaissance en code efficace
Bien que je comprenne un binôme courbe devrait donner un résultat exact, j'ai l'impression qu'il serait intrinsèquement inefficace. Une méthode rapide pour calculer un résultat approximatif suffirait.

Je tiens à souligner que ce calcul doit être rapide, et devrait idéalement être déterminés avec le standard de 64 ou 128 bits à virgule flottante calcul.

Je suis à la recherche d'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 toutes les réponses employer des pseudo-code ou code.

source d'informationauteur sanity

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

Normal Devis, bon pour n grand

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

6

Je suis sur un projet où nous en avions besoin pour être en mesure de calculer le binôme CDF dans un environnement qui n'ont pas une factorielle ou gamma fonction définie. Il m'a fallu quelques semaines, mais j'ai fini par venir avec l'algorithme suivant qui calcule la CDF exactement (c'est à dire pas de rapprochement nécessaire). Python est fondamentalement aussi bon que de pseudo, non?
```
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
```
De la Performance des gammes avec x. Pour de petites valeurs de x, cette solution est environ un ordre de grandeur plus rapide que scipy.stats.binom.cdfavec des performances similaires à autour de x=10 000 habitants.

Je n'entrerai pas dans une dérivation de cet algorithme, car stackoverflow ne prend pas en charge MathJax, mais la poussée est d'abord 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)])
Qui nous pouvons réécrire comme:
- sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k
ou dans l'espace du 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 le CDF est une somme de Pmf, on peut utiliser cette formule pour calculer le coefficient binomial (le journal de ce qui est b dans la fonction ci-dessus) pour PMF_{x=i} et le coefficient de nous avons calculé pour PMF_{x=i-1}. Cela signifie que nous pouvons tout faire à l'intérieur d'une seule boucle à l'aide d'accumulateurs, et nous n'avons pas besoin de calculer les factorielles!

La raison pour laquelle la plupart des calculs sont faits dans le journal de l'espace est d'améliorer la stabilité numérique du polynôme, c'est à dire p^x et (1-p)^(1-x) ont le potentiel d'être très grand ou très petit, ce qui peut provoquer des erreurs de calcul.
4

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

Il y a une belle mise en œuvre à l'aide d'une fraction continue de la représentation en "Numerical Recipes In C", chapitre 6: "les Fonctions Spéciales'.

Je ne peux pas totalement se porter garant de l'efficacité, mais Scipy a un pour ce module

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

2

Efficace et, surtout, numérique algorithme stable existe dans le domaine de des Courbes de Bézier utilisé dans la Conception Assistée par Ordinateur. Il est appelé de l'algorithme de Casteljau utilisé pour évaluer la Polynômes de Bernstein utilisé pour définir des Courbes de Bézier.

Je crois que je suis seulement permis à un lien par réponse donc de commencer avec Wikipédia - Bernstein Pour Les Polynômes

Avis de la relation très étroite entre la Distribution Binomiale et les Polynômes de Bernstein. Ensuite, cliquez sur le lien de l'algorithme de Casteljau.

Permet de dire que je sais que la probabilité de jeter une tête avec une pièce de monnaie est P.
Quelle est la probabilité de me lancer
la pièce en T de fois et d'obtenir au moins
S têtes?
- Ensemble n = T
- Ensemble bêta[i] = 0 pour i = 0, ... S - 1
- Ensemble bêta[i] = 1 pour i = S, ... T
- Ensemble t = p
- Évaluer B(t) à l'aide de Casteljau
ou à la plupart des têtes?
- Ensemble n = T
- Ensemble bêta[i] = 1 pour i = 0, ... S
- Ensemble bêta[i] = 0 pour i = S + 1, ... T
- Ensemble t = p
- Évaluer B(t) à l'aide de Casteljau
Le code Open source existe sans doute déjà. Les Courbes NURBS (Non-Uniform Rational B-spline) sont une généralisation des Courbes de Bézier et sont largement utilisés en CAO. Essayez openNurbs (la licence est très libéral) ou, à défaut, Ouvrez la CASCADE (un peu moins libéral et opaque de licence). Les deux kits sont en C++, bien que, autant que je me souvienne, .NET liaisons existent.

Si vous utilisez Python, pas besoin de code vous-même. Scipy que vous êtes couvert:

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

1

De la partie de votre question "obtenir au moins S chefs" vous voulez le cumulée de la distribution binomiale fonction. Voir http://en.wikipedia.org/wiki/Binomial_distribution de l'équation, qui est décrite en termes de "régularisés en fonction bêta incomplète" (comme déjà répondu). Si vous voulez juste pour calculer la réponse sans avoir à mettre en œuvre l'ensemble de la solution vous-même, la GNU Scientific Library fournit la fonction: gsl_cdf_binomial_P et gsl_cdf_binomial_Q.
1

La DCDFLIB Projet a C# fonctions (wrappers autour de code C) à évaluer de nombreux CDF fonctions, y compris la distribution binomiale. Vous pouvez trouver l'original C et FORTRAN code ici. Ce code est bien testé et précis.

Si vous souhaitez écrire votre propre code pour éviter d'être dépendant d'une bibliothèque externe, vous pouvez utiliser l'approximation normale de la loi binomiale mentionné dans d'autres réponses. Voici quelques notes sur la qualité de l'approximation est dans des circonstances différentes. Si vous aller dans cette voie et ont besoin de code pour calculer la normale CDF, voici Le code Python pour le faire. Il ne s'agit que d'une dizaine de lignes de code et pourrait facilement être porté à toute autre langue. Mais si vous voulez une grande exactitude et efficacité du code, il est préférable d'utiliser un code tiers comme DCDFLIB. Plusieurs années de travail consacrés à la production de cette bibliothèque.
0

Essayer cette uneutilisé dans des BPF. Une autre référence est cette.

Vous devez vous connecter pour publier un commentaire.