R_: nls with constraints on a sum of parameters

Question

I am trying to use the minpack.lm package to fit a non linear model to data with a constrain on the sum of two parameters. Here A_b+A_s<1, A_b>0, A_s>0 and k>0 are the constraints I want to optimize for, thus the log(A_b + A_s) and log(1 - (A_b + A_s)) functions.

I tried to reproduce the example given here.

However, the nlsLM throws an error, and I suspect it is because it is not able to put constraint on a sum of parameters, as opposed to a constraint to the value of a parameter only.

library(minpack.lm)
data <- data.frame(x = c(0.63, 0.67, 1.08, 1.37, 1.17, 1.34, 1.33, 1.48, 1.40, 1.58, 1.02, 1.48, 0.90, 0.55, 0.90), y = c(0.58, 0.61, 0.50, 0.51, 0.50, 0.47, 0.49, 0.51, 0.47, 0.48, 0.53, 0.43, 0.57, 0.68, 0.49))

n <- nrow(data)
mydata0 <- rbind(data, 0 * data[1:2,])

nlsLM(y ~ c(A_b + (A_s * k / x[1:n]) * log((x[1:n] + k) / k) ,
            log(A_b + A_s),
            log(1 - (A_b + A_s))),
      mydata0,
      start = list(A_b = 0.5, A_s = 0.5, k = 1))

Answer 1

Indeed, you can't include inequality constraints with minpack.lm. You can with the cobyla function of the nloptr package:

library(nloptr)

x <- c(0.63, 0.67, 1.08, 1.37, 1.17, 1.34, 1.33, 1.48, 1.40, 1.58, 1.02, 1.48, 0.90, 0.55, 0.90) 
y <- c(0.58, 0.61, 0.50, 0.51, 0.50, 0.47, 0.49, 0.51, 0.47, 0.48, 0.53, 0.43, 0.57, 0.68, 0.49)

# quantity to be minimized: sum of squared residuals
fn <- function(pars) {
  A_b <- pars[1]
  A_s <- pars[2]
  k <- pars[3]
  target <- A_b + (A_s * k / x) * log((x + k) / k)
  sum((y - target)^2)
}

# initial values
pars0 <- c(0.2, 0.8, 0.3)

# inequality constraints: each component >= 0
# note: we could use only h[1] and set lower = c(0, 0, 0) in the call to cobyla
ineqs <- function(pars) {
  # A_b+A_s<1, A_b>0, A_s>0 and k>0
  A_b <- pars[1]
  A_s <- pars[2]
  k <- pars[3]
  h <- numeric(4)
  h[1] <- 1 - A_b - A_s
  h[2] <- A_b
  h[3] <- A_s
  h[4] <- k
  h
}

# optimization
cobyla(
  x0 = pars0, fn = fn, hin = ineqs, 
  control = list(xtol_rel = 1e-8, maxeval = 10000),
  nl.info = TRUE
)