### Finding the MLE for the Hardy-Weinberg Law

x <- c(342,500,187)

log.lik <- function(the){
	(2*x[1] + x[2])*log(1-the) + (x[2] + 2*x[3])*log(the) + x[2]*log(2)
}

theta <- seq(0,1,0.01)

zz <- log.lik(theta)

# plot the log liklihood function of theta

plot(theta,zz,type="l")

# find the maximum of the log likelihood

theta.mle <- optimize(log.lik, c(0,1), maximum=TRUE)
theta.mle <- theta.mle$maximum

# So we have ...

theta.mle


# plot the log liklihood function of theta (Thanks to Michael Bissell)

plot(theta,zz,type="l", ylab='log likelihood of theta', main='log likelihood of theta')
lines(c(theta.mle,theta.mle),c(log.lik(theta.mle),min(zz[zz!=-Inf])),lty='dotted')
abline(h=min(zz[zz!=-Inf]))
points(c(theta.mle,theta.mle),c(log.lik(theta.mle),min(zz[zz!=-Inf])),cex=2)
title(sub=paste('y = max log lik theta = ', round(log.lik(theta.mle),4), ' at x = MLE =', round(theta.mle,4)))

p1.mle <- (1-theta.mle)^2
p1.mle
p2.mle <- 2*theta.mle*(1-theta.mle)
p2.mle
p3.mle <- theta.mle^2
p3.mle

# parametric bootstrap 

B <- 1000
N <- 1029
p <- c(p1.mle, p2.mle, p3.mle)

theta.mle.boot <- numeric(B)

for(i in 1:B){
	x <- rmultinom(1,N,p)
	theta.mle.star <- optimize(log.lik, c(0,1), maximum=TRUE)
	theta.mle.boot[i] <- theta.mle.star$maximum
}

hist(theta.mle.boot)

plot(density(theta.mle.boot, n=50, window="g"),type="l",xlab="x",ylab="density")

theta.mle.boot.m <- mean(theta.mle.boot)
theta.mle.boot.s <- sd(theta.mle.boot)
theta.mle.boot.ci <- quantile(theta.mle.boot,c(0.025,0.975))

p1.mle.boot <- (1-theta.mle.boot)^2
p2.mle.boot <- 2*theta.mle.boot*(1-theta.mle.boot)
p3.mle.boot <- theta.mle.boot^2

par(mfrow=c(2,2))
plot(density(p1.mle.boot, n=50, window="g"),type="l",main="p1",xlab="x",ylab="density")
plot(density(p2.mle.boot, n=50, window="g"),type="l",main="p2",xlab="x",ylab="density")
plot(density(p3.mle.boot, n=50, window="g"),type="l",main="p3",xlab="x",ylab="density")

p1.mle.boot.m <- mean(p1.mle.boot)
p1.mle.boot.m
p1.mle.boot.s <- sd(p1.mle.boot)
p1.mle.boot.s
p1.mle.boot.ci <- quantile(p1.mle.boot,c(0.025,0.975))
p1.mle.boot.ci

p2.mle.boot.m <- mean(p2.mle.boot)
p2.mle.boot.m
p2.mle.boot.s <- sd(p2.mle.boot)
p2.mle.boot.s
p2.mle.boot.ci <- quantile(p2.mle.boot,c(0.025,0.975))
p2.mle.boot.ci

p3.mle.boot.m <- mean(p3.mle.boot)
p3.mle.boot.m
p3.mle.boot.s <- sd(p3.mle.boot)
p3.mle.boot.s
p3.mle.boot.ci <- quantile(p3.mle.boot,c(0.025,0.975))
p3.mle.boot.ci