## 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

Subscribe to:
Post Comments (Atom)

## 2 comments:

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 ?

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

Post a Comment