x=c(31,29,19,18,31,28,
34,27,34,30,16,18,
26,27,27,18,24,22,
28,24,21,17,24)
lambda.hat=mean(x)
print(lambda.hat)
## [1] 24.91304
par(mfcol=c(1,1))
hist(x, xlim=c(0,1.5*max(x)), probability=TRUE)

hist(x,nclass=15, xlim=c(0,1.5*max(x)), probability=TRUE)

hist(x,nclass=10, xlim=c(0,1.5*max(x)), probability=TRUE)
x.grid=seq(0,1.5*max(x),1)
x.probs=dpois(x.grid, lambda=lambda.hat)
lines(x.grid,x.probs, type="h", col='blue',lwd=4)

lambda.hat.sterror=sqrt(lambda.hat/length(x))
print(lambda.hat.sterror)
## [1] 1.040757
lambda.CI.Limits=lambda.hat + c(-1,1)*qnorm(.95)*lambda.hat.sterror
print(lambda.CI.Limits)
## [1] 23.20115 26.62494
qnorm(.95)
## [1] 1.644854
qnorm(.05)
## [1] -1.644854
nu0= 15/25
alpha0=nu0*15
lambda.grid=seq(0.1,30,.1)
lambda.grid.priorpmf=dgamma(lambda.grid, shape=alpha0, rate=nu0)
priorpmf.grid =dgamma(lambda.grid, shape=alpha0, rate=nu0)
plot(lambda.grid, priorpmf.grid,
col='red', type="l",
main="Prior Density")

likelihood.grid=dpois(sum(x),lambda=length(x)*lambda.grid)
plot(lambda.grid,likelihood.grid,
main="Likelihood of lambda", col='green',type="l")

alpha1= alpha0 + sum(x)
nu1 = nu0 + length(x)
posteriorpmf.grid=dgamma(lambda.grid,shape=alpha1, rate=nu1)
plot(lambda.grid, posteriorpmf.grid,
main="Posterior Density", col='blue', type="l")

par(mfcol=c(1,1))
plot(lambda.grid, posteriorpmf.grid,ylab="density",
main="Densities \n Gamma Prior/Posterior (red/blue)\n",
type="n")
lines(lambda.grid, priorpmf.grid, col='red')
lines(lambda.grid, posteriorpmf.grid, col='blue')
lines(lambda.grid, (1+0*priorpmf.grid)*.1/30, col='orange')
posteriorpmf.uniformprior.grid=dgamma(lambda.grid,shape=1+ sum(x), rate= length(x))
lines(lambda.grid, posteriorpmf.uniformprior.grid, col='green')
title(main="\n\nUniform Prior/Posterior (orange/green)")

posterior.mean=alpha1/nu1
posterior.stdev=sqrt(alpha1/(nu1*nu1))
posterior.mode= (alpha1-1)/nu1
posterior.llimit=qgamma(.05, shape=alpha1, rate=nu1)
posterior.ulimit=qgamma(.95, shape=alpha1, rate=nu1)
bayes1.attributes=data.frame(
mode=posterior.mode,
mean=posterior.mean,
stdev=posterior.stdev,
ulimit=posterior.ulimit,
llimit=posterior.llimit)
alpha1=sum(x)
nu1=length(x)
posterior.mean=alpha1/nu1
posterior.stdev=sqrt(alpha1/(nu1*nu1))
posterior.mode= (alpha1-1)/nu1
posterior.llimit=qgamma(.05, shape=alpha1, rate=nu1)
posterior.ulimit=qgamma(.95, shape=alpha1, rate=nu1)
bayes2.attributes=data.frame(
mode=posterior.mode,
mean=posterior.mean,
stdev=posterior.stdev,
ulimit=posterior.ulimit,
llimit=posterior.llimit)
mle.attributes=data.frame(
mode=lambda.hat,
mean=NA,
stdev=sqrt(lambda.hat/length(x)),
ulimit=lambda.hat + qnorm(.95)*sqrt(lambda.hat/length(x)),
llimit=lambda.hat + qnorm(.05)*sqrt(lambda.hat/length(x)))
estimates.table=data.frame(cbind(Bayes1=t(bayes1.attributes), Bayes2=t(bayes2.attributes),
MLE=t(mle.attributes))
)
dimnames(estimates.table)[[2]]<-c("Bayes 1", "Bayes 2", "Maximum Likelihood")
print(estimates.table)
## Bayes 1 Bayes 2 Maximum Likelihood
## mode 24.618644 24.869565 24.913043
## mean 24.661017 24.913043 NA
## stdev 1.022232 1.040757 1.040757
## ulimit 26.366182 26.649296 26.624937
## llimit 23.004027 23.226222 23.201150