# This program is written by Weizhen Wang to compute the intervals in "On Exact Interval 
# Estimation for the Odds Ratio in the Subject-specific Table" by Wang (2017). All rights reserved.

# 1. Here X \sim Bin(n,p). 

# 2.  Here, n is known, we follow Wang (2014, Statistica Sinica) and estimate p based on X.
# i) one-sided CI for p, [L(X),1]
# ii) one-sided CI for p, [0,U(X)]
# iii) two-sided CI for p, [L(X),U(X)]
# iv) U(X)=1-L(n-X)

# 3. We provide Bino(n,X,conf.level,CIType,precision)
# For example, Bino(20,3,0.95,"lower",0.0001), #the smallest 95% lower one-sided interval for x=3
# Bino(20,3,0.95,"upper",0.0001),  #the smallest 95% upper one-sided interval for x=3 
# Bino(20,3,0.95,"two.sided",0.0001),  # It contains x and 
# an admissible 95% two-sided interval (2nd and 3rd) for x=3.

# 4. We convert Bino(n,X,conf.level,CIType,precision) to oddsr(nstar,n01,conf.level,CIType,precision),
# which is the confidence interval (16) in the paper for r, the odds ratio in (3). 
# The function output includes the observation n_01 and its interval for r, 
# and all n01-values and their intervals for r (n^*+1 intervals in total).

# 5. The first line of  oddsr(53,37,0.95,"two.sided",0.0001) is 37 1.255127e+00 4.346806e+00
#    Run oddsr(53,37,0.8744,"two.sided",0.0001) for Example 1.

# 6. For example,
#    oddsr(600,37,0.95,"two.sided",0.0001) takes 300 minutes
#    oddsr(400,37,0.95,"two.sided",0.0001) takes 100 minutes
#    oddsr(800,37,0.95,"two.sided",0.0001) takes 600 minutes

# 7. One-sided intervals (12) and (13) for r can also be obtained by changing 
# the value of CIType to "upper" or "lower".   

Bino<-function(n,X,conf.level,CIType,precision){    #compute CI for pm
n=n
alpha=1-conf.level
prec1=precision
prec2= 0.000000001 # precision for the minimum coverage 

lcp<- function(x,n){   # the lower limit of 1-alpha CP interval using F
aa=1:length(x)
for (i in aa){
if (x[i]<0.5){aa[i]=0}
else {aa[i]=(1+(n[i]-x[i]+1)/(x[i]*qf(alpha/2,2*x[i],2*(n[i]-x[i]+1))))^(-1)}
}
aa
}

ucp<- function(x,n){ # the upper limit of 1-alpha CP interval
1-lcp(n-x,n)
}

lcp1<- function(x,n){   # the lower limit of 1-2alpha CP interval using F
aa=1:length(x)
for (i in aa){
if (x[i]<0.5){aa[i]=0}
else {aa[i]=(1+(n[i]-x[i]+1)/(x[i]*qf(alpha,2*x[i],2*(n[i]-x[i]+1))))^(-1)}
}
aa
}

ucp1<- function(x,n){ # the upper limit of 1-2alpha CP interval
1-lcp1(n-x,n)
}

#indicator fucntion of interval [a,b]
ind<- function(x,a,b){
(x>=a)*(x<=b)
}

# the coverage probability function for the interval [lcin,ucin]. 
cpcig<- function(p,lcin,ucin){
kk=1:length(p)
for (i in kk){
xx<- 0:n
indp=xx
uu=1 
while (uu<length(xx)+0.5) {
indp[uu]=ind(p[i],lcin[uu],ucin[uu])*dbinom(xx[uu],n,p[i])
uu=uu+1}
kk[i]=sum(indp)
}
kk
}

xx=seq(0,n,by=1)
nn=n+0*xx
lcpn=lcp(xx,nn)   #the lower limit of CP interval
ucpn=ucp(xx,nn)   #the upper limit of CP interval
lcp1n=lcp1(xx,nn) #the lower limit of one-sided CP interval >=lcpn
ucp1n=ucp1(xx,nn) #the upper limit of one-sided CP interval <=ucpn

lcpw=lcpn         #the lower limit of CP refinement, should be between lcpn and lcp1n
ucpw=ucpn         #the upper limit of CP refinement, should be between ucpn and ucp1n

if(CIType=="upper"){
BinoL=0
BinoU=ucp1n[X+1]
result=cbind(xx,BinoL+0*xx,ucp1n)
result=rbind(cbind(X,BinoL,BinoU),result)
#result=c(X,BinoL,BinoU)
}else if(CIType=="lower"){
BinoL=lcp1n[X+1]
BinoU=1
result=cbind(xx,lcp1n, BinoU+0*xx)
result=rbind(cbind(X,BinoL,BinoU),result)
#result=c(X,BinoL,BinoU)
}else if(CIType=="two.sided"){  #*

### compute lcpw and ucpw at X=n/2+1

if(n/2==floor(n/2)){#**
ss=n/2+1
a0=lcpn[ss]
a1=lcp1n[ss]  # a1>a0. the final lower limit is between a1 and a0.
while (a1-a0>prec1){
aa=(a0+a1)/2
lcpw[ss]=aa
ucpw[ss]=1-aa
pp=c((lcpw-prec2)*((lcpw-prec2)>0),ucpw[1:n]+prec2)
cp=cpcig(pp,lcpw,ucpw)
minic=min(cp)
if (minic>=1-alpha){
a0=aa
}
else{a1=aa}
}
lcpw[ss]=a0
ucpw[ss]=1-a0


###

#### compute lcpw and ucpw at X=n/2 down to 0
ss=n/2
while (ss>.5){      # construct the lower and upper limit from X=n/2 to 0

a0=lcpw[ss]
a1=lcp1n[ss]  # a1>a0. the final lower limit is between a1 and a0.
b0=ucp1n[ss]
b1=ucpw[ss]   # b1>b0. the final upper limit is between b1 and b0.

  while (b1-b0>prec1){ # solve the upper limit
   bb=(b0+b1)/2
   ucpw[ss]=bb
   lcpw[n+2-ss]=1-bb
pp=c((lcpw-prec2)*((lcpw-prec2)>0),ucpw[1:n]+prec2)   ####
cp=cpcig(pp,lcpw,ucpw)
  minic=min(cp)
  if (minic>=1-alpha){
  b1=bb
  }
 else {b0=bb}
}
ucpw[ss]=b1
lcpw[n+2-ss]=1-b1

  while (a1-a0>prec1){ # solve the lower limit
   aa=(a0+a1)/2
   lcpw[ss]=aa
   ucpw[n+2-ss]=1-aa
pp=c((lcpw-prec2)*((lcpw-prec2)>0),ucpw[1:n]+prec2)
cp=cpcig(pp,lcpw,ucpw)
  minic=min(cp)
  if (minic>1-alpha){
  a0=aa
  }
  else{a1=aa}
}
lcpw[ss]=a0
ucpw[n+2-ss]=1-a0
ss=ss-1
}
}else{ #**
ss=floor(n/2)+1
while (ss>.5){      # construct the lower and upper limit from X=n/2 to 0

a0=lcpw[ss]
a1=lcp1n[ss]  # a1>a0. the final lower limit is between a1 and a0.
b0=ucp1n[ss]
b1=ucpw[ss]   # b1>b0. the final upper limit is between b1 and b0.

  while (b1-b0>prec1){ # solve the upper limit
   bb=(b0+b1)/2
   ucpw[ss]=bb
   lcpw[n+2-ss]=1-bb
pp=(lcpw-prec2)*((lcpw-prec2)>0)
#pp=c((lcpw-prec2)*((lcpw-prec2)>0),ucpw[1:n]+prec2)
cp=cpcig(pp,lcpw,ucpw)
  minic=min(cp)
  if (minic>=1-alpha){
  b1=bb
  }
 else {b0=bb}
}
ucpw[ss]=b1
lcpw[n+2-ss]=1-b1

  while (a1-a0>prec1){ # solve the lower limit
   aa=(a0+a1)/2
   lcpw[ss]=aa
   ucpw[n+2-ss]=1-aa
pp=(lcpw-prec2)*((lcpw-prec2)>0)
#pp=c((lcpw-prec2)*((lcpw-prec2)>0),ucpw[1:n]+prec2)
cp=cpcig(pp,lcpw,ucpw)
  minic=min(cp)
  if (minic>1-alpha){
  a0=aa
  }
  else{a1=aa}
}
lcpw[ss]=a0
ucpw[n+2-ss]=1-a0

ss=ss-1
}
}#**
BinoLcp=lcpn[X+1]
BinoUcp=ucpn[X+1]
BinoLw=lcpw[X+1]
BinoUw=ucpw[X+1]
result=cbind(xx,lcpw, ucpw)
result=rbind(cbind(X,BinoLw,BinoUw),result)
#result=c(X,BinoLw,BinoUw)
#result=cbind(X,BinoLcp,BinoUcp,BinoLw,BinoUw)
}else{  #*
result="Type error!" 
}
result
}

oddsr<-function(nstar,n01,conf.level,CIType,precision){  # compute CI for r
  aa=Bino(nstar,n01,conf.level,CIType,precision)
  bb=aa/(1-aa)
  cbind(aa[,1],bb[,2],bb[,3])
}

tim=proc.time()
oddsr(53,37,0.95,"two.sided",0.0001)
oddsr(400,37,0.95,"two.sided",0.0001)
(proc.time()-tim)/60