#Gehans' Two stage design:

sscal.gehan=function(theta0=0.2,alpha0=0.05,a=1,b=1,n.min=2,n.max=5000){
condl.prob=1
ntry=n.min
while(condl.prob>alpha0){
condl.prob=pbeta(theta0,a,b+ntry,lower.tail=F)
ntry=ntry+1
if(ntry>n.max){return(c("required sample size>",n.max))}
}
n.opt=ntry-1
cat(c("n.opt=",n.opt,"(theta0=",theta0,"alpha0=",alpha0,")"),fill=T)
#return(n.opt)
}

#Simons two-stage design:

prob.fn=function(theta,r1,n1,r,n){
if(r1>n1){return("r1<n1 violated")}
x=(r1+1):n1
out=sum(pbinom(r-x,size=n-n1,prob=theta,lower.tail=F)*dbinom(x,size=n1,prob=theta))
return(out)
}

prob.success=function(theta,r10,n10,r0,n0){
out=sapply(theta,prob.fn,r1=r10,n1=n10,r=r0,n=n0)
return(out)
}

prob.failure=function(theta,r10,n10,r0,n0){
out=1-prob.success(theta,r10=r10,n10=n10,r0=r0,n0=n0)
return(out)
}

FNR=function(r1=0,n1=9,r=2,n=24,theta0=0.05){
num=integrate(prob.success,lower=0,upper=theta0,r10=r1,n10=n1,r0=r,n0=n)$value
den=integrate(prob.success,lower=0,upper=1,r10=r1,n10=n1,r0=r,n0=n)$value
out=num/den
return(out)
}
#Compare FNP to alpha=0.05

FPR=function(r1=0,n1=9,r=2,n=24,theta1=0.25){
num=integrate(prob.failure,lower=theta1,upper=1,r10=r1,n10=n1,r0=r,n0=n)$value
den=integrate(prob.failure,lower=0,upper=1,r10=r1,n10=n1,r0=r,n0=n)$value
out=num/den
return(out)
}
#Compare FPR to beta=0.2