# STA 2102/450 Spring 2000 - S Functions for Assignment #2, Part 2 - R. M. Neal
# LOG LIKELIHOOD FUNCTION FOR The DECAY MODEL. The arguments are a vector
# of times for observations, a vector of counts at those times, and the
# two parameters, a and b. Constant factors are omitted from the likelihood.
decay.log.likelihood <- function(ti,ni,a,b)
{
sum (ni*(a-b*ti) - exp(a-b*ti))
}
# CALCULATE THE GRADIENT OF THE LOG LIKELIHOOD. Arguments are as above.
decay.gradient <- function(ti,ni,a,b)
{
c( sum (ni - exp(a-b*ti)), sum (-ni*ti + ti*exp(a-b*ti)) )
}
# CALCULATE THE HESSIAN OF THE LOG LIKELIHOOD. Arguments are as above.
decay.hessian <- function(ti,ni,a,b)
{
H <- matrix(0,2,2)
H[1,1] <- - sum (exp(a-b*ti))
H[2,2] <- - sum ((ti^2)*exp(a-b*ti))
H[1,2] <- H[2,1] <- sum (ti*exp(a-b*ti))
H
}
# FIND THE MLE BY THE NEWTON-RAPHSON METHOD. The arguments are a vector
# of times for observations, a vector of counts at those times, initial
# guesses for the two parameters, a and b, and the maximum number of iterations
# of Newton-Raphson iteration to do. The initial guesses default to the
# log of the maximum count in the data for a, and zero for b. The default
# for the maximum number of iterations is 100. A final "debug" argument can
# also be given, which if set to T causes intermediate values to be printed.
#
# The result is a list containing the maximum likelihood values for a and b,
# the covariance matrix for the estimates (found from the observed information),
# and the number of iterations required for convergence.
decay.mle <- function (ti,ni,a=log(max(ni)),b=0,max.iter=100,debug=F)
{
converged <- F
n <- 0
if (debug) cat(n,a,b,decay.log.likelihood(ti,ni,a,b),"\n")
while (!converged && n