[R] Help on simple problem with optim

[Cristian Montes]

Did you check if the data in "da" has any NA in the dependent or the independent data?
Remember that your function llk.mar is going to evaluate dnorm for each pair.  If any of those
pairs has an NA value, your function will return an NA at the end (sum(c(NA,1,2,3)) = NA)

I would check if the llk.mar function is fine for the whole domain of your data.  My suggestion
is to add an if right before llk<-sum(dnorm....), where you evaluate for NAs in that vector.

If you find one, get rid of it before returning the function!.  Try that before optim, then
let it do the solving.

Cheers,

Cristián Montes.

-----Mensaje original-----
De: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] En nombre de Zhang,Yanwei
Enviado el: Jueves, 09 de Septiembre de 2010 02:54 p.m.
Para: r-help at r-project.org
Asunto: [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.