# R program to fit Poisson distribution to frequency data and
# simulate distribution of ML estimate (and obtain confidence interval)
# using parametric bootstrapping.
#
# Enter frequency data in the following two vectors:
xx=c(0,1,2,3,4,5)
ff=c(213,128,37,18,3,1)  #  Rice, third ed., problem 8.3 data set 1
nsim=1000
alpha=.05

n=sum(ff)
lambda.ml=sum(xx*ff)/n


lambda.star=numeric(nsim)
for (i in 1:nsim) {
  lambda.star[i]=mean(rpois(n,lambda.ml))
}
hist(lambda.star,main="",xlab="lambda values")
expected=n*exp(-lambda.ml+xx*log(lambda.ml)-lfactorial(xx))
observed=ff
cbind(xx,observed,expected)

lambda.star=sort(lambda.star)
lambda.lo=lambda.star[floor((alpha/2)*nsim)]
lambda.hi=lambda.star[ceiling((1-alpha/2)*nsim)]
lambda.ci=c(lambda.lo,lambda.hi)
lambda.ml
lambda.ci