[R] R_using non linear regression with constraints
David Winsemius
dwinsemius at comcast.net
Mon Jun 19 00:09:32 CEST 2017
Dear Dr Nash;
> On Jun 18, 2017, at 11:52 AM, J C Nash <profjcnash at gmail.com> wrote:
>
> I ran the following script. I satisfied the constraint by
> making a*b a single parameter, which isn't always possible.
> I also ran nlxb() from nlsr package, and this gives singular
> values of the Jacobian. In the unconstrained case, the svs are
> pretty awful, and I wouldn't trust the results as a model, though
> the minimum is probably OK. The constrained result has a much
> larger sum of squares.
I have version 2017.2.19 of nlsr on my Mac. Not sure what version is throwing errors.
>
> Notes:
> 1) nlsr has been flagged with a check error by CRAN (though it
> is in the vignette, and also mentions pandoc a couple of times).
> I'm working to purge the "bug", and found one on our part, but
> not necessarily all the issues.
> 2) I used nlxb that requires an expression for the model. nlsLM
> can use a function because it is using derivative approximations,
> while nlxb actually gets a symbolic or automatic derivative if
> it can, else squawks.
I think of such problems as potentially searching for the maximum along a constraint boundary as a problem at an intersection of a family 2 surfaces in the "a cross b" space, in this case the surface ( a*b = 1000 ) and the family of curves that are implicitly defined by a*b*(1-exp(-b*r*t)) with r fixed and t <- mydata$x.
The intersection of these two geometric structures ( a constraint surface and a family parametrized by a and b) is in my construction the problem domain.
I can see where an algorithm might erroneously return a "satisfied" results when it was following an iterative path to the edge of the constraint boundary but was prematurely "seeing the task as complete". Can you comment on the logical difficulties and solutions to surmounting such error after such an event occurs?
Sincerely;
David Winsemius
>
> JN
>
> # Here's the script #
> #
> # Manoranjan Muthusamy <ranjanmano167 at gmail.com>
> #
>
> library(minpack.lm)
> 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))
> summary(myfit)
> library(nlsr)
> r <- 2
> myfitj=nlxb(y~a*b*(1-exp(-b*r*x)),data=mydata,start=list(a=2000,b=0.05), trace=TRUE)
> summary(myfitj)
> print(myfitj)
>
> myfitj2<-nlxb(y~ab*(1-exp(-b*r*x)),data=mydata,start=list(ab=2000*0.05,b=0.05), trace=TRUE)
> summary(myfitj2)
> print(myfitj2)
>
> myfitj2b<-nlxb(y~ab*(1-exp(-b*r*x)),data=mydata,start=list(ab=2000*0.05,b=0.05),
> trace=TRUE, upper=c(1000, Inf))
> summary(myfitj2b)
> print(myfitj2b)
> # End of script #
>
> On 2017-06-18 01:29 PM, Bert Gunter wrote:
>> https://cran.r-project.org/web/views/Optimization.html
>> (Cran's optimization task view -- as always, you should search before posting)
>> In general, nonlinear optimization with nonlinear constraints is hard,
>> and the strategy used here (multiplying by a*b < 1000) may not work --
>> it introduces a discontinuity into the objective function, so
>> gradient based methods may in particular be problematic. As usual, if
>> either John Nash or Ravi Varadhan comment, heed what they suggest. I'm
>> pretty ignorant.
>> Cheers,
>> Bert
>> Bert Gunter
>> "The trouble with having an open mind is that people keep coming along
>> and sticking things into it."
>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>> On Sun, Jun 18, 2017 at 9:43 AM, David Winsemius <dwinsemius at comcast.net> wrote:
>>>
>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <ranjanmano167 at gmail.com> wrote:
>>>>
>>>> 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.
>>>
>>> At the moment your coefficients do satisfy that constraint so that dataset is not suitable for testing. A slight modification of the objective function to include the logical constraint as an additional factor does not "break" that particular solution.:
>>>
>>> myfun2=function(a,b,r,t){
>>> prd=a*b*(1-exp(-b*r*t))*(a*b<1000)
>>> return(prd)}
>>>
>>>
>>> myfit=nlsLM(y~myfun2(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
>>> lower = c(1000,0), upper = c(3000,1))
>>>
>>> #------------------
>>> myfit
>>> Nonlinear regression model
>>> model: y ~ myfun2(a, b, r = 2, t = x)
>>> data: mydata
>>> a b
>>> 3.000e+03 2.288e-02
>>> residual sum-of-squares: 38.02
>>>
>>> Number of iterations to convergence: 8
>>> Achieved convergence tolerance: 1.49e-08
>>> #--
>>>
>>> prod(coef(myfit))
>>> #[1] 68.64909 Same as original result.
>>>
>>> How nlsLM will handle more difficult problems is not something I have experience with, but obviously one would need to keep the starting values within the feasible domain. However, if your goal was to also remove the upper and lower constraints on a and b, This problem would not be suitably solved by the a*b product without relaxation of the default maxiter:
>>>
>>>> myfit=nlsLM(y~myfun2(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
>>> + lower = c(0,0), upper = c(9000,1))
>>>> prod(coef(myfit))
>>> [1] 110.4382
>>>> coef(myfit)
>>> a b
>>> 9.000000e+03 1.227091e-02
>>>
>>>> myfit=nlsLM(y~myfun2(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
>>> + lower = c(0,0), upper = c(10^6,1))
>>> Warning message:
>>> In nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower, :
>>> lmdif: info = -1. Number of iterations has reached `maxiter' == 50.
>>>
>>> #---------
>>> myfit=nlsLM(y~myfun2(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
>>> lower = c(0,0), upper = c(10^6,1), control=list(maxiter=100))
>>> prod(coef(myfit))
>>>
>>> coef(myfit)
>>> #============
>>>
>>>
>>>> prod(coef(myfit))
>>> [1] 780.6732 Significantly different than the solution at default maxiter of 50.
>>>>
>>>> coef(myfit)
>>> a b
>>> 5.319664e+05 1.467524e-03
>>>>
>>>>
>>>
>>>
>>> --
>>> David.
>>>
>>>
>>>> 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?
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>> David Winsemius
>>> Alameda, CA, USA
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
David Winsemius
Alameda, CA, USA
More information about the R-help
mailing list