# Poisson-gamma S-plus code for Heart Valves problem, by B. Carlin
#  (last update:  7/10/02)

M <- c(.038, .0166, .0166); nprior <- length(M)
sd <- c(2*.038, .0166, .0166/2);  V <- sd^2
alpha <- M^2/V;  beta <- V/M

y <- c(1,3,20);  ndata <- length(y)
t <- 200;  OPC <- .038

a <- matrix(0,nprior,ndata);  b <- a

thet <- seq(from=.001,to=.2,length=200)
par(mfcol=c(nprior,ndata),oma=c(1,0,4,0))

#  HERE ARE DATA AND PRIOR LOOPS:

for(j in 1:ndata){
for(i in 1:nprior){

a[i,j] <- y[j]+alpha[i];  b[i,j] <- 1/(t + 1/beta[i])
  post <- dgamma(thet/b[i,j],a[i,j])/b[i,j]
  prior <- dgamma(thet/beta[i],alpha[i])/beta[i]
  top <- max(post,prior)
  plot(thet,post,type="l",xlab="",ylab="",ylim=c(0,top))
  lines(thet,prior,lty=2)
  legend(0.1,top,c("posterior","prior"),lty=1:2)
  title(sub=paste("M, sd =",round(M[i],3),round(sd[i],3),
    "; Y =",y[j]),cex=1.01)

  prob <- pgamma(2*OPC/b[i,j], a[i,j])
  cat(paste("\nFor prior ",i,", prob that theta < 2 OPC) = ",round(prob,3)))
  text(0.13,top/2,paste("P(theta < 2 OPC|y) = ",round(prob,3)),cex=1.05)
}
}

mtext(side=3,line=1,cex=1.4,outer=TRUE,
  "Priors and posteriors, Heart Valves Study A, Poisson-gamma model")
mtext(side=3,line=-2,cex=1.4,outer=TRUE,
  paste("for various prior (M, sd) and data (y) values; T =",t,
    ", 2 OPC =",2*OPC))

cat("\nThat's all!\n")