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