# Solution to STA 410/2102 Assignment #2, Spring 2003 - Program code. # COMPUTE THE LOG LIKELIHOOD FOR A LOGISTIC MODEL. Takes as arguments the # vector of data values, the mu parameter, and the omega parameter. The # omega parameter defaults to 1. logistic.log.lik <- function (x, mu, omega=1) { n <- length(x) v <- (x-mu)/omega - sum (log(1+exp(-v)) + log(1+exp(v))) - n*log(omega) } # COMPUTE PARTIAL DERIVATIVES OF THE LOG LIKELIHOOD. First and second # partial derivatives of the log likelihood with respect to the mu and # omega parameters are computed by these functions. All of them take # as arguments the vector of data values, the mu parameter, and the # omega parameter. The omega parameter defaults to 1 for the functions # that take derivatives only with respect to mu. logistic.log.lik.dmu <- function (x, mu, omega=1) { n <- length(x) v <- (x-mu)/omega (1/omega) * (n - sum (2/(1+exp(v)))) } logistic.log.lik.domega <- function (x, mu, omega) { n <- length(x) v <- (x-mu)/omega (1/omega) * (sum (v - 2*v/(1+exp(v))) - n) } logistic.log.lik.dmu2 <- function (x, mu, omega=1) { v <- (x-mu)/omega e <- exp(v) f <- 1/e - (1/omega^2) * sum (f/(1+f)^2 + e/(1+e)^2) } logistic.log.lik.domega2 <- function (x, mu, omega) { n <- length(x) v <- (x-mu)/omega e <- exp(v) f <- 1/e (1/omega^2) * (n - sum(v^2*f/(1+f)^2 + v^2*e/(1+e)^2 - 2*v/(1+e) + 2*v/(1+f))) } logistic.log.lik.dmu.domega <- function (x, mu, omega) { v <- (x-mu)/omega e <- exp(v) f <- 1/e (1/omega^2) * sum (- v*e/(1+e)^2 + 1/(1+e) - v*f/(1+f)^2 - 1/(1+f)) } # FIND APPROXIMATE DERIVATIVES NUMERICALLY. For use in testing the # derivative functions above. num.deriv <- function (f, a, delta=0.001) { (f(a+delta) - f(a-delta)) / (2*delta) } # FIND THE MAXIMUM LIKELIHOOD ESTIMATE FOR MU, WITH OMEGA FIXED AT 1. # # This functions solves Part I of the assignment. The arguments are # the vector of data values, the number of iterations of Newton-Raphson # to do, the intial value for mu, and the debug flag (default F). If # the debug flag is T, the values for mu and the log likelihood are printed # at each iteration. The default for the initial value of mu is the median # of the data. The value returned is a list with elements called "mle", # for the maximum likelihood estimate of mu, and "std.err", for the standard # error of the mle, based on the observed information. logistic.mle1 <- function (x, n, mu0=median(x), debug=F) { mu <- mu0 # Loop to do n Newton-Raphson iterations. i <- 0 repeat { if (debug) cat(i,": mu =",mu,"log.lik =",logistic.log.lik(x,mu),"\n") if (i>=n) break; # Do a Newton-Raphson update of the parameter mu. mu <- mu - logistic.log.lik.dmu(x,mu) / logistic.log.lik.dmu2(x,mu) i <- i+1 } list (mle=mu, std.err=sqrt(-1/logistic.log.lik.dmu2(x,mu))) } # FIND THE MAXIMUM LIKELIHOOD ESTIMATES FOR MU AND OMEGA. # # This functions solves Part II of the assignment. The arguments are # the vector of data values, the number of iterations of Newton-Raphson # to do, the intial values for mu and omega, and the debug flag (default F). # If the debug flag is T, the values for mu and omega and the log likelihood # are printed at each iteration. The default for the initial value of mu is # the median of the data. The default for the initial value of omega is # the interquartile range divided by log(9), which would produce the correct # value for omega if the interquartile range were that of the true distribution. # The value returned is a list with elements called "mle", for the maximum # likelihood estimate (a vector with the mles for mu and omega), and "std.err" # for the corresponding standard errors, based on the observed information. logistic.mle2 <- function (x, n, mu0=median(x), omega0=(quantile(x,0.75,names=F)-quantile(x,0.25,names=F))/log(9), debug=F) { mu <- mu0 omega <- omega0 H <- matrix(0,2,2) # Loop to do n Newton-Raphson iterations. i <- 0 repeat { if (debug) cat(i,": mu =",mu," omega =",omega, " log.lik =",logistic.log.lik(x,mu,omega),"\n") if (i>=n) break; # Do a Newton-Raphson update for the parameters mu and omega. These # parameters are regarded as a vector of length two. The first and # second partial derivatives of the log likelihood are put into the # gradient vector and the Hessian matrix, allowing the Newton-Raphson # update to be performed. The new values of mu and omega are then # extracted from the new vector. v <- c(mu,omega) g <- c(logistic.log.lik.dmu(x,mu,omega),logistic.log.lik.domega(x,mu,omega)) H[1,1] <- logistic.log.lik.dmu2(x,mu,omega) H[2,2] <- logistic.log.lik.domega2(x,mu,omega) H[1,2] <- H[2,1] <- logistic.log.lik.dmu.domega(x,mu,omega) v <- v - solve(H,g) mu <- v[1] omega <- v[2] i <- i+1 } # Create Hessian matrix one last time to find standard errors. H[1,1] <- logistic.log.lik.dmu2(x,mu,omega) H[2,2] <- logistic.log.lik.domega2(x,mu,omega) H[1,2] <- H[2,1] <- logistic.log.lik.dmu.domega(x,mu,omega) list (mle=c(mu,omega), std.err=sqrt(diag(solve(-H)))) }