On Oct 20, 2013, at 3:01 PM, Steven LeBlanc wrote: > 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. >
I recently noticed that teh dmvnorn function in package mixtools doe not throw an error when given a non-positive definite varaince-covariance matrix. Whether it should be doing so might be up for debate. Peter Dalgaard seemed to think that any self respecting such function _should_ be throwing NaN's back at you. The code it was using is rather compact: > mixtools::dmvnorm function (y, mu = NULL, sigma = NULL) { exp(logdmvnorm(y, mu = mu, sigma = sigma)) } <environment: namespace:mixtools> > mixtools::logdmvnorm function (y, mu = NULL, sigma = NULL) { if (is.vector(y)) y <- matrix(y, nrow = 1) d <- ncol(y) if (!is.null(mu)) y <- sweep(y, 2, mu, "-") if (is.null(sigma)) sigma <- diag(d) k <- d * 1.83787706640935 a <- qr(sigma) logdet <- sum(log(abs(diag(a$qr)))) if (nrow(y) == 1) mahaldist <- as.vector(y %*% qr.solve(a, t(y))) else mahaldist <- rowSums((y %*% qr.solve(a)) * y) -0.5 * (mahaldist + logdet + k) } I noticed that you were logging the dmvnorm values and so I wondered also if going directly to that function might save some time. (Note that I am way out of my mathematical dept here. Caveat emptor, maxima caveat emptor). -- David. > 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=list(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. David Winsemius Alameda, CA, USA ______________________________________________ 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.