Poisson confidence interval - why limit it to < 100 counts?
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
In the statistics toolbox, we use the function poissfit to calculate the parameter lambda (=mean =variance) of poisson-distributed data, together with a confidence interval. The confidence interval calculation is done by statpoisci, which uses an inverse chi^2 method. However, for data with more than 99 counts, statpoisci gives up and approximates the Poisson distribution with a normal distribution. The relevant code is:
% Chi-square exact method
lb(k) = chi2inv(alpha/2, 2*lsum(k))/2;
ub(k) = chi2inv(1-alpha/2, 2*(lsum(k)+1))/2;
where lsum is the total number of events from the Poisson process, k selects elements with <100 events, alpha is the confidence level (e.g. 0.95) and chi2inv is a built-in function. But for elements with >=100 events, it inverts k and uses:
% Normal approximation
lb(k) = norminv(alpha/2,lsum(k),sqrt(lsum(k)));
ub(k) = norminv(1 - alpha/2,lsum(k),sqrt(lsum(k)));
The normal approximation isn't bad for the lower bound, but the upper bound is noticeably underestimated. The strange thing is, if I remove the "<100" test, the chi-square exact method works up to about lsum <~ 5e14. Above that level, gaminv starts to fail (you get a convergence warning).
My questions are: why does statpoisci use an event limit that seems to be much lower than necessary? Why is a bog-standard normal approximation used above that limit, given that better approximations are available that still preserve the asymmetry of Poisson confidence limits (see e.g. Wikipedia or STSI, the latter of which gives a near-perfect upper interval for 1-sigma confidence)?
Thanks,
David.
0 commentaires
Réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!