Comment trouver un intervalle de confiance pour le nombre total d'événements

J'ai un détecteur qui détectera un événement avec une certaine probabilité p . Si le détecteur indique qu'un événement s'est produit, c'est toujours le cas, il n'y a donc pas de faux positifs. Après l'avoir exécuté pendant un certain temps, je reçois k événements détectés. Je voudrais calculer avec certitude 95% du nombre total d'événements qui se sont produits, détectés ou non, avec une certaine confiance.

Ainsi, par exemple, supposons que 13 événements soient détectés. J'aimerais pouvoir calculer qu'il y a eu entre 13 et 19 événements avec une confiance de 95% basée sur p .

Voici ce que j'ai essayé jusqu'à présent:

La probabilité de détecter k événements s'il y avait n total est:

binomial(n, k) * p^k * (1 - p)^(n - k)

La somme de cela sur n de k à l'infini est:

1/p

Ce qui signifie que la probabilité qu'il y ait n événements au total est:

f(n) = binomial(n, k) * p^(k + 1) * (1 - p)^(n - k)

Donc, si je veux être sûr à 95%, je devrais trouver la première somme partielle f(k) + f(k+1) + f(k+2) ... + f(k+m)qui est au moins 0,95 et la réponse est [k, k+m]. Est-ce la bonne approche? Existe-t-il également une formule fermée pour la réponse?

Je choisirais d'utiliser la distribution binomiale négative , qui renvoie la probabilité qu'il y aura X échecs avant le k_ième succès, lorsque la probabilité constante d'un succès est p.

Utiliser un exemple

k=17 # number of successes
p=.6 # constant probability of success

la moyenne et la sd des défaillances sont données par

mean.X <- k*(1-p)/p
sd.X <- sqrt(k*(1-p)/p^2) 

La distribution des échecs X, aura approximativement cette forme

plot(dnbinom(0:(mean.X + 3 * sd.X),k,p),type='l')

Ainsi, le nombre d'échecs sera (avec une confiance de 95%) approximativement entre

qnbinom(.025,k,p)
[1] 4

et

qnbinom(.975,k,p)
[1] 21

Donc, votre inertie serait [k + qnbinom (.025, k, p), k + qnbinom (.975, k, p)] (en utilisant les nombres de l'exemple [21,38])

En supposant que vous vouliez choisir une distribution pour n, p (n), vous pouvez appliquer la loi de Bayes.

Vous savez que la probabilité que k événements se produisent étant donné que n se sont réellement produits est régie par une distribution binomiale

$p(k|n) = {n \choose k} p^k (1-p)^{(n-k)}$

Ce que vous voulez vraiment savoir, c'est la probabilité que n événements se soient réellement produits, étant donné que vous avez observé k. Par Bayes était:

$p(n|k) = \frac{p(k|n)p(n)}{p(k)}$

En appliquant le théorème de la probabilité totale, nous pouvons écrire:

$p(n|k) = \frac{p(k|n)p(n)}{\sum_{n'} p(k|n')p(n')}$

Donc, sans plus d'informations sur la distribution de vous ne pouvez pas vraiment aller plus loin.$p(n)$

Cependant, si vous souhaitez choisir une distribution pour pour laquelle il existe une valeur supérieure à laquelle , ou suffisamment proche de zéro, vous pouvez faire un peu mieux. Par exemple, supposons que la distribution de soit uniforme dans la plage . ce cas:$p(n)$$n$$p(n) = 0$$n$$[0,n_{max}]$

$p(n) = \frac{1}{n_{max}}$

La formulation bayésienne se simplifie pour:

$p(n|k) = \frac{p(k|n)}{\sum_{n'} p(k|n')}$

En ce qui concerne la dernière partie du problème, je conviens que la meilleure approche consiste à effectuer une sommation cumulative sur , à générer la fonction de distribution de probabilité cumulée et à itérer jusqu'à ce que la limite de 0,95 soit atteinte.$p(n|k)$

Étant donné que cette question a migré de SO, un exemple de code de jouet en python est joint ci-dessous

import numpy.random

p = 0.8
nmax = 200

def factorial(n):
    if n == 0:
        return 1
    return reduce( lambda a,b : a*b, xrange(1,n+1), 1 )

def ncr(n,r):
    return factorial(n) / (factorial(r) * factorial(n-r))

def binomProbability(n, k, p):
    p1 = ncr(n,k)
    p2 = p**k
    p3 = (1-p)**(n-k)
    return p1*p2*p3

def posterior( n, k, p ):
    def p_k_given_n( n, k ):
        return binomProbability(n, k, p)
    def p_n( n ):
        return 1./nmax
    def p_k( k ):
        return sum( [ p_n(nd)*p_k_given_n(nd,k) for nd in range(k,nmax) ] )
    return (p_k_given_n(n,k) * p_n(n)) / p_k(k)


observed_k   = 80
p_n_given_k  = [ posterior( n, observed_k, p ) for n in range(0,nmax) ]
cp_n_given_k = numpy.cumsum(p_n_given_k)
for n in xrange(0,nmax):
    print n, p_n_given_k[n], cp_n_given_k[n]

Si vous mesurez événements et savez que votre efficacité de détection est vous pouvez automatiquement corriger votre résultat mesuré jusqu'au nombre "vrai" .$k$$p$$k_\mathrm{true} = k/p$

Votre question est alors de trouver la plage de où 95% des observations tomberont. Vous pouvez utiliser la méthode Feldman-Cousins pour estimer cet intervalle. Si vous avez accès à ROOT, il y a une classe pour faire ce calcul pour vous.$k_\mathrm{true}$

Vous calculeriez les limites supérieure et inférieure avec Feldman-Cousins ​​à partir du nombre non corrigé d'événements , puis les augmenteriez jusqu'à 100% avec . De cette façon, le nombre réel de mesures détermine votre incertitude, pas un certain nombre à l'échelle qui n'a pas été mesuré.$k$$1/p$

{
gSystem->Load("libPhysics");

const double lvl = 0.95;
TFeldmanCousins f(lvl);

const double p = 0.95;
const double k = 13;
const double k_true = k/p;

const double k_bg = 0;

const double upper = f.CalculateUperLimit(k, k_bg) / p;
const double lower = f.GetLowerLimit() / p;

std::cout << "["
  lower <<"..."<<
  k_true <<"..."<<
  upper <<
  "]" << std::endl;
}

Je pense que vous avez mal compris l'objectif des intervalles de confiance. Les intervalles de confiance vous permettent d'évaluer où se trouve la vraie valeur du paramètre. Donc, dans votre cas, vous pouvez construire un intervalle de confiance pour . Il n'est pas logique de construire un intervalle pour les données.$p$

Cela dit, une fois que vous avez une estimation de vous pouvez calculer la probabilité que vous observerez différentes réalisations telles que 14, 15, etc. en utilisant le binôme pdf.$p$