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

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 |

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

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