## Inference for a Poisson Rate

In 2006, Ryan Howard had 58 home runs in 581 at bats. Suppose the number of home runs y is distributed Poisson with mean ab lambda, where lambda is the true home run rate. If we assume the usual noninformative prior, the posterior for lambda is Gamma with shape alpha=y and rate ab.

To learn about lambda, we simulate 1000 draws from the Gamma posterior. We construct a density estimate of the simulated draws and then find a 90% interval estimate by finding the 5th and 95th percentiles of the draws.

> y=58

> t=581

> lambda=rgamma(1000,shape=y,rate=t)

> plot(density(lambda),xlab=”LAMBDA”,main=”POSTERIOR”)

> quantile(lambda,c(.05,.95))

5% 95%

0.07897515 0.12170184

We see that a 90% interval estimate for lambda is (.079, .122).

Next, suppose that Howard has 550 at-bats in 2007 — how many home runs will he hit? We learn about the future number of home runs ys by use of the posterior predictive distribution. In the R output, we (1) simulate 1000 draws from the posterior predictive distribution, (2) tabulate the values by use of the table function, (3) construct a graph of this distribution, and (4) construct a 90% interval estimate for ys by use of the discint function (from the LearnBayes package).

> ys=rpois(1000,500*lambda)

> T=table(ys)

> plot(T,main=”POSTERIOR PREDICTIVE”)

>

> ys=as.integer(names(T))

> freq=as.integer(T)

> freq=freq/sum(freq)

> dist=cbind(ys,freq)

> discint(dist,.9)

$prob

[1] 0.905

$set

[1] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[26] 61 63 64 65 66 67