[R] Help on simple problem with optim

[William Dunlap]

You can record all arguments and return values of the
calls that optim(par,fn) makes to fn with a function
like the following.  It takes your function and makes
a new function that returns the same thing but also
records information it its environment.  Thus, after
optim is done you can see its path to the optimum.

  trackFn <- function (fn) 
  {
      X <- NULL
      VALUE <- NULL
      force(fn)
      function(x) {
          X <<- rbind(X, x)
          val <- fn(x)
          VALUE <<- c(VALUE, val)
          val
      }
  }

E.g.,
> fn <- function(x)sum( (sin(x)-c(-1,0,1)) ^ 2 )
> optim(c(0,0,0), tfn <- trackFn(fn))
$par
[1] -1.583466e+00  6.726235e-05  1.558809e+00

$value
[1] 1.612853e-08

$counts
function gradient 
     146       NA 

$convergence
[1] 0

$message
NULL

> with(environment(tfn), length(VALUE)) # agrees w/ counts above
[1] 146
> with(environment(tfn), plot(VALUE, log="y"))
> with(environment(tfn), matplot(X, VALUE, log="y", type="l"))

You could also record warning messages by including a call
to withCallingHandlers() that stashed the warning in a list.

Recall trackFn each time you call optim().

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of Zhang,Yanwei
> Sent: Thursday, September 09, 2010 11:54 AM
> To: r-help at r-project.org
> Subject: [R] Help on simple problem with optim
> 
> Dear all,
> 
> I ran into problems with the function "optim" when I tried to 
> do an mle estimation of a simple lognormal regression. Some 
> warning message poped up saying NANs have been produced in 
> the optimization process. But I could not figure out which 
> part of my code has caused this. I wonder if anybody would 
> help. The code is in the following and the data is in the attachment.
> 
> 
> da <- read.table("da.txt",header=TRUE)
> 
> # fit with linear regression using log transformation of the 
> response variable
> fit <- lm(log(yp) ~ as.factor(ay)+as.factor(lag),data=da)
> 
> # define the log likelihood to be maximized over
> llk.mar <- function(parm,y,x){
>         # parm is the vector of parameters
>         # the last element is sigma
>         # y is the response
>         # x is the design matrix
>         l <- length(parm)
>         beta <- parm[-l]
>         sigma <- parm[l]
>         x <- as.matrix(x)
>         mu <- x %*% beta
>         llk <- sum(dnorm(y, mu, sigma,log=TRUE))
>         return(llk)
> }
> 
> # initial values
> parm <- c(as.vector(coef(fit)),summary(fit)$sigma)
> y <- log(da$yp)
> x <- model.matrix(fit)
> 
> op <- optim(parm, llk.mar, 
> y=y,x=x,control=list(fnscale=-1,maxit=100000))
> 
> 
> After running the above code, I got the warning message:
> Warning messages:
> 1: In dnorm(x, mean, sd, log) : NaNs produced
> 2: In dnorm(x, mean, sd, log) : NaNs produced
> 
> 
> I would really appreciate if anybody would help to point out 
> the problem with this code or tell me how to trace it down 
> (using "trace"?)?
> Many thanks in advance.
> 
> 
> 
> 
> 
> 
> 
> Wayne (Yanwei) Zhang
> Statistical Research
> CNA
> 
> 
> 
> 
> 
> NOTICE:  This e-mail message, including any attachments and 
> appended messages, is for the sole use of the intended 
> recipients and may contain confidential and legally 
> privileged information.
> If you are not the intended recipient, any review, 
> dissemination, distribution, copying, storage or other use of 
> all or any portion of this message is strictly prohibited.
> If you received this message in error, please immediately 
> notify the sender by reply e-mail and delete this message in 
> its entirety.
>