library(MASS) # for truehist

# function to calculate binomial log likelihood	
# p: prob of success
# n: number of trials
# x: successes observed
binomial.loglik <- function(p, n, x) {
	# catch illegal tries, returning very small likelihood
	if (p <= 0 || p >= 1)
		return (-9999999999);

	return(x * log(p)  + (n - x) * log(1 - p));
}

# A function to fit an MLE to an observed binomial experiment.  Pass
# in the number of trials and number of successes observed.
binom.mle <- function(N, successes) {
	mle <- optim(c(.5), # starting value for p
		binomial.loglik,  # the function to maximize
		method="BFGS",    # optimization algorithm
		hessian=T,        # we want the hessian
		control=list(fnscale=-1), # maximize, don't minimize
		x=successes, n=N) # pass in observed data
	return(mle)
}

# Simulate binomial data with N trials and probability of success p
N <- 100
p <- .65

data <- runif(N) < p   # create a vector of weighted coin flips
successes <- sum(data) # count the number of heads
cat("Successes\n")
print(successes)

# Plot the log-likelihood and compare it to p
ruler <- seq(0.01, .99, .01)
plot(ruler, binomial.loglik(ruler, x=successes, n=N), type='l',
		 ylab="logL(p)", xlab="p")
abline(v = p, col="red")

# Fit the mle, using the function above
mle <- binom.mle(N, successes)

# Print MLE estimates
cat("MLE of p\n")
print(mle$par)
cat("Hessian at MLE\n")
print(mle$hessian)
cat("MLE Standard Error\n")
SE <- sqrt(-1 * (1 / mle$hessian))
print(SE)
cat("MLE CI\n")
CI <- c(mle$par - 1.96 * SE, mle$par + 1.96 * SE)
print(CI)

# Now a non-parametric bootstrap for the parameter p
NBOOTS <- 1000

# Iterate NBOOTS times, resampling and fitting each time
p.boots <- double(NBOOTS) # set up vector for results
for (i in 1:NBOOTS) {
	data.star <- sample(data, N, replace =T)
	successes.star <- sum(data.star)
	mle.star <- binom.mle(N, successes.star)
	p.boots[i] <- mle.star$par
}

# Plot the bootstrap distribution of p
X11()  # opens a new plot window
truehist(p.boots)
# Plot a line marking the mle estimate
abline(v = mle$par, col="red")
# Calcualte 95% CIs, print, and plot
CI.boots <- quantile(p.boots, probs=c(.025, .975))
# Plot the bootstrap CIs and print them to the terminal
abline(v = CI.boots, col="red", lty=2)
cat("Non-parametric bootstrap CIs\n")
print(CI.boots)

# Now plot traditional CIs
abline(v = CI, col="blue", lty=2)