##someFunctions
##binomial likelihood normalized
like <- function(p, x, n, uselog = FALSE){
  ll <- x * log(p) + (n - x) * log(1 - p) - x * log(x / n) - (n - x) * log(1 - x / n)
  if (uselog) return(ll)
  else return(exp(ll))
}
##binomial likelihood plotting function
lPlot <- function(x, n, plotPoints = 1000, ...){
  pSeq <- seq(0, 1, length = plotPoints)
  lSeq <- like(pSeq, x, n)
  plot(pSeq, lSeq, type = "l", xlab = "p", ylab = "likelihood", ...)
  vals <- range(pSeq[lSeq >= 1 / 8])
  lines(vals, c(1/8,1/8))
  vals <- range(pSeq[lSeq >= 1 / 32])
  lines(vals, c(1/32,1/32))
}

##proble 4.b
noResamples <- 1000
y <- c(4,2,4,7,1,5,3,2,2,4,5,2,5,3,1,4,3,1,1,3)
Means <- apply(matrix(sample(y, length(y) * noResamples, replace = TRUE), noResamples), 1, mean)
quantile(Means, c(.025, .975))

##problem 5.a
lPlot(12, 20)

##problem 6.c
lPlot(12, 85)
##restrict the range
lPlot(12, 85, xlim = c(0, .5))

##problem 7.
lPlot(20, 1000)
##restrict the range
lPlot(20, 1000, plotPoints = 10000, xlim = c(0, .1))

##problem 8
##a useful function
ciFunc <- function(x, n, type = "wald", alpha = .05){
  if (type == "wald") phat <- x / n
  else if (type == "score") phat <- (x + 2) / (n + 4)
  phat + c(-1, +1) * qnorm(1 - alpha / 2) * sqrt(phat * (1 - phat) / n)
}
##a second useful function
mcCoverage <- function(p, n, type){
  xs <- rbinom(1000, n, p)
  cis <- t(sapply(xs, ciFunc, n = n, type = type))
  mean(cis[,1] < p & cis[,2] > p)
}

##a. and c.
mcCoverage(.3, 10, "wald")

##d.
mcCoverage(.3, 10, "score")

##e.
mcCoverage(.1, 10, "wald")
mcCoverage(.1, 10, "score")
mcCoverage(.5, 10, "wald")
mcCoverage(.5, 10, "score")

##extra credit
pSeq <- seq(0, 1, length = 100)
coverageWald <- sapply(pSeq, mcCoverage, n = 10, type = "wald")
coverageScore <- sapply(pSeq, mcCoverage, n = 10, type = "score")

plot(c(0, 1), c(0, 1), type = "n", xlab = "P", ylab = "coverage")
lines(pSeq, coverageWald, type = "l", lty = 1)
lines(pSeq, coverageScore, type = "l", lty = 2)
abline(h = .95)
legend(.4, .4, c("wald", "approx score"), lty = 1 : 2)