[R] Error in maximum likelihood estimation.

From: Dong-hyun Oh <r.arecibo_at_gmail.com>
Date: Mon, 16 Jun 2008 16:41:04 +0200


Dear UseRs,

I wrote the following function to use MLE.



mlog <- function(theta, nx = 1, nz = 1, dt){

   beta <- matrix(theta[1:(nx+1)], ncol = 1)    delta <- matrix(theta[(nx+2):(nx+nz+1)], ncol = 1)    sigma2 <- theta[nx+nz+2]
   gamma <- theta[nx+nz+3]

   y <- as.matrix(dt[, 1], ncol = 1)
   x <- as.matrix(data.frame(1, as.matrix(dt[, 2:(nx+1)], ncol = 2)))
   z <- as.matrix(dt[, (nx+2):(nx+nz+1)], ncol = nz)

   d <- z %*% delta / (gamma * sigma2)^.5    mustar <- (1-gamma) * z %*% delta - gamma * ( y - x %*% beta)    sigmastar <- (gamma * (1-gamma) * sigma2)^.5    dstar <- mustar / sigmastar

   loglik <- (-0.5 * nrow(x) *(log(2*pi) + log(sigma2))
              -0.5 * sum(( y - x %*% beta + z %*% delta)^2/sigma2)
              -sum(log(pnorm(d))) + sum(log(pnorm(dstar))))
   return(-loglik)
}

Loglikelihood function is from page 21of Battese and Coelli (1993). (You can download this article at http://www.une.edu.au/economics/publications/econometrics/emwp69.PDF   )

To test the above function with an artificial data set, I created the following data.frame.


x1 <- abs(rnorm(100))*100
x2 <- abs(rnorm(100))*10
z1 <- abs(rnorm(100))*5
z2 <- abs(rnorm(100))*7

y <- abs(0.3 + 0.3* log(x1) + 0.7* log(x2)) dat <- data.frame(log(y), log(x1), log(x2), z1, z2)

The starting value I set up is as follows:



theta.start <- c(0.09, 0.008, 0.008, 0.023, 0.0008, 0.008, 0.008)

Applying nlm() function to the above function with the starting values gives the following error message.



out <- nlm(mlog, theta.start, nx = 2 , nz = 2, dt = dat, print.level = 2, hessian\

iteration = 0
Step:
[1] 0 0 0 0 0 0 0
Parameter:
[1] 0.0900 0.0080 0.0080 0.0230 0.0008 0.0080 0.0080 Function Value
[1] 6116.2
Gradient:
[1] -10704.8 -46465.7 -22536.4 47168.2 54542.9 -776512.5 579.9

Error in nlm(mlog, theta.start, nx = 2, nz = 2, dt = dat, print.level = 2, :

   (converted from warning) NA/Inf replaced by maximum positive value


My questions are
1. Is my loglikelihood function set up appropriately? 2. How to set up starting values efficiently? I read a thread shown at https://stat.ethz.ch/pipermail/r-help/2005-September/079617.html   , but I cannot figure out how to set up starting values with expand.grid() of parameters I used (beta, delta, gamma, sigma2).

Thank you in advance.

Sincerely,
Dong-hyun Oh



R-help_at_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. Received on Mon 16 Jun 2008 - 15:02:30 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 16 Jun 2008 - 15:30:43 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive