CODEWITHSUNDEEP
rnls
Adela
Adela

Reputation: 1797

Constraints for nls coefficients

I'm trying to fit data with nls() function where the nature of data gives me bounds for one coefficient and for sum of two coefficients. Let me introduce short example to see where is the problem. I want parameter b1 to be between 0 and 1 and I want sum of parameters b1 and b2 to be between 0 and 1 as well.

set.seed(123)

# example where everything is OK
x <- 1:200
g <- rbinom(200, 1, 0.5)
y <- 3 + (0.7 + 0.2 * g) * x 
yeps <- y + rnorm(length(y), sd = 0.1)  

# both parameter b1 and sum of parameters b1 and b2 are between 0 and 1
nls(yeps ~ a + (b1 + b2 * g) * x, start = list(a = 0.12345, b1 = 0.54321, b2 = 0.4213))

# using more extreme values
x <- 1:200
g <- rbinom(200, 1, 0.5)
y <- 3 + (0.9 - 0.99 * g) * x 
yeps <- y + rnorm(length(y), sd = 15) 

# b1 is OK, but b1 + b2 < 0
nls(yeps ~ a + (b1 + b2 * g) * x, 
    start = list(a = 0.12345, b1 = 0.54321, b2 = 0.4213))

# trying constraints, not good, sum is still out of range
nls(yeps ~ a + (b1 + b2 * g) * x, 
    start = list(a = 0.12345, b1 = 0.54321, b2 = 0.4213),
    lower = list(a = -Inf, b1 = 0, b2 = -1),
    upper = list(a = Inf, b1 = 1, b2 = 1),
    algorithm = "port")

What I'm looking for is something like that (does not work):

nls(yeps ~ a + (b1 + b2 * g) * x, 
    start = list(a = 0.12345, b1 = 0.54321, b2 = 0.4213),
    lower = list(a = -Inf, b1 = 0, b2 = -b1),
    upper = list(a = Inf, b1 = 1, b2 = 1 - b1),
    algorithm = "port")

Is it possible to set constraints with other parameters in nls() function? Thanks for any suggestions!

Upvotes: 1

Views: 1433

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269431

Let B2 = b1+b2 so b2 = B2-b1 and substituting B2-b1 for b2 we get a problem in terms of a, b1 and B2 of which the latter two are between 0 and 1 so:

fm <- nls(yeps ~ a + (b1 + (B2-b1) * g) * x, lower = c(-Inf, 0, 0), upper = c(Inf, 1, 1),
    start = list(a = 0.1, b1 = 0.5, B2 = 0.1), alg = "port")

giving the following (hence b2 = B2 - b1 = 0 - 0.9788 = -0.9788)

> fm
Nonlinear regression model
  model: yeps ~ a + (b1 + (B2 - b1) * g) * x
   data: parent.frame()
      a      b1      B2 
-5.3699  0.9788  0.0000 
 residual sum-of-squares: 42143

Algorithm "port", convergence message: both X-convergence and relative convergence (5)

and plotting:

plot(yeps ~ x)
points(fitted(fm) ~ x, pch = 20, col = "red")

screenshot

Upvotes: 2

Related Questions