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

From: Ajay Narottam Shah <ajayshah_at_mayin.org>
Date: Tue 31 May 2005 - 19:17:09 EST


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.
-- 
Ajay Shah                                                   Consultant
ajayshah@mayin.org                      Department of Economic Affairs
http://www.mayin.org/ajayshah           Ministry of Finance, New Delhi

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Tue May 31 19:20:47 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:16 EST