Poisson confidence interval - why limit it to < 100 counts?
2 views (last 30 days)
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)?