## Poisson
x <- seq(0.9, 2.5, len=199)
plot(x, dgamma(x, 217.5, 111), col=pal(2)[1], xlab=expression(theta), ylab="Posterior density", type="l", lwd=3)
lines(x, dgamma(x, 66.5, 44), col=pal(2)[2], lwd=3)
toplegend(c("Less than Bachelor's", "Bachelor's or higher"), ncol=2, col=pal(2), lwd=3)

## Poisson predictive
Y <- cbind(dnbinom(0:6, 0.5+217, 111/112),
           dnbinom(0:6, 0.5+66, 44/45))
X <- cbind(rep(-0.03,7), rep(0.03,7))+0:6
matplot(X, Y, col=pal(2), lwd=3, type="h", lty=1, xlab="Number of children", ylab="Probability")
toplegend(c("Less than Bachelor's", "Bachelor's or higher"), ncol=2, col=pal(2), lwd=3)

## Transplant study
y <- c(2, 3, 4, 4, 6, 7, 10, 12)
Pi <- rbeta(10000, 9, 3)
Lam <- rgamma(10000, 8, sum(y))
Omega <- Pi/Lam
mean(Omega)
mean(Omega > 2)
Success <- rbinom(10000, size=1, prob=Pi)
Surv <- rexp(10000, rate=Lam)
TotSurv <- Success*Surv
mean(TotSurv)
mean(TotSurv > 2)