0

I want to perform a non-linear least square estimation using the function nlsLM from the package minpack.lm.

I want to impose an upper and a lower bound to the estimation to force the algorithm to find the solution within a specific support.

Something along these lines:

library(minpack.lm)
mydata <- data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45))
myfit <- nlsLM(formula(y ~ a*x), data = mydata, start=list(a = 2.5), lower = c(a = 0), upper = c(a = 5))
summary(myfit)

My question is:

is it possible to apply a penalty function to nlsLM, so to avoid that the algorithm returns a corner solution? e.g. in my example a is not equal to 0 or 5.

Caveats: I am exclusively looking at a solution with the package minpack.lm.

shazz
  • 387
  • 1
  • 13

1 Answers1

2

To implement this barrier function which approaches infinity as a approaches 0 from the right or 5 from the left:

log(a)^2 + log(5-a)^2

as a penalty we form the objective function

(y[1] - a*x[1])^2 + ... + (y[n] - a*x[n])^2 + (0 - log(a))^2 + (0 - log(5-a))^2

like this

n <- nrow(mydata)
mydata0 <- rbind(mydata, 0 * mydata[1:2, ])
nlsLM(y ~ c(a * x[1:n], log(a), log(5 - a)), mydata0, start = list(a = 1))
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • Thank you for your help. It is not all clear to me. First, why I need to create mydata0? Second, it looks to me that the barrier is at log(a), and log(5-a), not log(a)^2 and log(5-a)^2, am I correct? Third, can you explain how e.g. log(5-a) creates a penalty? Suppose for instance "a" goes to 4.8, this means that we are applying a penalty of log(5-4.8) = -1.61? And this penalty to what is applied, the RSS? Fourth, can you clarify how bounds are used? To create e.g. a barrier function between -1 and 1, should I replace log(a) and log(5-a) with e.g. "(-1-a-0.05)" and "(1-a-0.05)" – shazz Dec 11 '20 at 08:35
  • Have added explanation. Be sure to read the link. – G. Grothendieck Dec 11 '20 at 10:23
  • The added explanation is extremely helpful and clear, and the link is a great reference. Thank you very very much. Still, it is unclear to me what is the logic behind mydata0. Why I need two new rows in my dataset, and why these have to contain values equal to zero? – shazz Dec 11 '20 at 10:52
  • So that we can form the loast two terms in the objective function. – G. Grothendieck Dec 11 '20 at 10:53
  • Ok now it is all clear. I have accepted your answer. Indeed, thank you very much. – shazz Dec 11 '20 at 11:03