# SOLUTION TO STA 410/2102, SPRING 2002, ASSIGNMENT #3. # # Radford Neal, March 2002 # THE LOGISTIC FUNCTION. This is the "g" function in the assignment writeup. logistic <- function (z) 1/(1+exp(-z)) # LOG LIKELIHOOD FUNCTION. Takes as arguments the vector of values # for the explanatory variable, x, the vector of values for the # response variable, y (which are either -1 or +1), the alpha parameter, # and the vector of the (beta0, beta1) parameters. log.like <- function (x,y,alpha,beta) { g <- logistic (y * (beta[1]+beta[2]*x)) p <- alpha/2 + (1-alpha)*g sum(log(p)) } # FIND THE MLE USING THE EM ALGORITHM. Takes arguments as follows: # # x Vector of values for the explanatory variable. # y Vector of values for the response variable (-1 or +1). # iters Number of EM iterations to do (defaults to 100). # m.limit Limit on nlm iterations in the M step (default 1000000). # p.info Should information be printed at each iteration? (default T) # init.alpha Initial guess for alpha (defaults to 0.1) # init.beta Initial value for (beta0,beta1), defaults to the logistic # regression coefficients found by glm. # # If p.info is T, the parameter values and the log likelihood are printed # for each iteration. # # The value returned is a list containing elements alpha and beta, which # are the final estimate for the MLE. em <- function (x, y, iters=100, m.limit=100000, p.info=T, init.alpha=0.1, init.beta = coef(glm ((y+1)/2 ~ x, family=binomial))) { alpha <- init.alpha beta <- init.beta if (p.info) { cat("Iteration",0,":",round(log.like(x,y,alpha,beta),3), ":",alpha,beta,"\n"); } for (i in 1:iters) { # E Step. Finds r, the vector of probabilities that u is 1. p <- logistic(y*(beta[1]+beta[2]*x)) r <- (alpha/2) / (alpha/2 + (1-alpha)*p) # M Step. Finds new alpha from r, and new betas using nlm, starting # with the old betas as the initial guess. alpha <- mean(r) beta <- nlm ( function (beta) - sum((1-r)*log(logistic(y*(beta[1]+beta[2]*x)))), beta, iterlim=m.limit ) $ estimate if (p.info) { cat("Iteration",i,":",round(log.like(x,y,alpha,beta),3), ":",alpha,beta,"\n"); } } list (alpha=alpha, beta=beta) }