Fast optimizer

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>