[R] Solved: linear regression example using MLE using optim()
Douglas Bates
bates at stat.wisc.edu
Tue May 31 15:49:21 CEST 2005
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.
More information about the R-help
mailing list