Reparametrize to remove constraints for optimization (R)

Question

I am teaching myself how to run some Markov models in R, by working through the textbook "Hidden Markov Models for Time Series: An Introduction using R". I am a bit stuck on how to go about implementing something that is mentioned in the text.

So, I have the following function:

f <- function(samples,lambda,delta) -sum(log(outer(samples,lambda,dpois)%*%delta))

Which I can optimize with respect to, say, lambda using:

optim(par, fn=f, samples=x, delta=d)
where "par" is the initial guess for lambda, for some x and d.

Which works perfectly fine. However, in the part of the text corresponding to the example I am trying to reproduce, they note: "The parameters delta and lambda are constrained by sum(delta_i)=1 for i=1,...m, delta_i>0, and lambda_i>0. It is therefore necessary to reparametrize if one wishes to use an unconstrained optimizer such as nlm". One possibility is to maximize the likelihood with respect to the 2m-1 unconstrained parameters".

The unconstrained parameters are given by

eta<-log(lambda) 
tau<-log(delta/(1-sum(delta))

I don't entirely understand how to go about implementing this. How would I write a function to optimize over this transformed parameter space?

It's unclear to me the exact shape of the objects involved here. A clear [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input data would be very helpful. — MrFlick, Feb 27 '15 at 01:27
Incorporate the unconstrained parameters to your target function, i.e., express it in terms of `eta` and `tau`. As @MrFlick suggested, provide a reproducible example. — Khashaa, Feb 27 '15 at 01:40

score 1 · Answer 1 · answered Feb 27 '15 at 05:47

When using optim() without parmater transfromations like so:

simpleFun  <-  function(x)
                   (x-3)^2
out = optim(par=5,
            fn=simpleFun)

the set of parmaters estimates would be obtained via out$par which is 3 in the case, as you might expect. Alternatively, you can wrap your function f in a transformation the parameters like so:

out = optim(par=5,
            fn=function(x)
                # apply the transformation x -> x^3
                simpleFun(x^3))

and now the trick to get the correct set of optimal parmeters to your function you need to apply the same transfromation to the parameter estimates as in:

(out$par)^3
#> 2.99741

(and yes, the parameter estimate is slightly different. For this contrived example, you could set method="BFGS" for a slightly better estimate. Anyhow, this goes to show that the choice of transformation does matter in some cases, but that's for another discussion...)

To complete the answer, It sounds like you a want to use a wrapper like so

# the function to be optimized
f <- function(samples,lambda,delta) 
    -sum(log(outer(samples,lambda,dpois)%*%delta))

out <- optim(# par it now a 2m vector
             par = c(eta1 = 1,
                     eta2 = 1,
                     eta3 = 1,
                     tau1 = 1,
                     tau2 = 1,
                     tau3 = 1),

             # a wrapper that applies the constraints
             fn=function(x,samples){
                 # exp() guarantees that the values of lambda are > 0
                 lambda = exp(x[1:3]) 
                 # delta is also > 0
                 delta = exp(x[4:6])
                 # and now it sums to 1
                 delta = delta / sum(delta)
                 f(samples,lambda,delta)
             },
             samples=samples)

The above guarantees that the the parameters passed to f()have the correct constraints, and as long as you apply the same transformation to out$par, optim() will estimate an optimal set of parameters for f().

Reparametrize to remove constraints for optimization (R)

1 Answers1