[R] nlminb() - how do I constrain the parameter vector properly?
William Dunlap
wdunlap at tibco.com
Mon Oct 21 16:52:11 CEST 2013
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
> -----Original Message-----
> From: Steven LeBlanc [mailto:oreslag at gmail.com]
> Sent: Sunday, October 20, 2013 9:35 PM
> To: William Dunlap
> Subject: Re: [R] nlminb() - how do I constrain the parameter vector properly?
>
>
> On Oct 20, 2013, at 6:41 PM, William Dunlap <wdunlap at 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
>
> Hi Bill,
>
> I've clipped out the superfluous code to leave the objective function and call to nlminb()
> below.
>
> I'm using the first two parameters to construct the vector of means 'u' for a bivariate
> normal, and the final three parameters to construct the corresponding covariance matrix
> 'sigma'. Both are required by dmvnorm(). In summary, nlminb() required a vector of
> parameters, so I supplied the number of parameters I needed nlminb() to optimize and
> simply built the required formats within the function.
>
> It didn't occur to me until I saw the error that nlminb() would have no way of knowing the
> proper boundaries of the parameter space unless there is some way to communicate the
> constraints. nlminb() implements simple box constraints, but I don't see a way to
> communicate "parameters 3, 4, and 5 must satisfy 3*4 - 5^2 > 0.
>
> Regarding your suggestion, I don't think I understand. Might you elaborate?
>
> Thanks & Best Regards,
> Steven
>
> > exact<-function(theta,complete,deleted){
> >>
> >> 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))
More information about the R-help
mailing list