LISTINGS OF SPLUS FUNCTIONS, AND SPLUS LOG, FOR CALCULATION 
  AND PLOTTING OF BINOMIAL-PROBABILITY CONFIDENCE INTERVAL 
                COVERAGE PROBABILITIES

### START WITH 
> options(gui="motif")
> motif()


####                          NOW COME FUNCTION-LISTINGS

> ScorInt
function (y,n,alpha, upr=T) 
{
### Score-statistic upper/lower confidence limit
   phat <- y/n
   fr <- n/(n+qnorm(1-alpha/2)^2)
   phat*fr + 0.5*(1-fr)+(2*as.numeric(upr)-1)*sqrt(
       (1-fr)*(phat*(1-phat)*fr+0.25*(1-fr)))
}
> WaldInt
function (y,n,alpha, upr=T) 
{
### Upper or Lower Wald-statistic confidence limit
   phat <- y/n
   phat + (2*as.numeric(upr)-1)*qnorm(1-alpha/2)*
      sqrt(phat*(1-phat)/n)
}
> ClPint
function (y,n,alpha, upr=T) 
{
### Clopper-Pearson upper/lower confidence limits
   if(y==0) y <- 1.e-8 else if(y==n) y<- n-1.e-8
   1/(1+ if(upr) 
(n-y)/((y+1)*qf(1-alpha/2,2*(y+1),2*(n-y)))
  else (n-y+1)/(y*qf(alpha/2,2*y,2*(n-y+1)))) 
}

> BinCovrg
function (n, p, alpha) 
{
## Calculation of coverage probabilities under the
##   Binom(n,p) model for the three confidence intervals
##   on  p  calculated within the function  BinomInt.
## First step is to find the values y (not nec integer)
##   solving  ConfLim(y)=p  
  Scory <- Waldy <- numeric(2) -> ClPy
  Waldy[1] <- if(WaldInt(0,n,alpha,F)>p) -1 else
    if(WaldInt(n,n,alpha,F)<p) n+1 else uniroot(
       function(y,n,alpha,pp) WaldInt(y,n,alpha,F)-pp,
          c(0,n), n=n, alpha=alpha, pp=p)$root
  Waldy[2] <- if(WaldInt(0,n,alpha,T)>p) -1 else
    if(WaldInt(n,n,alpha,T)<p) n+1 else uniroot(
       function(y,n,alpha,pp) WaldInt(y,n,alpha,T)-pp,  
          c(0,n), n=n, alpha=alpha, pp=p)$root
  Scory[1] <- if(ScorInt(0,n,alpha,F)>p) -1 else
    if(ScorInt(n,n,alpha,F)<p) n+1 else uniroot(
       function(y,n,alpha,pp) ScorInt(y,n,alpha,F)-pp,
          c(0,n), n=n, alpha=alpha, pp=p)$root
  Scory[2] <- if(ScorInt(0,n,alpha,T)>p) -1 else
    if(ScorInt(n,n,alpha,T)<p) n+1 else uniroot(
       function(y,n,alpha,pp) ScorInt(y,n,alpha,T)-pp,
          c(0,n), n=n, alpha=alpha, pp=p)$root
  ClPy[1] <- if(ClPint(0,n,alpha,F)>p) -1 else
    if(ClPint(n,n,alpha,F)<p) n+1 else uniroot(
       function(y,n,alpha,pp) ClPint(y,n,alpha,F)-pp, 
          c(0,n), n=n, alpha=alpha, pp=p)$root
  ClPy[2] <- if(ClPint(0,n,alpha,T)>p) -1 else
    if(ClPint(n,n,alpha,T)<p) n+1 else uniroot(
       function(y,n,alpha,pp) ClPint(y,n,alpha,T)-pp,
          c(0,n), n=n, alpha=alpha, pp=p)$root
  list(ClP.prob = pbinom(ceiling(ClPy[1])-1,n,p)-
        pbinom(ceiling(ClPy[2])-1,n,p),
     Score.prob = pbinom(ceiling(Scory[1])-1,n,p)-
        pbinom(ceiling(Scory[2])-1,n,p), Wald.prob = 
        pbinom(ceiling(Waldy[1])-1,n,p)-pbinom(ceiling(
        Waldy[2])-1,n,p), Waldy=Waldy, Scory=Scory, ClPy=ClPy)
}

> CritVecs
function (n, alpha) 
{
### This function calculates points of discontinuity for the
### 3 intervals under study, and puts their union in vector
  ys <- 1:(n-1)
  props <- ys/n
  znom <- qnorm(1-alpha/2)
  fr <- n/(n+znom^2)
  Newvec <- function(Lovec, Hivec) {
      tmpvec <- c(Lovec, Hivec)
      sort(tmpvec[!duplicated(tmpvec)])
  }
  WaldHi <- pmin(1,props+znom*sqrt(props*(1-props)/n))
  WaldLo <- pmax(0, props-znom*sqrt(props*(1-props)/n))
  Waldvec <- Newvec(WaldLo,WaldHi)
  ScorLo <- props*fr + 0.5*(1-fr)-sqrt(
       (1-fr)*(props*(1-props)*fr+0.25*(1-fr)))
  ScorHi <- props*fr + 0.5*(1-fr)+sqrt(
       (1-fr)*(props*(1-props)*fr+0.25*(1-fr)))
  Scorvec <- Newvec(ScorLo, ScorHi)
  ClPLo <- 1/(1+(n-ys+1)/(ys*qf(alpha/2,2*ys,2*(n-ys+1))))
  ClPHi <- 1/(1+ (n-ys)/((ys+1)*qf(1-alpha/2,2*(ys+1),2*(n-ys))))
  ClPvec <- Newvec(ClPLo, ClPHi) 
  list(Waldvec=Waldvec, Scorvec=Scorvec, ClPvec=ClPvec)
}

#### Now can use these functions    with n=40

> auxv <- CritVecs(40,.05)
   xv <- union(seq(.005,1,.005), union(auxv$Waldvec, 
     union(auxv$Scorvec, auxv$ClPvec)))  ## now length 429
> xv <- sort(xv)[2:428]
   ypv <- array(0, dim=c(427,3))
> for (i in 1:427) ypv[i,] <- unlist(BinCovrg(40,xv[i],.05)[1:3])
> plot(xv, ypv[,1], ylim=c(.7,1), type="l", xlab="True p", 
       ylab="Coverage Probability",  
       main="Binomial Confidence Interval Coverage, n=40")
   lines(xv, ypv[,2], lty = 3)
   lines(xv, ypv[,3],lty=5)