*! Date    : 08 May 2018
*! Version : 1.0

/*
 08/05/18 v1.0 Initial version complete
*/

program define ci_fixed, rclass
version 11.0
syntax, [n(integer 25) pi(numlist) alpha(real 0.05) method(string) ///
         SUMmary(integer 0) PLot(string) *]

local xopt `"`options'"'
preserve

if ("`method'" == "") {
  local method "all"
}

///// Perform checks on input variables ////////////////////////////////////////

if (`n' <= 0) {
  di "{error} n must be an integer in (0,Inf)."
  exit(198)
}
if ("`pi'" ~= "") {
  local lenpi:list sizeof pi
  forvalues i = 1/`lenpi' {
    local pii:word `i' of `pi'
    if (`pii' > 1 | `pii' < 0) {
      di "{error} Elements in pi must belong to [0,1]."
      exit(198)
    }
  }
}
if ((`alpha' <= 0) | (`alpha' >= 1)) {
  di "{error} alpha must be a real strictly between 0 and 1."
  exit(198)
}
if ("`method'" ~= "all" & "`method'" ~= "agresti_coull" & "`method'" ~= "clopper_pearson" & "`method'" ~= "jeffreys" & "`method'" ~= "mid_p" & "`method'" ~= "wald" & "`method'" ~= "wilson") {
  di "{error} method must be one of: all, agresti_coull, clopper_pearson, jeffreys, mid_p, wald, or wilson."
  exit(198)
}
// Set up matrices to pass to mata
if ("`pi'" ~= "") {
  local matapi ""
  foreach i of local pi{
    if "`matapi'" == "" local matapi "`i'"
    else local matapi "`matapi',`i'"
  }
  mat pi = (`matapi')
}
else {
  mat pi = .
}
if ("`plot'" ~= "" ) {
  if ("`plot'" ~= "coverage" & "`plot'" ~= "length") {
    di "{error} plot must be one of: coverage and length."
    exit(198)
  }
}

///// Compute and Output ///////////////////////////////////////////////////////

mata: CIFixed(`n', `alpha', "`method'", `summary', "`plot'", `"`xopt'"')
matrix colnames ci = s m Method "clow(s,m)" "cupp(s,m)" "l(s,m)"
return mat ci = ci
matrix colnames perf = pi Method "bar(L)" "max(L)" "E(L|pi)" "Var(L|pi)" "Cover(C|pi)"
return mat perf = perf
return mat n = n
return mat pi = pi
return mat alpha = alpha
return mat summary = summary

restore

end

///// Mata  ////////////////////////////////////////////////////////////////////

mata:

void CIFixed(n, alpha, method, summary, plot, xopt)
{

  pi = st_matrix("pi")
  if (pi == .) {
    pi = 0.01*(0::100)
  }

  ///// Print Summary //////////////////////////////////////////////////////////



  ///// Main Computations //////////////////////////////////////////////////////

  if (method == "all") {
    int_method = 1::6
  }
  else if (method == "agresti_coull") {
    int_method = 1
  }
  else if (method == "clopper_pearson") {
    int_method = 2
  }
  else if (method == "jeffreys") {
    int_method = 3
  }
  else if (method == "mid_p") {
    int_method = 4
  }
  else if (method == "wald") {
    int_method = 5
  }
  else {
    int_method = 6
  }
  s  = 0::n
  n  = J(n + 1, 1, n)
  ci = (vec(J(length(int_method), 1, s')), vec(J(length(int_method), 1, n')), J(n[1] + 1, 1, int_method), J(length(int_method)*(n[1] + 1), 2, 0))
  for (i = 1; i <= rows(ci); i++) {
	if (ci[i, 3] == 1) {
	  ci[i, 4::5] = ci_fixed_agresti_coull(ci[i, 1], ci[i, 2], alpha)
	}
	else if (ci[i, 3] == 2) {
	  ci[i, 4::5] = ci_fixed_clopper_pearson(ci[i, 1], ci[i, 2], alpha)
	}
	else if (ci[i, 3] == 3) {
	  ci[i, 4::5] = ci_fixed_jeffreys(ci[i, 1], ci[i, 2], alpha)
	}
	else if (ci[i, 3] == 4) {
	  ci[i, 4::5] = ci_fixed_mid_p(ci[i, 1], ci[i, 2], alpha)
	}
	else if (ci[i, 3] == 5) {
	  ci[i, 4::5] = ci_fixed_wald(ci[i, 1], ci[i, 2], alpha)
	}
	else if (ci[i, 3] == 6) {
	  ci[i, 4::5] = ci_fixed_wilson(ci[i, 1], ci[i, 2], alpha)
	}
	if (summary == 1 & mod(i, 100) == 0) {
	  printf("{txt}...{res}%g{txt} confidence intervals determined...\n", i)
	}
  }
  ci = (ci, ci[, 5] - ci[, 4])
  dbinomial_pi = J(n[1] + 1, length(pi), 0)
  for (i = 1; i <= length(pi); i++) {
    dbinomial_pi[, i] = binomialp(n[1], 0::n[1], pi[i])
  }
  perf = (J(length(int_method), 1, pi), vec(J(length(pi), 1, int_method')), J(length(pi)*length(int_method), 5, 0))
  for (i = 1; i <= rows(perf); i++) {
    ci_i     = select(ci, ci[, 3] :== perf[i, 2])
	pi_index   = mm_which(pi :== perf[i, 1])
	perf[i, 3] = sum(dbinomial_pi[, pi_index]:*ci_i[, 6])
	perf[i, 4] = sum(dbinomial_pi[, pi_index]:*(ci_i[, 6]:^2)) - perf[i, 3]^2
	perf[i, 5] = mean(ci_i[, 6])
	perf[i, 6] = max(ci_i[, 6])
	coverage   = select(ci_i, ci_i[, 4] :<= perf[i, 1] :& ci_i[, 5] :>= perf[i, 1])
	perf[i, 7] = sum(dbinomial_pi[coverage[, 1] :+ 1, pi_index])
	if (summary == 1 & mod(i, 100) == 0) {
	  printf("{txt}...performance for {res}%g{txt} elements of pi evaluated...\n", i)
	}
  }
  if (plot == "length") {
	if (method ~= "all") {
	  st_matrix("perf", perf)
	  stata("qui svmat perf")
	  stata(`"twoway line perf3 perf1, xtitle({&pi}) ytitle(E({it:L}|{&pi}))"'+ xopt)
	}
	else {
	  perf1 = select(perf, perf[, 2] :== 1)
	  perf2 = select(perf, perf[, 2] :== 2)
	  perf3 = select(perf, perf[, 2] :== 3)
	  perf4 = select(perf, perf[, 2] :== 4)
	  perf5 = select(perf, perf[, 2] :== 5)
	  perf6 = select(perf, perf[, 2] :== 6)
      st_matrix("perf1", perf1)
      st_matrix("perf2", perf2)
      st_matrix("perf3", perf3)
      st_matrix("perf4", perf4)
      st_matrix("perf5", perf5)
	  st_matrix("perf6", perf6)
	  stata("qui svmat perf1")
	  stata("qui svmat perf2")
	  stata("qui svmat perf3")
	  stata("qui svmat perf4")
	  stata("qui svmat perf5")
	  stata("qui svmat perf6")
	  stata(`"twoway (line perf13 perf11, color(black)) (line perf23 perf21, color(blue)) (line perf33 perf31, color(red)) (line perf43 perf41, color(green)) (line perf53 perf51, color(yellow)) (line perf63 perf61, color(orange)), xtitle({&pi}) ytitle(E({it:L}|{&pi})) legend(lab(1 "Agresti-Coull") lab(2 "Clopper-Pearson") lab(3 "Jeffreys") lab(4 "Mid-p") lab(5 "Wald") lab(6 "Wilson"))"'+ xopt)
	}
  }
  else if (plot == "coverage") {
    if (method ~= "all") {
	  st_matrix("perf", perf)
	  stata("qui svmat perf")
	  stata(`"twoway line perf7 perf1, xtitle({&pi}) ytitle(E({it:L}|{&pi}))"'+ xopt)
	}
	else {
	  perf1 = select(perf, perf[, 2] :== 1)
	  perf2 = select(perf, perf[, 2] :== 2)
	  perf3 = select(perf, perf[, 2] :== 3)
	  perf4 = select(perf, perf[, 2] :== 4)
	  perf5 = select(perf, perf[, 2] :== 5)
	  perf6 = select(perf, perf[, 2] :== 6)
      st_matrix("perf1", perf1)
      st_matrix("perf2", perf2)
      st_matrix("perf3", perf3)
      st_matrix("perf4", perf4)
      st_matrix("perf5", perf5)
	  st_matrix("perf6", perf6)
	  stata("qui svmat perf1")
	  stata("qui svmat perf2")
	  stata("qui svmat perf3")
	  stata("qui svmat perf4")
	  stata("qui svmat perf5")
	  stata("qui svmat perf6")
	  stata(`"twoway (line perf17 perf11, color(black)) (line perf27 perf21, color(blue)) (line perf37 perf31, color(red)) (line perf47 perf41, color(green)) (line perf57 perf51, color(yellow)) (line perf67 perf61, color(orange)), xtitle({&pi}) ytitle({it:Cover}({it:C}|{&pi})) legend(lab(1 "Agresti-Coull") lab(2 "Clopper-Pearson") lab(3 "Jeffreys") lab(4 "Mid-p") lab(5 "Wald") lab(6 "Wilson"))"'+ xopt)
	}
  }

  ///// Return /////////////////////////////////////////////////////////////////

  if (summary == 1) {
    printf("...outputting.")
  }

  st_matrix("ci", ci)
  st_matrix("perf", perf)
  st_matrix("n", n[1])
  st_matrix("pi", pi')
  st_matrix("alpha", alpha)
  st_matrix("summary", summary)
}

// Function for determining Agresti-Coull CI in a fixed design
real rowvector ci_fixed_agresti_coull(real scalar s, real scalar m,
                                      real scalar alpha)
{
  if (s == 0) {
    m_tilde  = m + invnormal(alpha)^2
    pi_tilde = (s + 0.5*invnormal(alpha)^2)/m_tilde
    piL      = 0
    piU      = pi_tilde - invnormal(alpha)*sqrt(pi_tilde*(1 - pi_tilde)/m_tilde)
    if (piU < 0) {
      piU    = 0
    }
	else if (piU > 1) {
      piU    = 1
    }
  }
  else if (s == m) {
    m_tilde  = m + invnormal(alpha)^2
    pi_tilde = (s + 0.5*invnormal(alpha)^2)/m_tilde
    piL      = pi_tilde + invnormal(alpha)*sqrt(pi_tilde*(1 - pi_tilde)/m_tilde)
    piU      = 1
    if (piL < 0) {
      piL    = 0
    }
	else if (piL > 1) {
      piL    = 1
    }
  }
  else {
    m_tilde  = m + invnormal(alpha/2)^2
    pi_tilde = (s + 0.5*invnormal(alpha/2)^2)/m_tilde
    piL      = pi_tilde + invnormal(alpha/2)*sqrt(pi_tilde*(1 - pi_tilde)/m_tilde)
    piU      = pi_tilde - invnormal(alpha/2)*sqrt(pi_tilde*(1 - pi_tilde)/m_tilde)
    if (piL < 0) {
      piL    = 0
    }
	else if (piL > 1) {
      piL    = 1
    }
    if (piU < 0) {
      piU    = 0
    }
	else if (piU > 1) {
      piU    = 1
    }
  }
  return((piL, piU))
}

// Function for determining Clopper-Pearson CI in a fixed design
real rowvector ci_fixed_clopper_pearson(real scalar s, real scalar m,
                                        real scalar alpha)
{
  if (s == 0) {
    piL = 0
    piU = 1 - (alpha/2)^(1/m)
  }
  else if (s == m) {
    piL = (alpha/2)^(1/m)
    piU = 1
  }
  else {
    piL = invibeta(s, m - s + 1, alpha/2)
	piU = invibeta(s + 1, m - s, 1 - alpha/2)
  }
  return((piL, piU))
}

// Function for determining Jeffreys CI in a fixed design
real rowvector ci_fixed_jeffreys(real scalar s, real scalar m,
                                 real scalar alpha)
{
  if (s == 0) {
    piL = 0
    piU = invibeta(s + 0.5, m - s + 0.5, 1 - alpha)
  } else if (s == m) {
    piL = invibeta(s + 0.5, m - s + 0.5, alpha)
    piU = 1
  } else {
    piL = invibeta(s + 0.5, m - s + 0.5, alpha/2)
    piU = invibeta(s + 0.5, m - s + 0.5, 1 - alpha/2)
  }
  return((piL, piU))
}

real scalar lower_lim(real scalar pi, real scalar s, real scalar m,
                      real scalar alpha)
{
    return(0.5*binomialp(m, s, pi) + binomialtail(m, s + 1, pi) - alpha/2)
}

real scalar upper_lim(real scalar pi, real scalar s, real scalar m,
                      real scalar alpha)
{
    return(0.5*binomialp(m, s, pi) + binomial(m, s - 1, pi) - alpha/2)
}

// Function for determining mid-p CI in a fixed design
real rowvector ci_fixed_mid_p(real scalar s, real scalar m, real scalar alpha)
{
  if (s == 0) {
    piL   = 0
    piU   = 1 - alpha^(1/m)
  }
  else if (s == m) {
    piL   = alpha^(1/m)
    piU   = 1
  }
  else {
    output = mm_root(piL = ., &lower_lim(), 0, s/m, 0, 1000, s, m, alpha)
	if (output ~= 0) {
	  error("Root finding algorithm did not converge")
	}
	output = mm_root(piU = ., &upper_lim(), s/m, 1, 0, 1000, s, m, alpha)
    if (output ~= 0) {
	  error("Root finding algorithm did not converge")
	}
	if (piL < 0) {
      piL = 0
    }
	else if (piL > 1) {
      piL = 1
    }
    if (piU < 0) {
      piU = 0
    }
	else if (piU > 1) {
      piU = 1
    }
  }
  return((piL, piU))
}

// Function for determining Wald CI in a fixed design
real rowvector ci_fixed_wald(real scalar s, real scalar m, real scalar alpha)
{
  pi_hat  = s/m
  if (s == 0){
    piL   = 0
    piU   = 0
  } else if (s == m){
    piL   = 1
    piU   = 1
  } else {
    piL   = pi_hat - invnormal(1 - alpha/2)*sqrt(pi_hat*(1 - pi_hat)/m)
    piU   = pi_hat + invnormal(1 - alpha/2)*sqrt(pi_hat*(1 - pi_hat)/m)
    if (piL < 0) {
      piL = 0
    } else if (piL > 1) {
      piL = 1
    }
    if (piU < 0) {
      piU = 0
    } else if (piU > 1) {
      piU = 1
    }
  }
  return((piL, piU))
}

// Function for determining Wilson CI in a fixed design
real rowvector ci_fixed_wilson(real scalar s, real scalar m, real scalar alpha)
{
  pi_hat  = s/m
  if (s == 0){
    z     = invnormal(1 - alpha)
    z2    = z^2
    piL   = 0
    piU   = (pi_hat + z2/(2*m) + z*sqrt(pi_hat*(1 - pi_hat)/m + z2/(4*m^2)))/(1 + z2/m)
    if (piU < 0) {
      piU = 0
    } else if (piU > 1) {
      piU = 1
    }
  } else if (s == m){
    z     = invnormal(1 - alpha)
    z2    = invnormal(1 - alpha)^2
    piL   = (pi_hat + z2/(2*m) - z*sqrt(pi_hat*(1 - pi_hat)/m + z2/(4*m^2)))/(1 + z2/m)
    piU   = 1
    if (piL < 0) {
      piL = 0
    } else if (piL > 1) {
      piL = 1
    }
  } else {
    z     = -invnormal(alpha/2)
    z2    = z^2
    piL   = (pi_hat + z2/(2*m) - z*sqrt(pi_hat*(1 - pi_hat)/m + z2/(4*m^2)))/(1 + z2/m)
    piU   = (pi_hat + z2/(2*m) + z*sqrt(pi_hat*(1 - pi_hat)/m + z2/(4*m^2)))/(1 + z2/m)
    if (piL < 0) {
      piL = 0
    } else if (piL > 1) {
      piL = 1
    }
    if (piU < 0) {
      piU = 0
    } else if (piU > 1) {
      piU = 1
    }
  }
  return((piL, piU))
}

end