# 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" )