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 <wdunlap@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:oreslag@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 <wdunlap@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]]

