[R] nlminb() - how do I constrain the parameter vector properly?
William Dunlap
wdunlap at tibco.com
Mon Oct 21 18:28:38 CEST 2013
[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 at 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 at 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
More information about the R-help
mailing list