## Prediction problem

In class, I described a simple Bayesian prediction problem. One observes the number of fire calls in a small California city for 30 days — the total number of fire calls in this period is 158. An administrator is interested in predicting the total number of calls in the next seven days.

We assume that y, the number of calls in a day, is Poisson with mean . Suppose our prior beliefs about are described by a normal density with mean 10 and standard deviation 2.

I first write a simple function compute.post that computes and plots the posterior (using the likelihood times prior recipe) on a grid between lo and hi.

`compute.post=function(lo, hi,...){`

lambda = seq(lo, hi, length=100)

prior = dnorm(lambda, 10, 2)

like = lambda^158 * exp(-30*lambda)

post = prior * like

plot(lambda, post,...)

list(L = lambda, P = post)

}

With some trial and error, the interval (3.3, 7.5) seems to bracket the posterior well.

`S=compute.post(3.3, 7.5, type="l", col="red")`

Next I take a simulated sample from the posterior by using the sample function.

`sim.lambda = sample(S$L, 1000, replace=TRUE, prob=S$P)`

Last, I simulate a sample from the posterior predictive density of z, the number of fire calls in seven days. Since [z | ] is Poisson(7 , we simulate z by (1) simulating from the posterior and (2) simulating z from a Poisson(7 ) distribution. This is easy to do by using the vector properties of the random simulation functions.

`sim.z = rpois(1000, 7*sim.lambda)`

Here is a graph of the predictive distribution (by tabulating and plotting the simulated values of z).