################ #### The code included here is used to produce Figure 2 in the paper. ################ ################ #### Define all the functions that we need for the main body ################ #### Compute lower limits and upper limits PD_limit<-function(lambda,n,alpha,df1,df2){ tempcp=0 temparea=rep(0,(n+1)) cv=qchisq(1-alpha,1) L=0;U=0 for(x in 0:n){ phat=x/n;qhat=1-phat if(lambda==0){ if(x==0){ L[x+1]=0;U[x+1]=1-exp(-cv/(2*n)) } else if(x==n){ L[x+1]=exp(-cv/(2*n));U[x+1]=1 } else { powerd=function(p) 2*n*(phat*log(phat/p)+(1-phat)*log((1-phat)/(1-p)))-cv #### Find lower limits dis=1;maxiter=100000;FuncValue=1 steps=1;LeftP=0;RightP=phat while(dis>1e-15&steps1e-15){ bisec=0.5*(LeftP+RightP) value=powerd(bisec) if(value>0) LeftP=bisec else RightP=bisec steps=steps+1 dis=abs(LeftP-RightP) FuncValue=abs(value) } if(abs(bisec)<1e-10) bisec=0 L[x+1]=max(bisec,0) #### Find upper limits dis=1;maxiter=100000;FuncValue=1 steps=1;LeftP=phat;RightP=1 while(dis>1e-15&steps1e-15){ bisec=0.5*(LeftP+RightP) value=powerd(bisec) if(value<0) LeftP=bisec else RightP=bisec steps=steps+1 dis=abs(LeftP-RightP) FuncValue=abs(value) } if(abs(1-bisec)<1e-10) bisec=1 U[x+1]=min(bisec,1) } } else if(lambda==1){ L[x+1]=max((x+0.5*cv)/(n+cv)-sqrt(0.25*cv^2+n*phat*qhat*cv)/(n+cv),0) U[x+1]=min((x+0.5*cv)/(n+cv)+sqrt(0.25*cv^2+n*phat*qhat*cv)/(n+cv),1) } else { ALambda=lambda*(lambda+1)/(2*n)*cv+1 powerd=function(p) 2*n/(lambda^2+lambda)*((phat^(lambda+1))/(p^lambda)+((1-phat)^(lambda+1))/((1-p)^lambda)-1)-cv if(x==0){ L[x+1]=0 U[x+1]=1-ALambda^(-1/lambda) } else if(x==n){ L[x+1]=ALambda^(-1/lambda) U[x+1]=1 } else { #### Find lower limits dis=1;maxiter=100000;FuncValue=1 steps=1;LeftP=0;RightP=phat while(dis>1e-15&steps1e-15){ bisec=0.5*(LeftP+RightP) value=powerd(bisec) if(value>0) LeftP=bisec else RightP=bisec steps=steps+1 dis=abs(LeftP-RightP) FuncValue=abs(value) } if(abs(bisec)<1e-10) bisec=0 L[x+1]=max(bisec,0) #### Find upper limits dis=1;maxiter=100000;FuncValue=1 steps=1;LeftP=phat;RightP=1 while(dis>1e-15&steps1e-15){ bisec=0.5*(LeftP+RightP) value=powerd(bisec) if(value<0) LeftP=bisec else RightP=bisec steps=steps+1 dis=abs(LeftP-RightP) FuncValue=abs(value) } if(abs(1-bisec)<1e-10) bisec=1 U[x+1]=min(bisec,1) } } } return(data.frame(L,U)) } Bayes_limit<-function(n,alpha,a,b,df1,df2){ tempcp=0 L=0;U=0 for(x in 0:n){ phat=x/n;qhat=1-phat if(x==0){ L[x+1]=0 U[x+1]=qbeta(1-(alpha/2), a, n+b, ncp = 0, lower.tail = TRUE, log.p = FALSE) } else if(x==n){ L[x+1]=qbeta(alpha/2, n+a, b, ncp = 0, lower.tail = TRUE, log.p = FALSE) U[x+1]=1 } else { L[x+1]=qbeta(alpha/2, x+a, n-x+b, ncp = 0, lower.tail = TRUE, log.p = FALSE) U[x+1]=qbeta(1-(alpha/2), x+a, n-x+b, ncp = 0, lower.tail = TRUE, log.p = FALSE) } } return(data.frame(L,U)) } #### Compute the coverage probability and length PD_CP<-function(lambda,n,alpha){ cp<-0;length<-0 for(j in 1:length(truep)){ temp1<-rep(NA,n+1);temp2<-rep(NA,n+1) for(x in 0:n){ temp1[x+1]=dbinom(x, n, truep[j])*(truep[j]>limit_PD$L[x+1])*(truep[j]limit_LR$L[x+1])*(truep[j]limit_Bayes$L[x+1])*(truep[j]