Confidence Intervals for Poisson
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