# 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