# Using the Newton-Raphson method to find the MLE for the Hardy-Weinberg Law

# the data

x <- c(342,500,187)

# log likelihood.  
# Note the log likelihood is proportional to the function.

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

# plot the log likelihood.

theta = seq(0,1,0.01)
zz = log.lik(theta)
plot(theta,zz,type = "l")

# the derivative of the log likelihood.  

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

# plot the derivative of the log likelihood.

zz = log.lik.d(theta)
plot(theta,zz,type = "l")

### Newton-Raphson

N = 10			# number of interations
theta = numeric(N)
l.lik.d = numeric(N)
l.lik.dd = numeric(N)
theta[1] = 0.4
for ( k in 2:N){
	l.lik.d[k-1] = -(2*x[1]+x[2])/(1-theta[k-1])+(x[2]+2*x[3])/theta[k-1]
	l.lik.dd[k-1] = -(2*x[1]+x[2])/(1-theta[k-1])^2 - (x[2]+2*x[3])/theta[k-1]^2
	theta[k] <- theta[k-1] - (l.lik.d[k-1])/(l.lik.dd[k-1])
}
l.lik.d
l.lik.dd
theta

# We would set an epsilon if we were writing a general algorithm.
# We can see that there is convergence for the NR Method here.