Bill: I didn't look at the code but I think the OP means that during the nlminb algorithm, the variance covariance parameters hit values such that the covariance matrix estimate becomes negative definite.
Again, I didn't look at the details but one way to deal with this is to have the likelihood function return -Inf whenever the covariance matrix becomes not positive definite. so, the likelihood should check for positive definiteness first before it actually calculates anything. If PD is not true, the -Inf value should push nlminb towards values that obtain a positive definite matrix. But there could be something more subtle going on that I'm not understanding. I don't know even what algorithm nlminb is using ( probably quasi-newton ) but this is one thing the OP could try. On Sun, Oct 20, 2013 at 9:41 PM, William Dunlap <wdun...@tibco.com> wrote: > Do you mean that your objective function (given to nlminb) parameterized > a positive definite matrix, P, as the elements in its upper (or lower) > triangle? > If so, you could reparameterize it by the non-zero (upper triangular) > elements > of the Choleski decomposition, C, of P. Compute P as crossprod(C), compute > the initial estimate of C as chol(P). > > Bill Dunlap > Spotfire, TIBCO Software > wdunlap tibco.com > > > > -----Original Message----- > > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] > On Behalf > > Of Steven LeBlanc > > Sent: Sunday, October 20, 2013 3:02 PM > > To: r-help@R-project.org > > Subject: [R] nlminb() - how do I constrain the parameter vector properly? > > > > Greets, > > > > I'm trying to use nlminb() to estimate the parameters of a bivariate > normal sample and > > during one of the iterations it passes a parameter vector to the > likelihood function > > resulting in an invalid covariance matrix that causes dmvnorm() to throw > an error. Thus, > > it seems I need to somehow communicate to nlminb() that the final three > parameters in > > my parameter vector are used to construct a positive semi-definite > matrix, but I can't see > > how to achieve this using the constraint mechanism provided. Additional > details are > > provided in the code below. > > > > Suggestions? > > > > Best Regards, > > Steven > > > > Generate the data set I'm using: > > > > mu<-c(1,5) > > sigma<-c(1,2,2,6) > > dim(sigma)<-c(2,2) > > set.seed(83165026) > > sample.full<-sample.deleted<-mvrnorm(50,mu,sigma) > > delete.one<-c(1,5,13,20) > > delete.two<-c(3,11,17,31,40,41) > > sample.deleted[delete.one,1]<-NA > > sample.deleted[delete.two,2]<-NA > > missing<-c(delete.one,delete.two) > > complete<-sample.deleted[!(is.na(sample.deleted[,1]) | is.na > (sample.deleted[,2])),] > > deleted<-sample.deleted[missing,] > > > > Try to estimate the parameters of the data set less the deleted values: > > > > exact<-function(theta,complete,deleted){ > > one.only<-deleted[!(is.na(deleted[,1])),1] > > two.only<-deleted[!(is.na(deleted[,2])),2] > > u<-c(theta[1],theta[2]) > > sigma<-c(theta[3],theta[5],theta[5],theta[4]) > > dim(sigma)<-c(2,2) > > -sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))- > > sum(log(dnorm(one.only,u[1],sigma[1,1])))- > > sum(log(dnorm(two.only,u[2],sigma[2,2]))) > > } > > > nlminb(start=c(0,0,1,1,0),objective=exact,complete=complete,deleted=deleted,control=l > > ist(trace=1)) > > > > Escape and it stops at Iteration 9 on my machine. > > ______________________________________________ > > 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. > > ______________________________________________ > 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. > [[alternative HTML version deleted]] ______________________________________________ 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.
