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

[Douglas Bates]

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.