Monday, September 10, 2007

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

2 comments:

slowmo said...

Hi Jim,
For the posterior predictive distribution you have:

ys=rpois(1000,500*lambda)

where lambda was previously calculated as:

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

I'm a bit unsure about this. Isn't the rpois call only using the first element of the lambda vector ? I would have thought to do something like this:
ys <- numeric(1000)
for ( i in 1:1000 ) ys[i] <- rpois(1, sample(lambda,1))

Am I missing something ?

slowmo said...

PS. sorry, obviously that should have been rpois(1, 500*sample(lambda, 1))