## code to accompany the heart Poisson example from lecture 4
## of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy adapted from Peter Hoff


## the prior (estimated from 10 hospitals)
a <- 16
b <- 15175

## observations for hospital A
ya <- 1
ea <- 66

## observations for hospital B
yb <- 4
eb <- 1767

## grid of lambda values
lambda <- seq(0.0001, 0.0025, length=1000)

## prior and posterior
prior <- dgamma(lambda, a, b)
pa <- dgamma(lambda, a+ya, b+ea)
pb <- dgamma(lambda, a+yb, b+eb)
plot(lambda, prior, type="l", lwd=2,
     ylab="density", xlab="lambda", bty="n")
lines(lambda, pa, col=2, lty=2, lwd=2)
lines(lambda, pb, col=3, lty=3, lwd=2)
legend("topright", c("prior", "post-A", "post-B"),
       col=1:3, lty=1:3, lwd=2)

## which hospital is worse
T <- 10000
mean(rgamma(T, a+ya, b+ea) < rgamma(T, a+yb, b+eb))