# 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")