[R] Question with nlm

From: <adschai_at_optonline.net>
Date: Fri, 15 Jun 2007 01:41:41 +0000 (GMT)


Hi,

I would really appreciate if I could get some help here. I'm using nlm to minimize my negative log likelihood function. What I did is as follows:

My log likelihood function (it returns negative log likelihood) with 'gradient' attribute defined inside as follows:

# ==========Method definition======================
logLikFunc3 <- function(sigma, object, totalTime) {

    y <- as.matrix(object_at_data$output[1:totalTime,1]);     x <- as.matrix(object_at_data$input[1:totalTime,]);
# compute necessary matrices

    M <- as.matrix(object_at_model$M);
    P <- diag(sigma*sigma);
    A <- AMatrix(totalTime, M, object_at_data$input[1:totalTime,]);
    Q <- IMatrix(totalTime)+A %*% outerM(IMatrix(totalTime-1),P) %*% t(A);
    invQ <- solve(Q,IMatrix(dim(Q)[1]));     xM <- matrix(rep(0, dim(M)[2]*totalTime), ncol=dim(M)[2], nrow=totalTime);     for (i in 1:totalTime) {

       xM[i,] <- x[i,] %*% powerM(M, -totalTime+i);     }
    tmp <- solve((t(xM) %*% invQ %*% xM), IMatrix(dim(xM)[2]));     Bt <- (tmp %*% t(xM)) %*% (invQ %*% y);     N <- IMatrix(totalTime)-(xM %*% tmp %*% t(xM) %*% invQ);     

    sigma2 <- (1/totalTime) * t(y- xM %*% Bt)%*% invQ %*% (y- xM %*% Bt);
# log likelihood function

    loglik <- -0.5*log(abs(det(diag(rep(sigma2,totalTime)))))-0.5*log(abs(det(Q)))-

       (0.5/sigma2)* (t(y- (xM%*% Bt)) %*% invQ %*% (y-(xM %*% Bt)));

    sgm <- sigma;
# gradients eq. (4.16)

    gr <- function(sgm) {

       gradVecs <- c();
       # sgm <- c(sigma1, sigma2);
       sgm <- sgm*sgm;
       for (i in 1:length(sgm)) {
          Eij <- matrix(rep(0, length(sgm)^2), nrow=length(sgm), ncol=length(sgm));
          Eij[i,i] <- 1.0;
          # trace term
          term1 <- -sum(diag((invQ %*% A) %*% outerM(IMatrix(totalTime-1),Eij) %*% t(A)));
          # very long term
          term2 <- (1/totalTime)*solve((t(y) %*% t(N) %*% invQ %*% y), IMatrix(dim(y)[2]));
          term3 <- (t(y) %*% t(N) %*% invQ %*% A) %*% outerM(IMatrix(totalTime-1),Eij) %*% (t(A) %*% invQ %*% N %*% y);
          gradVecs <- -1*c(gradVecs, term1+ (term2 %*% term3));
       } # end for
       print(paste("Gradient has length:", length(gradVecs)));
       return(gradVecs);

    }
    res <- -loglik;
    attr(res, "gradient") <- gradVecs;
    return(res);
}
#=========end method definition=====================================

Then when I call the nlm on this function, i.e.

nlm(f=logLikFunc3, p=as.numeric(c(1,1)), object=this, totalTime=200, print.level=2)

It complains that my analytic gradient returns vector length different from number of my unknowns. In this case, I tried print the length of gradient vector that I returned (as you could see in the code). It has the same length as my input parameter vectors. Did I do anything wrong here?

Also, I would like to be able to put some constraints on this optimization as well. I tried constrOptim with:

ui <- diag(c(1,1));
ci <- matrix(rep(0,2), ncol=1, nrow=2);

using the same parameters passed to nlm above. constrOptim gives me an error that initial value is in infeasible region which I don't quite understand. As my constraints simply says that the two parameters must be greater than zero. My assigned initial values are both 1. So it should be ok. Any help would be really appreciated. Thank you.


R-help_at_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 and provide commented, minimal, self-contained, reproducible code. Received on Fri 15 Jun 2007 - 01:50:03 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 Fri 15 Jun 2007 - 02:37:36 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.