[R] Confidence Intervals for Poisson

[David Winsemius]

On Jan 15, 2009, at 1:10 AM, John Kerpel wrote:

> Hi folks!
> I run the following code to get a CI for a Poisson with lambda=12.73
>
> library(MASS)
>
> set.seed(125)
>
> x <- rpois(100,12.73)
>
> lambda_hat<-fitdistr(x, dpois, list(lambda=12))$estimate
>
> #Confidence Intervals - Normal Approx.
>
> alpha<-c(.05,.025,.01)
>
> for(n in 1:length(alpha)) {
>
> LowerCI<-mean(x)-(qnorm(1-alpha[n]/2, mean = 0, sd =
> 1)*sqrt(var(x)/length(x)))
>
> UpperCI<-mean(x)+(qnorm(1-alpha[n]/2, mean = 0, sd =
> 1)*sqrt(var(x)/length(x)))
>
> cat("For
> Alpha
> =
> ",alpha
> [n
> ],"LowerCI
> =",LowerCI,"<","Lambda=",mean(x),"<","UpperCI=",UpperCI,"\n")
>
>
> }
>
>
> When I do something like:
>
> qpois(.975, 12.73, lower.tail = TRUE, log.p = FALSE)
> [1] 20
>> qpois(.025, 12.73, lower.tail = TRUE, log.p = FALSE)
> [1] 6
>
> I get quite a different result.  Is this the difference between the  
> normal
> approx and an (almost) exact Poisson CI?

Yes and no.

Using Byar's approximation, which is reasonably accurate at this  
expected value, I get 6.2 and 20.9 so R's qpois seems pretty sensible.

Your results don't look like a proper creation of a Normal approx,  
however. Weren't you worried that your code might not be performing as  
desired when the upper CL for your alpha= 0.05, and 0.01 results were  
only different by 0.3?

I would have thought a (much more simple) Normal approximation for the  
Poisson 0.05 CL around an expected of E might be
    E +/- 1.96* E^(.5).  So  12 +/- 2* 3.4 (5.2, 18.8)  might be an  
eyeball estimate.

 > 12 + 1.96*sqrt(12)
[1] 18.78964
 > 12 - 1.96*sqrt(12)
[1] 5.210361


-- 
David Winsemius