# MAXIMUM LIKELIHOOD ESTIMATION FOR POISSON DATA WITH ZEROS UNOBSERVED.
# 
# Implements the three methods discussed in lecture.  For all three functions,
# the arguments are the vector of observations, the number of iterations to
# do, and the initial guess (defaulting to the sample mean).  The result is 
# a data frame with one row for each iteration (including the initial guess),
# with the columns being the estimate for lambda at that iteration and the
# log likelihood for that value of lambda (minus the factorial terms that
# don't involve lambda).


# COMPUTE LOG LIKELIHOOD.  The arguments are the data vector and a value for
# lambda (or vector of values).  The result is the log probability of the data 
# given that value for lambda, omitting the factorial terms (or a vector of
# log probabilities if lambda is a vector).

nzp.log.likelihood <- function (n, lambda)
{
  sum(n) * log(lambda) - length(n) * (lambda + log(1-exp(-lambda)))
}


# FIND MLE BY SIMPLE ITERATION.

nzp.simple.iteration <- function (n, r, lambda0=mean(n))
{ 
  mean.n <- mean(n)

  lambda <- rep(0,r+1)
  lambda[1] <- lambda0

  for (i in 1:r)
  { lambda[i+1] <- mean.n * (1 - exp(-lambda[i]))
  }

  data.frame (lambda=lambda, log.lik=nzp.log.likelihood(n,lambda))
}


# FIND MLE BY NEWTON-RAPHSON ITERATION.

nzp.newton.raphson <- function (n, r, lambda0=mean(n))
{ 
  mean.n <- mean(n)

  lambda <- rep(0,r+1)
  lambda[1] <- lambda0

  for (i in 1:r)
  { e <- exp(-lambda[i])
    lambda[i+1] <- lambda[i] - 
       (mean.n/lambda[i] - 1/(1-e)) / (e/(1-e)^2 - mean.n/lambda[i]^2)
  }

  data.frame (lambda=lambda, log.lik=nzp.log.likelihood(n,lambda))
}


# FIND MLE BY THE METHOD OF SCORING.

nzp.method.of.scoring <- function (n, r, lambda0=mean(n))
{ 
  mean.n <- mean(n)

  lambda <- rep(0,r+1)
  lambda[1] <- lambda0

  for (i in 1:r)
  { e <- exp(-lambda[i])
    lambda[i+1] <- lambda[i] - 
       (mean.n/lambda[i] - 1/(1-e)) / ((e/(1-e) - 1/lambda[i]) / (1-e))
  }

  data.frame (lambda=lambda, log.lik=nzp.log.likelihood(n,lambda))
}