# Re: [R] Solved: linear regression example using MLE using optim()

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Tue 31 May 2005 - 23:49:21 EST

Ajay Narottam Shah wrote:
> Thanks to Gabor for setting me right. My code is as follows. I found
> it useful for learning optim(), and you might find it similarly
> useful. I will be most grateful if you can guide me on how to do this
> better. Should one be using optim() or stats4::mle?
>
> set.seed(101) # For replicability
>
> # Setup problem
> X <- cbind(1, runif(100))
> theta.true <- c(2,3,1)

> y <- X %*% theta.true[1:2] + rnorm(100)
>
> # OLS --
> d <- summary(lm(y ~ X[,2]))
> theta.ols <- c(d\$coefficients[,1], d\$sigma)
>
> # Switch to log sigma as the free parameter
> theta.true[3] <- log(theta.true[3])
> theta.ols[3] <- log(theta.ols[3])
>
> # OLS likelihood function --
> ols.lf <- function(theta, K, y, X) {
> beta <- theta[1:K]
> sigma <- exp(theta[K+1])
> e <- (y - X%*%beta)/sigma
> logl <- sum(log(dnorm(e)/sigma))
> return(logl)
> }
>
> # Experiment with the LF --
> cat("Evaluating LogL at stupid theta : ", ols.lf(c(1,2,1), 2, y, X), "\n")
> cat("Evaluating LogL at true params : ", ols.lf(theta.true, 2, y, X), "\n")
> cat("Evaluating LogL at OLS estimates: ", ols.lf(theta.ols, 2, y, X), "\n")
>
> optim(c(1,2,3), # Starting values
> ols.lf, # Likelihood function
> control=list(trace=1, fnscale=-1), # See ?optim for all controls
> K=2, y=y, X=X # "..." stuff into ols.lf()
> )
> # He will use numerical derivatives by default.
>

It looks as if you are dividing e by sigma twice in your ols.lf function. Did you intend that? Also when you find yourself writing log(d<dist>(...)) it is usually better to write d<dist>(..., log = TRUE). Finally, you can avoid defining and passing K if you make log(sigma) the first element of the parameter vector theta. Combining all these would give you a likelihood function of the form

ols.lf <- function(theta, y, X) {

sum(dnorm((y - X %*% theta[-1])/exp(theta[1]), log = TRUE))) }

or, perhaps,

ols.lf <- function(theta, y, X) {

sum(dnorm(y, mean = X %*% theta[-1], sd = exp(theta[1]), log = TRUE)) }

The use of log = TRUE in the dnorm function is more than cosmetic. Frequently it is easier to evaluate the log-density than to evaluate the density. Also the log-density can be evaluated accurately over a wider range than can the density followed by the logarithm.

R-help@stat.math.ethz.ch mailing list