# Help producing the table for the example in class23-prep-a ----
n = 12
alpha = .1

doPrintTables = FALSE
doDrawFigure = FALSE
doLargeSampleCI = TRUE

findNonrejectionRegion = function(p,alpha, theta)
{
  # Works in from the ends choosing the smallest value to accumulate into rejection region    
  favorRight = TRUE  #ties go the the right side 
  if (theta > .5)
  {
    favorRight = FALSE
  }
  n = length(p) - 1  #matches n in binom(n,theta)
  leftInd = 0
  rightInd = n
  totsig = 0.0
  while (totsig <= alpha)
  {
    a = p[leftInd+1]
    b = p[rightInd+1]
    if ((a < b) || (a==b && !favorRight))
    {
      if ((totsig + a) > alpha)
      {
        break
      }
      leftInd = leftInd + 1
      totsig = totsig + a
    }
    else 
    {
      if ((totsig + b) > alpha)
      {
        break
      }
      rightInd =rightInd - 1
      totsig = totsig + b
    }
    if (leftInd >= rightInd)
    {
      break  #This should never happen
    }
  }
  return(c(leftInd:rightInd)) #range indices are 0-based
}

printTable = function(n, thetas, alpha, nrowsToColor, highlightXValue = -1)
{
  # nRowsToColor also indicates the rows to print the significance
  cline = sprintf("\\cline{1-%d}\n",n+2)

  cat("% Created with class23-prep.r\n")
  cat("% Needs package colortbl\n")
  cat("\\def\\r{\\cellcolor{rejcolor}}\n")
  cat("\\def\\nr{\\cellcolor{nonrejcolor}}\n")
  cat("\\def\\nrh{\\cellcolor{nonrejhighlightcolor}}\n")
  cat("\\definecolor{rejcolor}{rgb}{1,.65,.30}\n")
  cat("\\definecolor{nonrejcolor}{rgb}{0,.64,.84}\n")
  cat("\\definecolor{nonrejhighlightcolor}{rgb}{.16,.8,1}\n")
  
  s = "\\small{\n\\begin{tabular}{|c|"
  for (j in 0:n)
  {
    s = sprintf("%sc|",s)
  }
  s = sprintf("%sc}\n",s)
  cat(s)
  cat(cline)

  s = "$\\theta$/$x$"
  for (j in 0:n)
  {
    s = sprintf("%s & %d", s, j)
  }
  s2 = ""
  if (nrowsToColor > 0)
  {
    s2 = "significance"
  }
  s = sprintf("%s & %s\\\\\n",s,s2)
  cat(s)
  cat(cline)

  for (k in 1:length(thetas))
  {
    theta = thetas[k]
    colorCells = FALSE # This also determines if significance is printed
    if (k <= nrowsToColor)
    {
      colorCells = TRUE
    }
    p = dbinom(0:n,n,theta) #0-based
    nonrej = findNonrejectionRegion(p,alpha,theta) #0-based
    significance = 1 - sum(p[nonrej+1]) #0-based
    s = sprintf("%.1f", theta)
    for (j in 0:n)
    {
      cellcolstart = ""
      cellcolend = ""
      if (colorCells)
      {
        cellcolstart = "\\r"
        if (any(nonrej == j))
        {
          if (highlightXValue == j)
          {
            cellcolstart = "\\nrh{\\bf"
            cellcolend = "}"
          }
          else
          {
            cellcolstart = "\\nr"
          }
        }
      }
      s = sprintf("%s & %s %.2f%s", s, cellcolstart, p[j+1], cellcolend)
    }
    s2 = ""
    if (colorCells)
    {
      s2 = sprintf("%.3f", significance) 
    }
    s = sprintf("%s & %s\\\\\n",s, s2)
    cat(s)
    cat(cline)
  }
  s = "\\end{tabular}\n}\n\n"
  cat(s)
}

if (doPrintTables)
{
  thetas = seq(1,0,-.1)
  printTable(n, thetas, alpha, 0)
  printTable(n, thetas, alpha, 1)
  printTable(n, thetas, alpha, 2)
  printTable(n, thetas, alpha, length(thetas), highlightXValue=8)
}

maxsig = 0.0
maxtheta = 0.0
drawfig = function(n,alpha, CIForX = -1)
{
  #CIForX = -1 is a flag not to draw this
  
  thetas = seq(0,1,.005)
  nonrejcol = rgb(0,.64,.84)
  nonrejhighlightcol = rgb(.16,.80,1)
  rejcol = rgb(1,.65,.30)
  lwd = 3
  drawverticallines = FALSE
  plot(c(-.5,n+.5),c(0,1), type='n', xlab='x', ylab=expression(theta))
  s = sprintf("n=%d, \u03B1=%.3f",n,alpha)
  title(s)
  ciStart = -1 #-1 and 100 are flags
  ciEnd = 100
  for (theta in thetas)
  {
    p = dbinom(0:n,n,theta) #0-based
    th = c(theta,theta)
    nonrej = findNonrejectionRegion(p,alpha,theta) #0-based
    sig = 1.0 - sum(p[nonrej+1])
    if (sig > maxsig)
    {
      maxsig = sig
      maxtheta = theta
    }
    m = length(nonrej)
    a = nonrej[1]-.5
    b = nonrej[length(nonrej)]+.5 
    lines(c(-.5,a), th, col=rejcol,lwd=lwd)
    lines(c(a,b), th, col=nonrejcol,lwd=lwd)
    lines(c(b,n+.5),th, col=rejcol,lwd=lwd)
    if (any(nonrej == CIForX))
    {
      if (ciStart == -1)
      {
        ciStart = theta
      }
    }
    else
    {
      if (ciStart != -1 && ciEnd == 100)
      {
        ciEnd = theta
      }
    }
  }
  if (drawverticallines)
  {
    for (j in 0:(n+1))
    {
      abline(v=j-.5)
    }
  }
  if (ciStart >= 0)
  {
    rect(CIForX-.5, ciStart,CIForX+.5, ciEnd, col=nonrejhighlightcol,
         border='blue')
    abline(h=ciStart,col='blue',lty='dashed')
    abline(h=ciEnd,col='blue',lty='dashed')
    print(c(ciStart,ciEnd, maxsig,maxtheta))
  }
}

if  (doDrawFigure)
{
  drawfig(1,.5)
  drawfig(1,.1)
  drawfig(2,.1)
  drawfig(10,.1)
  drawfig(50,.1)
  drawfig(100,.1)
  drawfig(12,.1,CIForX=8)
}

largeSampleCI = function(n,conf, ntrials, mu, rgenerator, useTCritVal)
{
  alpha = 1- conf
  if (useTCritVal)
  {
    # Note everything on the web indicates that people use z_alpha/2 for
    # large sample confidence intervals. The use to t_alpha/2 is just experimenting.
    c_alpha2 = -qt(alpha/2, n-1)
  }
  else
  {
    c_alpha2 = -qnorm(alpha/2,mean=0,sd=1) 
  }
  ncorrect = 0
  for (j in 1:ntrials)
  {
    #x = rexp(n,lambda)
    x = rgenerator(n)
    xbar = mean(x)
    sxbar = sqrt(var(x)/n)
    if (abs((xbar-mu)/sxbar) <= c_alpha2)
    {
      ncorrect = ncorrect + 1
    }
  }
  s = sprintf("c=%.3f, n=%d, ntrials=%d, correct=%.3f\n", 
              conf,n,ntrials, ncorrect/ntrials)
  cat(s)
}

if (doLargeSampleCI)
{
  ntrials = 10000

  lambda = 1
  mu = 1/lambda
  rgenerator = function(n) {return(rexp(n,lambda)) }
  cat("\nRunning with Exp(1) and z_alpha/2\n")
  for (n in c(20,50,100,400))
  {   for (conf in c(.95,.90,.80))
      { largeSampleCI(n,conf,ntrials, mu, rgenerator, FALSE)
      }
  }

  lambda = 1
  mu = 1/lambda
  rgenerator = function(n) {return(rexp(n,lambda)) }
  cat("\nRunning with Exp(1) and t_alpha/2\n")
  for (n in c(20,50,100,400))
  {   for (conf in c(.95,.90,.80))
      { largeSampleCI(n,conf,ntrials, mu, rgenerator, TRUE)
      }
  }
 

  mu = 0
  rgenerator = function(n) {return(rnorm(n)) }
  cat("\nRunning with norm(0,1) and z_alpha/2\n")
  for (n in c(20,50,100,400))
  {   for (conf in c(.95,.90,.80))
  { largeSampleCI(n,conf,ntrials, mu, rgenerator, FALSE)
  }
  }

  mu = 0
  rgenerator = function(n) {return(rnorm(n)) }
  cat("\nRunning with norm(0,1) and t_alpha/2\n")
  for (n in c(20,50,100,400))
  {   for (conf in c(.95,.90,.80))
  { largeSampleCI(n,conf,ntrials, mu, rgenerator, TRUE)
  }
  }
  
}