[R] Fast optimizer

[Ravi Varadhan]

Hi,

Here is one approach to solve your likelihood maximization subject to constraints.  I have written a function called `constrOptim.nl' that can solve constrained optimization problems.  I will demonstrate this with a simulation.

# 1.  Simulate the data (x,y).

a <- 4
b <- 1
c <- 2.5
p <- 0.3
n <- 100
z <- rbinom(n, size=1, prob=p)
x <- ifelse(z, rpois(n, a), rpois(n, a+c))
y <- ifelse(z, rpois(n, b+c), rpois(n, b))

# 2.  Define the log-likelihood to be maximized

mloglik <- function(par, xx, y, ...){
# Note that I have named `xx' instead of `x' to avoid conflict with another function 
p <- par[1]
a <- par[2]
b <- par[3]
cc <- par[4]
sum(log(p*dpois(xx,a)*dpois(y,b+cc) + (1-p)*dpois(xx,a+cc)*dpois(y,b)))
}

# 3. Write the constraint function to define inequalities

hin <- function(par, ...) {
# All constraints are defined such that h[i] > 0 for all i.
h <- rep(NA, 7)
h[1] <- par[1]
h[2] <- 1 - par[1] 
h[3] <- par[2]
h[4] <- par[3]
h[5] <- par[4]
h[6] <- par[2] - par[4]
h[7] <- par[3] - par[4]
h
}

# Now, we are almost ready to maximize the log-likelihood.  We need a feasible starting value.  
# We construct a projection function that can convert any arbitrary real vector to a feasible starting vector.

# 5.  Finding a feasible starting vector of parameter values

project <- function(par, lower, upper) {
# a function to map a random vector to a feasible vector
slack <- 1.e-14
if (par[1] <= 0) par[1] <- slack
if (par[1] >= 1) par[1] <- 1 - slack
if (par[2] <= 0) par[2] <- slack
if (par[3] <= 0) par[3] <- slack
par[4] <- min(par[2], par[3], par[4]) - slack
par
} 

p0c <- runif(4, c(0,1,1,1), c(1, 5, 5, 2))  # a random candidate starting value
p0 <- project(p0c)  # feasible starting value

# 6.  Optimization

source("constrOptim_nl.R")  # we source the `constrOptim.nl' function
ans <- constrOptim.nl(par=p0, fn=mloglik, hin=hin, xx=x, y=y, control=list(fnscale=-1)) # Note: we are maximizing
ans


This concludes our brief tutorial.   

Hope this helps,
Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Ravi Varadhan <rvaradhan at jhmi.edu>
Date: Thursday, October 29, 2009 10:23 pm
Subject: Re: [R] Fast optimizer
To: R_help Help <rhelpacc at gmail.com>
Cc: r-help at r-project.org, r-sig-finance at stat.math.ethz.ch


> This optimization should not take you 1-2 mins.  My guess is that your 
> coding of the likelihood is inefficient, perhaps containing for loops. 
>  Would you mind sharing your code with us?
> 
> As far as incorporating inequality constraints, there are at least 4 approaches:
> 
> 1.  You could use `constrOptim' to incorporate inequality constraints. 
>  
> 
> 2.  I have written a function `constrOptim.nl' that is better than 
> `constrOptim', and it can be used here.
> 
> 3.  The projection method in spg() function in the "BB" package can be 
> used.
> 
> 4.  You can create a penalized objective function, where the 
> inequalities c < a and c < b can be penalized.  Then you can use 
> optim's L-BFGS-B or spg() or nlminb().
> 
> 
> Ravi.
> ____________________________________________________________________
> 
> Ravi Varadhan, Ph.D.
> Assistant Professor,
> Division of Geriatric Medicine and Gerontology
> School of Medicine
> Johns Hopkins University
> 
> Ph. (410) 502-2619
> email: rvaradhan at jhmi.edu
> 
> 
> ----- Original Message -----
> From: R_help Help <rhelpacc at gmail.com>
> Date: Thursday, October 29, 2009 9:21 pm
> Subject: Re: [R] Fast optimizer
> To: Ravi Varadhan <rvaradhan at jhmi.edu>
> Cc: r-help at r-project.org, r-sig-finance at stat.math.ethz.ch
> 
> 
> > Ok. I have the following likelihood function.
> > 
> > L <- p*dpois(x,a)*dpois(y,b+c)+(1-p)*dpois(x,a+c)*dpois(y,b)
> > 
> > where I have 100 points of (x,y) and parameters c(a,b,c,p) to
> > estimate. Constraints are:
> > 
> > 0 < p < 1
> > a,b,c > 0
> > c < a
> > c < b
> > 
> > I construct a loglikelihood function out of this. First ignoring the
> > last two constraints, it takes optim with box constraint about 1-2 min
> > to estimate this. I have to estimate the MLE on about 200 rolling
> > windows. This will take very long. Is there any faster implementation?
> > 
> > Secondly, I cannot incorporate the last two contraints using optim function.
> > 
> > Thank you,
> > 
> > rc
> > 
> > 
> > 
> > On Thu, Oct 29, 2009 at 9:02 PM, Ravi Varadhan <rvaradhan at jhmi.edu> 
> wrote:
> > >
> > > You have hardly given us any information for us to be able to help 
> 
> > you.  Give us more information on your problem, and, if possible, a 
> 
> > minimal, self-contained example of what you are trying to do.
> > >
> > > Ravi.
> > > ____________________________________________________________________
> > >
> > > Ravi Varadhan, Ph.D.
> > > Assistant Professor,
> > > Division of Geriatric Medicine and Gerontology
> > > School of Medicine
> > > Johns Hopkins University
> > >
> > > Ph. (410) 502-2619
> > > email: rvaradhan at jhmi.edu
> > >
> > >
> > > ----- Original Message -----
> > > From: R_help Help <rhelpacc at gmail.com>
> > > Date: Thursday, October 29, 2009 8:24 pm
> > > Subject: [R] Fast optimizer
> > > To: r-help at r-project.org, r-sig-finance at stat.math.ethz.ch
> > >
> > >
> > >> Hi,
> > >>
> > >> I'm using optim with box constraints to MLE on about 100 data points.
> > >> It goes quite slow even on 4GB machine. I'm wondering if R has any
> > >> faster implementation? Also, if I'd like to impose
> > >> equality/nonequality constraints on parameters, which package I should
> > >> use? Any help would be appreciated. Thank you.
> > >>
> > >> rc
> > >>
> > >> ______________________________________________
> > >> R-help at r-project.org mailing list
> > >>
> > >> PLEASE do read the posting guide
> > >> and provide commented, minimal, self-contained, reproducible code.
> > >
> 
> ______________________________________________
> R-help at r-project.org mailing list
> 
> PLEASE do read the posting guide 
> and provide commented, minimal, self-contained, reproducible code.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: constrOptim_nl.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091101/f5377cbf/attachment-0002.txt>