Steven: I'm not sure if it makes a difference but you might want to start off with the "square root" of sigmaStart because that will really start you off with sigmaStart. Essentially, what Bill is doing is a reparameterization using the correlation and the 2 standard deviations so compute those based on the sigmaStart and use those as the starting values. That might be a little faster.
On Mon, Oct 21, 2013 at 12:28 PM, William Dunlap <wdun...@tibco.com> wrote: > [I added r-help to the cc again. Please keep the replies on the list as > there are others who can contribute to or learn from them.] > > > I also learned you have to be very careful with the starting value, as > the simple identity > > matrix becomes singular under the transformation. > > That is why I suggested using the non-zero entries in chol(sigmaStart), > where > sigmaStart is your initial estimate of the variance matrix, as the initial > values of theta[3:5]. > > I should have said that crossProd(x) gives you a positive semidefinite > matrix, > not positive definite. When I said that variances near 0 would cause > problems > I meant that a singular covariance matrix would cause problems. > > > Bill Dunlap > Spotfire, TIBCO Software > wdunlap tibco.com > > > > -----Original Message----- > > From: Steven LeBlanc [mailto:ores...@gmail.com] > > Sent: Monday, October 21, 2013 9:21 AM > > To: William Dunlap > > Subject: Re: [R] nlminb() - how do I constrain the parameter vector > properly? > > > > > > On Oct 21, 2013, at 7:52 AM, William Dunlap <wdun...@tibco.com> wrote: > > > > > Try defining the function > > > theta345toSigma <- function(theta) { > > > cholSigma <- cbind(c(theta[3], 0), c(theta[4], theta[5])) > > > crossprod(cholSigma) # t(cholSigma) %*% cholSigma) > > > } > > > This creates a positive definite matrix for any theta (and > > > any positive definite matrix has such a representation, generally > > > more than one). It is like using the square root of a quantity > > > in the optimizer when you know the quantity must be non-negative. > > > > > > Then change your > > > sigma <- c(theta[3],theta[5],theta[5],theta[4]) > > > dim(sigma) <- c(2, 2) > > > to > > > sigma <- theta345toSigma(theta) > > > > > > If one of your variances is near 0 the optimizer may run into > > > trouble at saddlepoints. Others may be able to help better > > > with that issue. > > > > > > Bill Dunlap > > > Spotfire, TIBCO Software > > > wdunlap tibco.com > > > > Hi Bill, > > > > I tried your suggestion and the optimizer produces a result, but it > seems substantially far > > from the anticipated result and from the result obtained when I use Inf > as a return value > > for an invalid covariance matrix. Perhaps it would work if I made other > adjustments to > > account for the 'bias' (using this term very loosely) induced by > changing an invalid > > parameter vector into something valid but incorrect? > > > > I also learned you have to be very careful with the starting value, as > the simple identity > > matrix becomes singular under the transformation. > > > > > theta<-c(0,0,1,1,0) > > > cholSigma<-cbind(c(theta[3], 0), c(theta[4], theta[5])) > > > sigma<-crossprod(cholSigma) > > > sigma > > [,1] [,2] > > [1,] 1 1 > > [2,] 1 1 > > > > In any case, the Inf trick works for now. I was asking about > 'adjustments' strictly out of > > curiosity. Code and results are below. > > > > Best Regards, > > Steven > > > > > 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) > > + if(!(is.positive.semi.definite(sigma))){return(Inf)} > > + -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]))) > > + } > > > exact2<-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]) > > + cholSigma<-cbind(c(theta[3], 0), c(theta[4], theta[5])) > > + sigma<-crossprod(cholSigma) > > + -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=theta.hat.em,objective=exact,complete=complete,deleted=deleted)$par > > [1] 1.2289422 5.4995271 0.9395155 4.8069068 1.8009213 > > > > nlminb(start=theta.hat.em,objective=exact2,complete=complete,deleted=deleted)$par > > [1] 1.2289421 5.4995265 0.9692861 1.8579876 1.1639544 > > ______________________________________________ > 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.
