#############################################################
#
# Newton's method example for MLE
#
# Professor M. Zhu
# Current version  = Oct 2023
#
#############################################################
rm(list=ls(all=TRUE))                # clear workspace

myMLEfunc <- function(v, y, theta, max.iter=100, eps=1e-10)
{
# finds MLE for the censored Poisson model,
# given data (v,y), 
# using Newton's method ... 
# written by M. Zhu for teaching illustration
  i = 0; convg = FALSE; newtheta = theta
  while (i < max.iter && !convg)
  {
    theta = newtheta
    print(theta)                     # take a look
    u = exp(-v*theta)
    g = sum(y*v*u/(1-u) - (1-y)*v)   # first deriv
    H = -sum(y*v*v*u/(1-u)^2)        # second deriv
    newtheta = theta - g/H
    convg = abs(theta-newtheta)<eps
    i = i + 1
  }
  return(theta)
}

myLOGLIKfunc <- function(theta, v, y)
{
# computes the log-likelihood function for the
# censored Poisson model, given data (v,y) ...
# written by M. Zhu for teaching illustration
  u = exp(-v*theta)
  return(sum(y*log(1-u)-(1-y)*v*theta))
}

# Here are some data (Essential Statistics, Table 5.1)
v=c(8,8,8,8,8,8,8,8,8,8,
    4,4,4,4,4,4,4,4,4,4,
    2,2,2,2,2,2,2,2,2,2,
    1,1,1,1,1,1,1,1,1,1)
y=c(1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,0,0,
    1,1,1,1,1,1,1,0,0,0,
    1,1,1,0,0,0,0,0,0,0)

# calculate and plot the loglikelihood function for
# a number of different theta's ...
theta = seq(0.25,0.9,len=50)
plot(theta, sapply(theta, myLOGLIKfunc, v, y),
  type="l", 
  xlab=expression(theta),ylab='',
  main=expression(paste(log,group("[",L(theta),"]"))),
  cex.main=2,cex.lab=1.5,cex.axis=1.25)

# starting from theta=0.1
theta.opt = myMLEfunc(v,y,theta=0.1)

# now add the MLE onto the plot
abline(v=theta.opt, col='red')