# Beta-binomial design code, by Sparks and Carlin, 10/10/01
#   (last update:  2/12/02)
# Here the priors are all of the "weighted" form, 
#    ie Beta(110w,7w), 0 < w < 1

# NOTE:  New data is generated from the assumed prior!
#   (We *could* generate it from some other distribution...)

set.seed(99)
w <- c(1,.5,.1)
N <- c(20,50)
C <- c(.85,.86,.87,.88,.89,.90)
x <- seq(from=.8,to=.998,length=50);  Nrep <- 100

par(mfcol=c((length(N)+1),length(w)),oma=c(1,0,4,0))
pct <- matrix(0,length(C),length(w)*length(N)+1);  pct[,1] <- C

#  HERE IS CHRONICLE "A" DATA:
no.ae <- 110;  ae <- 7

#  NOW START MAIN SIMULATION LOOPS:

for(j in 1:length(w)){

  alpha <- no.ae*w[j]
  beta <- ae*w[j]
  cb.prior <- rbeta(Nrep,alpha,beta)

#  PLOT PRIORS SO WE CAN LOOK AT THEM:
  plot(x,dbeta(x,alpha,beta),type="l",ylim=c(0,80),xlab="",ylab="")
  title(sub=paste("Beta(110w,7w) prior with w =",w[j]),cex=1.05)

for(k in 1:length(N)){
	
cb.x <- rbinom(Nrep,N[k],cb.prior)

#  HERE IS THE CALCULATION OF THE LOWER .025 POINT OF THE POSTERIOR:
cb.post <- rep(NA,Nrep)
cb.post <- qbeta(0.025,cb.x + alpha,N[k] - cb.x + beta)

plot(x,dbeta(x,cb.x[1] + alpha,N[k] - cb.x[1] + beta),type="l",ylim=c(0,80),
  xlab="",ylab="")
for(i in 2:Nrep){
  lines(x,dbeta(x,cb.x[i] + alpha,N[k] - cb.x[i] + beta),lty=i)}
title(sub=paste("Nrep=",Nrep,", N=",N[k],", w =",w[j]),cex=1.05)


cat(paste("\n N = ", N[k], "   w = ", w[j]))
for(i in 1:length(C)){
  pct[i,1+(j-1)*length(N)+k] <- sum((cb.post > C[i])==T)/Nrep
  cat(paste("\n Percent of lower .025 points >",C[i]," is", 
    round(pct[i,1+(j-1)*length(N)+k],3)))
}
#cat("\n Hit enter for more...\n")
#readline()

}  # end of N[k] loop

}  # end of w[j] (prior) loop

mtext(side=3,line=1,cex=1.4,outer=TRUE,
  "Priors and simulated posteriors, Chronicle B study, beta-binomial design")
mtext(side=3,line=-2,cex=1.4,outer=TRUE,
  "using weighted Beta(a,b) priors (target->downweighted)")

cat(paste("\n\nHere is sample size table for Beta(110w,7w) prior:\n"))
print(pct)

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