From: Douglas Bates <bates_at_stat.wisc.edu>

Date: Tue 31 May 2005 - 23:49:21 EST

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 Wed Jun 01 00:00:42 2005

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

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jun 01 00:00:42 2005

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