Jak znaleźć przedział ufności dla całkowitej liczby zdarzeń

Mam detektor, który wykryje zdarzenie z pewnym prawdopodobieństwem p . Jeśli wykrywacz powie, że zdarzenie miało miejsce, zawsze tak jest, więc nie ma fałszywych trafień. Po uruchomieniu przez pewien czas wykryto k zdarzeń. Chciałbym obliczyć całkowitą liczbę zdarzeń, które miały miejsce, zostały wykryte lub w inny sposób, z pewną pewnością, powiedzmy 95%.

Załóżmy na przykład, że wykryto 13 zdarzeń. Chciałbym móc obliczyć, że miało miejsce od 13 do 19 zdarzeń z 95% pewnością na podstawie p .

Oto, co próbowałem do tej pory:

Prawdopodobieństwo wykrycia k zdarzeń, jeśli n było łącznie, wynosi:

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

Suma tego ponad n od k do nieskończoności wynosi:

1/p

Co oznacza, że ​​prawdopodobieństwo wystąpienia n zdarzeń ogółem wynosi:

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

Więc jeśli chcę mieć 95% pewności, powinienem znaleźć pierwszą sumę częściową, f(k) + f(k+1) + f(k+2) ... + f(k+m)która wynosi co najmniej 0,95, a odpowiedź brzmi [k, k+m]. Czy to jest właściwe podejście? Czy istnieje również zamknięta formuła odpowiedzi?

Wybrałbym użycie ujemnego rozkładu dwumianowego , który zwraca prawdopodobieństwo, że wystąpią awarie X przed k_tym sukcesem, gdy stałe prawdopodobieństwo sukcesu wynosi p.

Na przykładzie

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

średnią i sd dla awarii podano przez

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

Rozkład awarii X będzie miał w przybliżeniu ten kształt

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

Tak więc liczba awarii będzie (z 95% pewnością) w przybliżeniu pomiędzy

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

i

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

Więc twój inerval to [k + qnbinom (.025, k, p), k + qnbinom (.975, k, p)] (używając liczb z przykładu [21,38])

Zakładając, że chcesz wybrać rozkład dla n, p (n), możesz zastosować prawo Bayesa.

Wiesz, że prawdopodobieństwo wystąpienia k zdarzeń, biorąc pod uwagę fakt, że n rzeczywiście miało miejsce, jest regulowane rozkładem dwumianowym

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

Rzeczą, którą naprawdę chcesz wiedzieć, jest prawdopodobieństwo wystąpienia n zdarzeń, biorąc pod uwagę, że zaobserwowałeś k. By Bayes leżał:

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

Stosując twierdzenie o całkowitym prawdopodobieństwie, możemy napisać:

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

Bez dalszych informacji na temat rozkładu nie można tak naprawdę pójść dalej.$p(n)$

Jeśli jednak chcesz wybrać rozkład dla dla którego istnieje wartość większa niż która lub wystarczająco bliska zeru, możesz zrobić trochę lepiej. Załóżmy na przykład, że rozkład jest jednorodny w zakresie . ta sprawa:$p(n)$$n$$p(n) = 0$$n$$[0,n_{max}]$

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

Bayesowska formuła upraszcza:

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

Jeśli chodzi o ostatnią część problemu, zgadzam się, że najlepszym podejściem jest wykonanie sumowania skumulowanego nad , wygenerowanie funkcji skumulowanego rozkładu prawdopodobieństwa i iteracja do osiągnięcia limitu 0,95.$p(n|k)$

Biorąc pod uwagę, że to pytanie migrowało z SO, przykładowy kod zabawki w pythonie znajduje się poniżej

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]

Jeśli mierzysz zdarzeń i wiesz, że skuteczność wykrywania wynosi , możesz automatycznie skorygować zmierzony wynik do „prawdziwej” liczby .$k$$p$$k_\mathrm{true} = k/p$

Twoje pytanie dotyczy zatem znalezienia zakresu którym spadnie 95% obserwacji. Do oszacowania tego odstępu można użyć metody Feldmana-Cousinsa . Jeśli masz dostęp do ROOT, istnieje klasa, która wykona dla Ciebie te obliczenia.$k_\mathrm{true}$

Górną i dolną granicę obliczysz z Feldman-Cousins ​​na podstawie nieskorygowanej liczby zdarzeń a następnie skalujesz do 100% za pomocą . W ten sposób rzeczywista liczba pomiarów określa twoją niepewność, a nie jakaś skalowana liczba, która nie została zmierzona.$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;
}

Myślę, że źle zrozumiałeś cel przedziałów ufności. Przedziały ufności pozwalają ocenić, gdzie znajduje się prawdziwa wartość parametru. Tak więc w twoim przypadku możesz zbudować przedział ufności dla . Nie ma sensu konstruować przedziału dla danych.$p$

Powiedziawszy to, po oszacowaniu możesz obliczyć prawdopodobieństwo, że zaobserwujesz różne realizacje, takie jak 14, 15 itd., Używając dwumianowego pdf.$p$