### Finding the MLE for the Hardy-Weiberg 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 } X11() hist(theta.mle.boot) X11() 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