###############################################################
#Compute sample size that does not excced a given total error:#
###############################################################
#alpha= (Avg.) Total (type I + type II) error (TE)
#w, 1-w are weights for average type I and II errors, respectively.
#n.min=minimum (starting) value of sample size
#n.steps=number of steps required after n.min
#Copyright (2009) Sujit K Ghosh, NC State University [Bayesian Biostat course]
##############################################################################

##################################################
#This function computes the required sample size #
#for a given test, total error and weights       #
##################################################
sscal=function(alpha=0.25,w=0.1,n.min=10,n.steps=100){
n.max=n.min+n.steps
TE=1
n.try=max(n.min,2)
while(TE>alpha){
TE=bayes.test(n=n.try,w=w)
n.try=n.try+1
if(n.try>n.max){return(c("required sample size>",n.max))}
}
n.opt=n.try-1
cat(c("n.opt=",n.opt,"[alpha=",alpha,"w=",w,"TE=",round(TE,4),"]"),fill=T)
return(n.opt)
}

######################################################
#This function computes the total error given a value#
#of sample size, weights and prior parameters        #
######################################################
bayes.test=function(n=30,w=0.1,theta0=0.5,a=1,b=1){
#X ~ Bin(theta, n), theta ~ Beta(a,b)
#Null Hyp, H0: theta<=theta0
#Alt. Hyp, H1: theta>theta0

x.obs=0:n
log.prior.prob0=pbeta(theta0,a,b,log.p=T)
log.prior.prob1=pbeta(theta0,a,b,lower.tail=F,log.p=T)
log.post.prob0=pbeta(theta0,a+x.obs,b+n-x.obs,log.p=T)
log.post.prob1=pbeta(theta0,a+x.obs,b+n-x.obs,lower.tail=F,log.p=T)
log.marginal=lchoose(n,x.obs)+lbeta(a+x.obs,b+n-x.obs)-lbeta(a,b)
logBF=log.post.prob1+log.prior.prob0-log.post.prob0-log.prior.prob1

#Compute Bayes Avg Type I Error:
AE1.bayes=function(v){
ind=ifelse(logBF>v,1,0)
ae1=sum(ind*exp(log.marginal+log.post.prob0-log.prior.prob0))
return(ae1)
                     }
#Compute Bayes Avg Type II Error:
AE2.bayes=function(v){
ind=ifelse(logBF<=v,1,0)
ae2=sum(ind*exp(log.marginal+log.post.prob1-log.prior.prob1))
return(ae2)
                      } 
opt.fn=function(v){
(1-w)*AE1.bayes(v)+w*AE2.bayes(v)
                  }
#Set range of cut-off values:
a.min=min(logBF)-1
b.max=max(logBF)+1

#Compute the Bayes test cutoff value:
cutoff.bayes=optimize(opt.fn,interval=c(a.min,b.max))$minimum

#Compute the power functions:
power.bayes=function(theta,cutoff=0){
ind=ifelse(logBF>cutoff.bayes,1,0)
value=sum(ind*dbinom(x.obs,size=n,prob=theta))
return(value) } 

#Compute total Bayes test error:
type1.bayes=AE1.bayes(cutoff.bayes)
type2.bayes=AE2.bayes(cutoff.bayes)
TE.bayes=round(type1.bayes+type2.bayes,4)
cat(c("AE1=",round(type1.bayes,4),"AE2=",round(type2.bayes,4),"TE=",TE.bayes),fill=T)

#Plot the power functions:
theta=seq(0.01,0.99,l=100)
power.bayes.val=lapply(theta,power.bayes,cutoff=cutoff.bayes)
plot(theta,power.bayes.val,ylab="power function",ylim=c(0,1),type="l")
abline(v=theta0,col="gray")
par(cex=0.75)
title(paste("Bayes Power fn.: w=",round(w,2),",n=",n,"TE=",TE.bayes))

return(TE.bayes) 
}