6

I have been using the Excel solver to handle the following problem

solve for a b and c in the equation:

y = a*b*c*x/((1 - c*x)(1 - c*x + b*c*x))

subject to the constraints

0 < a < 100
0 < b < 100
0 < c < 100

f(x[1]) < 10
f(x[2]) > 20
f(x[3]) < 40

where I have about 10 (x,y) value pairs. I minimize the sum of abs(y - f(x)). And I can constrain both the coefficients and the range of values for the result of my function at each x.

I tried nls (without trying to impose the constraints) and while Excel provided estimates for almost any starting values I cared to provide, nls almost never returned an answer.

I switched to using optim, but I'm having trouble applying the constraints.

This is where I have gotten so far-

best = function(p,x,y){sum(abs(y - p[1]*p[2]*p[3]*x/((1 - p[3]*x)*(1 - p[3]*x + p[2]*p[3]*x))))}
p = c(1,1,1)
x = c(.1,.5,.9)
y = c(5,26,35)
optim(p,best,x=x,y=y)

I did this to add the first set of constraints-

optim(p,best,x=x,y=y,method="L-BFGS-B",lower=c(0,0,0),upper=c(100,100,100))

I get the error ""ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

and end up with a higher value of the error ($value). So it seems like I am doing something wrong. I couldn't figure out how to apply my other set of constraints at all.

Could someone provide me a basic idea how to solve this problem that a non-statistician can understand? I looked at a lot of posts and looked in a few R books. The R books stopped at the simplest use of optim.

joran
  • 169,992
  • 32
  • 429
  • 468
Jon Swanson
  • 324
  • 3
  • 6
  • 16
  • general insight into questions about optimization with constraints [here](http://stackoverflow.com/questions/9817001/optimization-with-constraints) – Chase Jul 08 '12 at 23:31

1 Answers1

14

The absolute value introduces a singularity: you may want to use a square instead, especially for gradient-based methods (such as L-BFGS).

The denominator of your function can be zero.

The fact that the parameters appear in products and that you allow them to be (arbitrarily close to) zero can also cause problems.

You can try with other optimizers (complete list on the optimization task view), until you find one for which the optimization converges.

x0 <- c(.1,.5,.9)
y0 <- c(5,26,35)
p <- c(1,1,1)
lower <- 0*p
upper <- 100 + lower
f <- function(p,x=x0,y=y0) sum( 
  (
    y - p[1]*p[2]*p[3]*x / ( (1 - p[3]*x)*(1 - p[3]*x + p[2]*p[3]*x) ) 
  )^2
)

library(dfoptim)
nmkb(p, f, lower=lower, upper=upper) # Converges

library(Rvmmin)
Rvmmin(p, f, lower=lower, upper=upper) # Does not converge

library(DEoptim)
DEoptim(f, lower, upper) # Does not converge

library(NMOF)
PSopt(f, list(min=lower, max=upper))[c("xbest", "OFvalue")] # Does not really converge
DEopt(f, list(min=lower, max=upper))[c("xbest", "OFvalue")] # Does not really converge

library(minqa)
bobyqa(p, f, lower, upper) # Does not really converge

As a last resort, you can always use a grid search.

library(NMOF)
r <- gridSearch( f, 
  lapply(seq_along(p), function(i) seq(lower[i],upper[i],length=200)) 
)
Vincent Zoonekynd
  • 31,893
  • 5
  • 69
  • 78
  • When I was using Excel, I tried minimizing both the sum of the absolute diffrences and the sum of the squares of the absolute differences. The data I am getting sometimes has a data point with high uncertainty and the square was trying too hard to fit it. What I was also able to do in Excel was use the individual estimated values as constraints. For example, y0 [3] is supposed to be 35. I could apply a constraint y0[3] < 40. I see how to set bounds on the parameters. I assume I need to modify the function I minimize to handle this other case, but don't know how I'd do it. – Jon Swanson Jul 10 '12 at 15:54
  • You can add the constraint breaches as penalties to the objective function ("soft constraints"), e.g., `pmax(y0[3] - 40, 0)^2`. You can also try other cost functions: the absolute value, a smoothed absolute value, `cost <- function(x) ifelse(abs(x)<.5,x^2,abs(x)-.25)` or functions that grow even more slowly, e.g., `cost <- function(x) ifelse(abs(x)<.5,x^2,sqrt(abs(x))-sqrt(.5)+.25)`. If you know which observations are noisiser, you can use weights in the objective function, to reduce their influence. – Vincent Zoonekynd Jul 11 '12 at 00:05