#Bayes tests: X ~ Bin(n, theta); H0: theta<=theta0 vs. H1:theta>theta0

compare.tests=function(n=20,w=0.1,alpha=0.05,theta0=0.5,a=1,b=1){

#Compute the Bayes test cutoff value:

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 frequentist test cutoff value:
cutoff.freq=qbinom(1-alpha,prob=theta0,size=n)

#Compute Freq Avg Type I Error:
AE1.freq=function(v){
ind=ifelse(x.obs>v,1,0)
ae1=sum(ind*exp(log.marginal+log.post.prob0-log.prior.prob0))
return(ae1)
                     }
#Compute Freq Avg Type II Error:
AE2.freq=function(v){
ind=ifelse(x.obs<=v,1,0)
ae2=sum(ind*exp(log.marginal+log.post.prob1-log.prior.prob1))
return(ae2)
                     } 
#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)                       }

power.freq=function(theta,cutoff=0){
ind=ifelse(x.obs>cutoff.freq,1,0)
value=sum(ind*dbinom(x.obs,size=n,prob=theta))
return(value)                       }

#Compute total frequentist test error:
type1.freq=AE1.freq(cutoff.freq)
type2.freq=AE2.freq(cutoff.freq)
TE.freq=round(type1.freq+type2.freq,4)

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

#print results on screen:
cat(c("Total Avg. Errors: [n=",n,"alpha=",alpha,"weight=",w,"]"),fill=T)
cat(c("TE.freq=",TE.freq," and TE.bayes=",TE.bayes),fill=T)

#plot the power functions:
theta=seq(0.01,0.99,l=100)
power.freq.val=lapply(theta,power.freq,cutoff=cutoff.freq)
power.bayes.val=lapply(theta,power.bayes,cutoff=cutoff.bayes)
plot(theta,power.bayes.val,ylab="power function",ylim=c(0,1),type="l")
lines(theta,power.freq.val,lty=2,col=2)
abline(v=theta0,col="gray")
abline(h=alpha,col="gray",lty=2)
par(cex=0.65)
legend(x=0.01,y=0.99,legend=c("bayes","freq"),lty=1:2,col=c(1,2))
title(paste("Comparing power fn.s: alpha=",round(alpha,2),",w=",round(w,2),",n=",n))
}

par(mfrow=c(2,2))
for(n in c(10,20,40,80)){
compare.tests(n=n)}

#Results:
#Total Avg. Errors: [n= 10 alpha= 0.05 weight= 0.1 ]
#TE.freq= 0.6387  and TE.bayes= 0.4688

#Total Avg. Errors: [n= 20 alpha= 0.05 weight= 0.1 ]
#TE.freq= 0.432  and TE.bayes= 0.3442

#Total Avg. Errors: [n= 40 alpha= 0.05 weight= 0.1 ]
#TE.freq= 0.2734  and TE.bayes= 0.2304

#Total Avg. Errors: [n= 80 alpha= 0.05 weight= 0.1 ]
#TE.freq= 0.1895  and TE.bayes= 0.1895