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.
"ay" "lag" "yp" "yc"
"1" 1 1 1376.384 33.81
"2" 1 2 1211.168 45.318
"3" 1 3 535.883 46.549
"4" 1 4 313.79 35.206
"5" 1 5 168.142 23.36
"6" 1 6 79.972 12.502
"7" 1 7 39.235 6.602
"8" 1 8 15.03 3.373
"9" 1 9 10.865 2.373
"10" 1 10 4.086 0.778
"11" 2 1 1576.278 37.663
"12" 2 2 1437.15 51.771
"13" 2 3 652.445 40.998
"14" 2 4 342.694 29.496
"15" 2 5 188.799 12.669
"16" 2 6 76.956 11.204
"17" 2 7 35.042 5.785
"18" 2 8 17.089 4.22
"19" 2 9 12.507 1.91
"20" 3 1 1763.277 40.63
"21" 3 2 1540.231 56.318
"22" 3 3 678.959 56.182
"23" 3 4 364.199 32.473
"24" 3 5 177.108 15.828
"25" 3 6 78.169 8.409
"26" 3 7 47.391 7.12
"27" 3 8 25.288 1.125
"28" 4 1 1779.698 40.475
"29" 4 2 1498.531 49.697
"30" 4 3 661.401 39.313
"31" 4 4 321.434 24.044
"32" 4 5 162.578 13.156
"33" 4 6 84.581 12.595
"34" 4 7 53.449 2.908
"35" 5 1 1843.224 37.127
"36" 5 2 1573.604 50.983
"37" 5 3 613.095 34.154
"38" 5 4 299.473 25.455
"39" 5 5 176.842 19.421
"40" 5 6 106.296 5.728
"41" 6 1 1962.385 41.125
"42" 6 2 1520.298 53.302
"43" 6 3 581.932 40.289
"44" 6 4 347.434 39.912
"45" 6 5 238.375 6.65
"46" 7 1 2033.371 57.515
"47" 7 2 1430.541 67.881
"48" 7 3 633.5 86.734
"49" 7 4 432.257 18.109
"50" 8 1 2072.061 61.553
"51" 8 2 1458.541 132.208
"52" 8 3 727.098 20.923
"53" 9 1 2210.754 112.103
"54" 9 2 1517.501 33.25
"55" 10 1 2206.886 37.554
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
