# This function compute the confidence interval for Clopper_Pearson Interval CP <- function(al,x,n) { x1=n-x+1 x2=n-x if(x!=0) { f1=qf(al/2, 2*x,2*x1, lower.tail = TRUE, log.p = FALSE) L=(1+ x1/(x*f1))^-1 } if(x!=n) { f2=qf(1-al/2, 2*(x+1),2*x2, lower.tail = TRUE, log.p = FALSE) U=(1+ x2/((x+1)*f2))^-1 } if ( x==0) L=0 if (x==n) U=1 return(as.vector(c(L, U))) } #This part of the program create matrix of Clopper_Pearson limits mat1 <- function(n,al) { if(abs(al-0.05)<=0.05) alpha=seq(0.05,0.3,0.005) if(abs(al-0.01)<=0.01) alpha=seq(0.01,0.5,0.01) mm=length(alpha) MLc=matrix(0,n+1,mm) MUc=matrix(0,n+1,mm) aa=matrix(0,n+1,mm) for (k in 1:mm) { for (i in 1:(n+1)) { x=as.integer(i-1) out1=CP(alpha[k],x,n) MLc[i,k]=out1[1] MUc[i,k]=out1[2] } } tt=cbind(as.vector(MLc),as.vector(MUc)) return(tt) } #this function computes the lower and upper bound of model based confidence interval # by Desale Habtzghi, Ashish Das and Chand Midha New <- function(AL,x1, n) { tt=mat1(n,AL) m=10^10 ## Alpha values to fit intervals #alpha=seq(0.001,0.1,0.001) # alpha=seq(0.001,0.5,0.005) alpha=seq(0.05,0.3,0.005) #alpha=c(0.01,0.05) mm=length(alpha) MLc=matrix(tt[,1],n+1,mm) MUc= matrix(tt[,2],n+1,mm) #### Here we are computing intervals for different Ps and CIs p1=seq(0.01,0.99, 0.01) # probability of success nn=length(p1) summ2=0 fc=numeric(n+1) #p1=numeric(1) Lee=numeric(mm) Uee=numeric(mm) fail=matrix(0,mm,nn) succ=matrix(0,mm,nn) weight=matrix(0,mm,nn) aa=matrix(0,mm,nn) pp=matrix(0,mm,nn) nfail=numeric(mm) nsucc=numeric(mm) nfail=numeric(mm) ## the most outer loop fnew=numeric(n+1) sumL0=0 L=numeric(1) U=numeric(1) #aa=matrix(0,n+1) #pp1=matrix(0,n+1) ## Outer loop for our method ################################# # This loop computes the confidence interval for (k in 1:mm) { for (j in 1:(nn)) { summ2=0 ### fe= frequency of random numbers ### Confidence intervals ### x=number of successes #check coverage for a given p for ( i in 1:(n+1)) { x=as.integer(i-1) fc[i]=dbinom(x,n,p1[j])*m #fe[i]=as.integer(dbinom(x,n,p1[j])*m) if ((MLc[i,k] <= p1[j])&& (p1[j]<=MUc[i,k])) summ2=summ2+fc[i] } # end of the i loop fail[k,j]=m-summ2 # frequency of coverage for a given p succ[k,j]=summ2 # gets the lower bound for given x1 to find the coef of our estimate Lee[k]=MLc[x1+1,k] Uee[k]=MUc[x1+1,k] # frequency of coverage for a given p succ[k,j]=summ2 # gets the lower bound for given x1 to find the coef of our estimate aa[k,j]=alpha[k] # weights weight[k,j]=dbinom(x1,n,p1[j]) } #end of the j loop # computes the weighted coverage for all the possible p # values nsucc[k]=(succ[k,]%*%weight[1,])/sum(weight[1,]) } # end of the k loop delta=cbind(floor(nsucc),floor(m-nsucc)) #delta=nsucc/m fit1 <-glm(delta~Lee,family=binomial(logit)) fit2 <-glm(delta~Uee,family=binomial(logit)) fitp1 <-glm(delta~Lee,family=binomial(probit)) fitp2 <-glm(delta~Uee,family=binomial(probit)) B00=coef(fit1)[1] B11=coef(fit1)[2] UB00=coef(fit2)[1] UB11=coef(fit2)[2] B0=B00 B1=B11 UB0=UB00 UB1=UB11 # end the the most outer loop L[1]=(log((1-AL)/AL)-B0)/B1 U[1]=(log((1-AL)/AL)-UB0)/UB1 if( is.na(L[1])) L[1]=0 if( is.na(U[1])) U[1]=1 outt=as.vector(c(L[1],U[1])) return(outt) } # end of the main function function # New(alpha,x,n) # Arguments # 1. Alpha significan level # 2. x number of success # 3. n sample size ( number of trials) #Examples x=4; n=10; alpha=0.05 Limits=New(alpha,x,n) #print( "Model based confidence interval") CI=c(Limits[1],Limits[2]) CI