I am using nlsLM {minpack.lm}
to find the values of parameters a
and b
of function myfun
which give the best fit for the data set, mydata
.
mydata=data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45))
myfun=function(a,b,r,t){
prd=a*b*(1-exp(-b*r*t))
return(prd)
}
and using nlsLM
myfit=nlsLM(y~myfun(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
lower = c(1000,0), upper = c(3000,1))
It works. But now I would like to introduce a constraint which is a*b<1000
. I had a look at the option available in nlsLM
to set constraint via nls.lm.control
. But it's not much of help. can somebody help me here or suggest a different method to to this?
It seems that
nlsLM
doesn't accept constraints, one possible method is to useconstrOptim
function instead, but the drawback is that you have to turn the problem into the form thatconstrOptim
accepts. By the way, I'm currently working on an R package that gives a nice interface foroptim
andconstrOptim
and make necessary transformation automatically, but it is not ready to use yet. Here is how to useconstrOptim
for this problem.The first step is that
constrOptim
only accepts linear constraints,a * b < 1000
=>log(a) + log(b) < log(1000)
, so you need to reparametrize the problem usinglog(a)
andlog(b)
, the constraints becomes:-log(a) - log(b) > -log(1000)
,log(a) >= log(1000)
,-log(a) >= -log(3000)
,-log(b) >= 0
,ui
andci
are matrix form for these constraints. Andmyfun1
is the reparametrize version for your function, and then you can useconstrOptim
to obtain your solution forlog(a)
andlog(b)
.To my knowledge, in nlsLM of minpack.lm package lower and upper parameter bounds are the only available constrains.
One possible solution to your problem is based on the use of the package
nloptr
, an R interface to the free open-source library NLopt for nonlinear optimization.In the code below I use
cobyla
, an algorithm for derivative-free optimization with nonlinear inequality and equality constraints:The values of the optimal parameters are:
and the sum of squared errors is:
Hope it can help you.