Częstotliwość przedziału Poissona z numpy

Question

Próbuję wstawić paski błędu Poissona na histogramie, który robię z matplotlib, ale nie mogę znaleźć funkcji numpy, która da mi 95% przedział ufności przy założeniu poissonian dane. Idealnie rozwiązanie nie zależy od scipy, ale wszystko będzie działać. Czy taka funkcja istnieje? Znalazłem dużo o bootstrapowaniu, ale w moim przypadku wydaje się to trochę przesadzone.

Answer 1

Korzystanie scipy.stats.poisson oraz metodę interval:

>>> scipy.stats.poisson.interval(0.95, [10, 20, 30]) 
(array([ 4., 12., 20.]), array([ 17., 29., 41.])) 

Choć to tylko sprawia, że ​​ograniczony sens, aby obliczyć rozkład Poissona o wartości nie jest liczbą całkowitą, to dokładne przedziały ufności wnioskowane przez OP można obliczyć można to zrobić w następujący sposób:

>>> data = np.array([10, 20, 30]) 
>>> scipy.stats.poisson.interval(0.95, data) 
(array([ 4., 12., 20.]), array([ 17., 29., 41.])) 
>>> np.array(scipy.stats.chi2.interval(.95, 2 * data))/2 - 1 
array([[ 3.7953887 , 11.21651959, 19.24087402], 
     [ 16.08480345, 28.67085357, 40.64883744]]) 

---

możliwe jest również do korzystania z ppf metoda:

>>> data = np.array([10, 20, 30]) 
>>> scipy.stats.poisson.ppf([0.025, 0.975], data[:, None]) 
array([[ 4., 17.], 
     [ 12., 29.], 
     [ 20., 41.]]) 

Ale ponieważ dystrybucja jest dyskretne wartości powrotne będą liczbami całkowitymi, a przedział ufności nie obejmie 95% dokładnie:

>>> scipy.stats.poisson.ppf([0.025, 0.975], 10) 
array([ 4., 17.]) 
>>> scipy.stats.poisson.cdf([4, 17], 10) 
array([ 0.02925269, 0.98572239]) 

Answer 2

skończyło się pisać własne funkcje na podstawie some properties I found on Wikipedia.

def poisson_interval(k, alpha=0.05): 
    """ 
    uses chisquared info to get the poisson interval. Uses scipy.stats 
    (imports in function). 
    """ 
    from scipy.stats import chi2 
    a = alpha 
    low, high = (chi2.ppf(a/2, 2*k)/2, chi2.ppf(1-a/2, 2*k + 2)/2) 
    if k == 0: 
     low = 0.0 
    return low, high 

Powoduje zwrócenie ciągłych (a nie dyskretnych) ograniczeń, co jest bardziej standardowe w moim polu.

Answer 3

Ten problem pojawia się wiele w astronomii (moim polu!) I ten papier jest iść do odniesienia dla tych przedziałów ufności: Gehrels 1980

Ma dużo matematyki w nim dla dowolnego przedziału ufności z Statystyki Poissona, ale dla dwustronnego 95% przedziału ufności (odpowiadającego 2-sigma gaussowskiego przedziału ufności, lub S = 2 w kontekście tego artykułu) niektórych prostych analitycznych formuł dla górnych i dolnych granic ufności dla zdarzeń N są mierzone są

upper = N + 2. * np.sqrt(N + 1) + 4./3. 
lower = N * (1. - 1./(9. * N) - 2./(3. * np.sqrt(N))) ** 3. 

gdzie umieściłem je w formacie Python dla Ciebie alr eady. Wszystko czego potrzebujesz to numpy lub twój drugi ulubiony moduł pierwiastkowy. Należy pamiętać, że dają one górną i dolną granicę dla zdarzeń - a nie wartości +/-. Po prostu odejmij N od obu, aby je uzyskać.

Prosimy o zapoznanie się z papierem w celu uzyskania dokładności tych wzorów dla potrzebnego przedziału ufności, ale powinny one być więcej niż wystarczająco dokładne dla większości praktycznych zastosowań.