# Bernouli Probability Mass Function (pmf or pdf)
set.seed(20131014)
n = 40
x = rbinom( 1, n, 0.20)
x
pdf = function(x,n,p){ choose(n,x)*p^x*(1-p)^(n-x) }
pdf(x,n,0.20)

# The likelihood fixes the data
# and makes the parameter the variable of interest
likelihood = pdf
likelihood(x,n,0.20)
likelihood(x,n,0.10)
likelihood(x,n,0.05)

# try a bunch of values for p
delta = 0.001
ps = seq(delta, 1-delta, delta) 
ps

ls = likelihood(x,n,ps)
ls
round(ls,3)

plot( ps, ls, typ='l', lwd=3 )


# For reasons we'll learn these plots are more
# interpretable when scaled by the mle's likelihood
mls = max(ls)
ls = ls/mls
plot( ps, ls, typ='l', lwd=3 )


# Exact Interval
library(Hmisc)
binconf(x,n,method="exact")

lbe = binconf(x,n,method="exact")[2]
ube = binconf(x,n,method="exact")[3]
lines(
  c(lbe, ube),
  c(likelihood(x,n,lbe)/mls,likelihood(x,n,ube)/mls),
  lwd=3, col="red"
)

# Wilson Interval
binconf(x,n)

lb = binconf(x,n)[2]
ub = binconf(x,n)[3]
lines(
  c(lb,ub),
  c(likelihood(x,n,lb)/mls,likelihood(x,n,ub)/mls),
  lwd=3, col="green"
)

# Wald Interval
binconf(x,n,method="asymptotic")

lb = binconf(x,n,method="asymptotic")[2]
ub = binconf(x,n,method="asymptotic")[3]
lines(
  c(lb,ub),
  c(likelihood(x,n,lb)/mls,likelihood(x,n,ub)/mls),
  lwd=3, col="purple"
)

# 1/8th Interval
lowerIndices <- min(which(ls >= 1/8))
upperIndices <- max(which(ls >= 1/8))
lb <- ps[lowerIndices]
ub <- ps[upperIndices]
lines(
  c(lb, ub),
  c(1/8,1/8),
  lwd=3, col="blue"
)