R_using non linear regression with constraints

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

R_using non linear regression with constraints

ranjanmano167
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?

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

David Winsemius

> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

Bert Gunter-2
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 <[hidden email]> wrote:

>
>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>
>> ______________________________________________
>> [hidden email] 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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

J C Nash
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.

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.

JN

#  Here's the script #
#
# Manoranjan Muthusamy <[hidden email]>
#

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 <[hidden email]> wrote:
>>
>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>>
>>> ______________________________________________
>>> [hidden email] 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
>>
>> ______________________________________________
>> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

David Winsemius
Dear Dr Nash;

> On Jun 18, 2017, at 11:52 AM, J C Nash <[hidden email]> 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 <[hidden email]>
> #
>
> 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 <[hidden email]> wrote:
>>>
>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>> ______________________________________________
>> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

Jeff Newmiller
In reply to this post by J C Nash
I am not as expert as John, but I thought it worth pointing out that the
variable substitution technique gives up one set of constraints for
another (b=0 in this case). I also find that plots help me see what is
going on, so here is my reproducible example (note inclusion of library
calls for completeness). Note that NONE of the optimizers mentioned so far
appear to be finding the true best fit. The fact that myfun() yields 0
always if t=0 and that condition is within the data given seems likely to
be part of the problem. I don't know how to resolve this... perhaps John
will look at it again.

David: Your thinking makes fine sense if you are using a Monte Carlo or
brute force solution, but the fact that it creates a discontinuity in
the objective function will confuse any optimizer that uses analytic or
numerically estimated slopes.

##----------begin
library(minpack.lm)
library(ggplot2)

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 ) {
   a * b * ( 1 - exp( -b * r * t ) )
}

objdta <- expand.grid( a = seq( 1000, 3000, by=20 )
                      , b = seq( -0.01, 1, 0.01 )
                      , rowidx = seq.int( nrow( mydata ) )
                      )
objdta[ , c( "y", "t" ) ] <- mydata[ objdta$rowidx
                                    , c( "y", "x" ) ]
objdta$tf <- factor( objdta$t )
objdta$myfun <- with( objdta
                     , myfun( a = a, b = b, r = 2, t = t )
                     )
objdtass <- aggregate( ( objdta$myfun - objdta$y )^2
                      , objdta[ , c( "a", "b" ) ]
                      , FUN = function( x )
                               sum( x, na.rm=TRUE )
                      )
objdtassmin <- objdtass[ which.min( objdtass$x ), ]

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 )
               )
a <- as.vector( coef( myfit )[ "a" ] )
b <- as.vector( coef( myfit )[ "b" ] )

brks <- c( 500, 1e7, 2e7, 3e7, 4e7 )
ggplot( objdtass, aes( x=a, y=b, z = x, fill=x ) ) +
     geom_tile() +
     geom_contour( breaks= brks ) +
     geom_point( x=a, y=b, colour="red" ) +
     geom_point( x=objdtassmin$a
               , y=objdtassmin$b
               , colour="green" ) +
     scale_fill_continuous( name="SumSq", breaks = brks )
# Green point is brute-force solution
# Red point is optimizer solution for myfun

##############

myfun2 <- function( a, log1ab, r, t ) {
   ab <- 1000 - exp( log1ab )
   ab * ( 1 - exp( -ab/a * r * t ) )
}

objdta$log1ab <- with( objdta, log( 1000 - a * b ) )
objdta$myfun2 <- with( objdta
                      , myfun2( a = a
                              , log1ab = log1ab
                              , r = 2
                              , t = t
                              )
                      )
objdtass2 <- aggregate( ( objdta$myfun2 - objdta$y )^2
                       , objdta[ , c( "a", "b" ) ]
                       , FUN = function( x )
                                if ( all( is.na( x ) ) ) NA
                                else sum( x, na.rm=TRUE )
                       )
objdtass2min <- objdtass2[ which.min( objdtass2$x ), ]

myfit2 <- nlsLM( y ~ myfun2( a, log1ab, r = 2, t = x )
                , data = mydata
                , start = list( a = 2000
                              , log1ab = 4.60517
                              )
                , lower = c( 1000, 0 )
                , upper = c( 3000, 8.0063 )
                )
a2 <- as.vector( coef( myfit2 )[ "a" ] )
b2 <- ( 1000
       - exp( as.vector( coef( myfit2 )[ "log1ab" ] ) )
       ) / a

brks <- c( 500, 1e6, 2e6, 3e6, 4e6 )
ggplot( objdtass2, aes( x=a, y=b, z = x, fill=x ) ) +
     geom_tile() +
     geom_contour( breaks = brks ) +
     geom_point( x=a2, y=b2, colour="red" ) +
     geom_point( x=objdtass2min$a
               , y=objdtass2min$b
               , colour="green" ) +
     scale_fill_continuous( name="SumSq", breaks = brks )
# Green point is brute-force solution
# Red point is optimizer solution for myfun2

##----------end

On Sun, 18 Jun 2017, J C Nash 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.
>
> 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.
>
> JN
>
> #  Here's the script #
> #
> # Manoranjan Muthusamy <[hidden email]>
> #
>
> 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 <[hidden email]>
>> wrote:
>>>
>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy
>>>> <[hidden email]> 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]]
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>
> ______________________________________________
> [hidden email] 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.
>

---------------------------------------------------------------------------
Jeff Newmiller                        The     .....       .....  Go Live...
DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...
                                       Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
/Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

J C Nash
I've seen a number of problems like this over the years. The fact that the singular values of the Jacobian have
a ration larger than the usual convergence tolerances can mean the codes stop well before the best fit. That is
the "numerical analyst" view. David and Jeff have given geometric and statistical arguments. All views are useful,
but it takes some time to sort them all out and make sense of the problem as a whole. Right now I'm getting ready
to go to UseR!, so probably won't have time to spend sorting things out, though I'll have a go if I get time.

David: Did you get a crash with the Mac version of nlsr? (You have the latest version I uploaded to CRAN.)
I don't have Mac, just Linux and a very ancient Win XP. I suspect the latter is too old to take seriously. I use
Win-Builder for package checks. I don't think there's anything seriously wrong with nlsr functions, but there is
more to be done to the vignettes and better manual pages are always desirable. It is the vignettes that are
popping up warnings on CRAN checks, and I'll see if I can fix those.

Cheers, JN


On 2017-06-18 06:23 PM, Jeff Newmiller wrote:

> I am not as expert as John, but I thought it worth pointing out that the variable substitution technique gives up one
> set of constraints for another (b=0 in this case). I also find that plots help me see what is going on, so here is my
> reproducible example (note inclusion of library calls for completeness). Note that NONE of the optimizers mentioned so far
> appear to be finding the true best fit. The fact that myfun() yields 0 always if t=0 and that condition is within the
> data given seems likely to be part of the problem. I don't know how to resolve this... perhaps John will look at it again.
>
> David: Your thinking makes fine sense if you are using a Monte Carlo or brute force solution, but the fact that it
> creates a discontinuity in the objective function will confuse any optimizer that uses analytic or numerically estimated
> slopes.
>
> ##----------begin
> library(minpack.lm)
> library(ggplot2)
>
> 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 ) {
>    a * b * ( 1 - exp( -b * r * t ) )
> }
>
> objdta <- expand.grid( a = seq( 1000, 3000, by=20 )
>                       , b = seq( -0.01, 1, 0.01 )
>                       , rowidx = seq.int( nrow( mydata ) )
>                       )
> objdta[ , c( "y", "t" ) ] <- mydata[ objdta$rowidx
>                                     , c( "y", "x" ) ]
> objdta$tf <- factor( objdta$t )
> objdta$myfun <- with( objdta
>                      , myfun( a = a, b = b, r = 2, t = t )
>                      )
> objdtass <- aggregate( ( objdta$myfun - objdta$y )^2
>                       , objdta[ , c( "a", "b" ) ]
>                       , FUN = function( x )
>                                sum( x, na.rm=TRUE )
>                       )
> objdtassmin <- objdtass[ which.min( objdtass$x ), ]
>
> 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 )
>                )
> a <- as.vector( coef( myfit )[ "a" ] )
> b <- as.vector( coef( myfit )[ "b" ] )
>
> brks <- c( 500, 1e7, 2e7, 3e7, 4e7 )
> ggplot( objdtass, aes( x=a, y=b, z = x, fill=x ) ) +
>      geom_tile() +
>      geom_contour( breaks= brks ) +
>      geom_point( x=a, y=b, colour="red" ) +
>      geom_point( x=objdtassmin$a
>                , y=objdtassmin$b
>                , colour="green" ) +
>      scale_fill_continuous( name="SumSq", breaks = brks )
> # Green point is brute-force solution
> # Red point is optimizer solution for myfun
>
> ##############
>
> myfun2 <- function( a, log1ab, r, t ) {
>    ab <- 1000 - exp( log1ab )
>    ab * ( 1 - exp( -ab/a * r * t ) )
> }
>
> objdta$log1ab <- with( objdta, log( 1000 - a * b ) )
> objdta$myfun2 <- with( objdta
>                       , myfun2( a = a
>                               , log1ab = log1ab
>                               , r = 2
>                               , t = t
>                               )
>                       )
> objdtass2 <- aggregate( ( objdta$myfun2 - objdta$y )^2
>                        , objdta[ , c( "a", "b" ) ]
>                        , FUN = function( x )
>                                 if ( all( is.na( x ) ) ) NA
>                                 else sum( x, na.rm=TRUE )
>                        )
> objdtass2min <- objdtass2[ which.min( objdtass2$x ), ]
>
> myfit2 <- nlsLM( y ~ myfun2( a, log1ab, r = 2, t = x )
>                 , data = mydata
>                 , start = list( a = 2000
>                               , log1ab = 4.60517
>                               )
>                 , lower = c( 1000, 0 )
>                 , upper = c( 3000, 8.0063 )
>                 )
> a2 <- as.vector( coef( myfit2 )[ "a" ] )
> b2 <- ( 1000
>        - exp( as.vector( coef( myfit2 )[ "log1ab" ] ) )
>        ) / a
>
> brks <- c( 500, 1e6, 2e6, 3e6, 4e6 )
> ggplot( objdtass2, aes( x=a, y=b, z = x, fill=x ) ) +
>      geom_tile() +
>      geom_contour( breaks = brks ) +
>      geom_point( x=a2, y=b2, colour="red" ) +
>      geom_point( x=objdtass2min$a
>                , y=objdtass2min$b
>                , colour="green" ) +
>      scale_fill_continuous( name="SumSq", breaks = brks )
> # Green point is brute-force solution
> # Red point is optimizer solution for myfun2
>
> ##----------end
>
> On Sun, 18 Jun 2017, J C Nash 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.
>>
>> 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.
>>
>> JN
>>
>> #  Here's the script #
>> #
>> # Manoranjan Muthusamy <[hidden email]>
>> #
>>
>> 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 <[hidden email]> wrote:
>>>>
>>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] 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
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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.
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>>>
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>
> ---------------------------------------------------------------------------
> Jeff Newmiller                        The     .....       .....  Go Live...
> DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...
>                                        Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
> /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
> ---------------------------------------------------------------------------

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

David Winsemius

> On Jun 18, 2017, at 3:47 PM, J C Nash <[hidden email]> wrote:
>
> I've seen a number of problems like this over the years. The fact that the singular values of the Jacobian have
> a ration larger than the usual convergence tolerances can mean the codes stop well before the best fit. That is
> the "numerical analyst" view. David and Jeff have given geometric and statistical arguments. All views are useful,
> but it takes some time to sort them all out and make sense of the problem as a whole. Right now I'm getting ready
> to go to UseR!, so probably won't have time to spend sorting things out, though I'll have a go if I get time.
>
> David: Did you get a crash with the Mac version of nlsr? (You have the latest version I uploaded to CRAN.)

No. Everything proceeded without error.

> summary(myfitj)
$residuals
[1]  0.000000  0.229245  0.189927  0.130466  0.050909 -0.271920
attr(,"gradient")
              a      b
[1,] 0.0000e+00      0
[2,] 7.9014e-07  79784
[3,] 1.4207e-06 143370
[4,] 2.0498e-06 206741
[5,] 2.6774e-06 269898
[6,] 3.1473e-06 317126

$sigma
[1] 0.21341

$df
[1] 2 4

$cov.unscaled
   a  b
a NA NA
b NA NA

$param
    Estimate Std. Error t value Pr(>|t|)
a 1.4212e+07         NA      NA       NA
b 2.8129e-04         NA      NA       NA

$resname
[1] "myfitj"

$ssquares
[1] 0.18218

$nobs
[1] 6

$ct
[1] " " " "

$mt
[1] " " " "

$Sd
[1] 4.9301e+05 2.7070e-09

$gr
         [,1]
a -1.1226e-09
b -2.0297e-01

$jeval
[1] 525

$feval
[1] 656

attr(,"pkgname")
[1] "nlsr"
attr(,"class")
[1] "summary.nlsr"
> print(myfitj)
nlsr object: x
residual sumsquares =  0.18218  on  6 observations
    after  525    Jacobian and  656 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval  
a               14211701            NA         NA         NA  -1.123e-09      493005  
b            0.000281292            NA         NA         NA      -0.203   2.707e-09  

And (after considerable "verbosity":

> summary(myfitj2)
$residuals
[1]  0.000000  0.195004  0.149715  0.103277  0.055691 -0.230752
attr(,"gradient")
             ab       b
[1,] 0.00000000       0
[2,] 0.00016034  698106
[3,] 0.00028859 1256430
[4,] 0.00041682 1814611
[5,] 0.00054504 2372649
[6,] 0.00064119 2791083

$sigma
[1] 0.1785

$df
[1] 2 4

$cov.unscaled
   ab  b
ab NA NA
b  NA NA

$param
     Estimate Std. Error t value Pr(>|t|)
ab 6.9822e+04         NA      NA       NA
b  1.6035e-05         NA      NA       NA

$resname
[1] "myfitj2"

$ssquares
[1] 0.12746

$nobs
[1] 6

$ct
[1] " " " "

$mt
[1] " " " "

$Sd
[1] 4.3293e+06 6.2886e-08

$gr
          [,1]
ab -8.2085e-08
b  -2.6333e+02

$jeval
[1] 736

$feval
[1] 1024

attr(,"pkgname")
[1] "nlsr"
attr(,"class")
[1] "summary.nlsr"
> print(myfitj2)
nlsr object: x
residual sumsquares =  0.12746  on  6 observations
    after  736    Jacobian and  1024 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval  
ab               69821.8            NA         NA         NA  -8.208e-08     4329320  
b             1.6035e-05            NA         NA         NA      -263.3   6.289e-08  

The estimates were the same order of magnitude but off by a factor of 3 or so. I wondered whether some sort of scaling would have improved the estimates. I did try plotting the geometric picture I imagined using `persp`, but my R-fu was weak.

Best;
David.



> I don't have Mac, just Linux and a very ancient Win XP. I suspect the latter is too old to take seriously. I use Win-Builder for package checks. I don't think there's anything seriously wrong with nlsr functions, but there is
> more to be done to the vignettes and better manual pages are always desirable. It is the vignettes that are
> popping up warnings on CRAN checks, and I'll see if I can fix those.
>
> Cheers, JN
>
>
> On 2017-06-18 06:23 PM, Jeff Newmiller wrote:
>> I am not as expert as John, but I thought it worth pointing out that the variable substitution technique gives up one set of constraints for another (b=0 in this case). I also find that plots help me see what is going on, so here is my reproducible example (note inclusion of library calls for completeness). Note that NONE of the optimizers mentioned so far
>> appear to be finding the true best fit. The fact that myfun() yields 0 always if t=0 and that condition is within the data given seems likely to be part of the problem. I don't know how to resolve this... perhaps John will look at it again.
>> David: Your thinking makes fine sense if you are using a Monte Carlo or brute force solution, but the fact that it creates a discontinuity in the objective function will confuse any optimizer that uses analytic or numerically estimated slopes.
>> ##----------begin
>> library(minpack.lm)
>> library(ggplot2)
>> 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 ) {
>>   a * b * ( 1 - exp( -b * r * t ) )
>> }
>> objdta <- expand.grid( a = seq( 1000, 3000, by=20 )
>>                      , b = seq( -0.01, 1, 0.01 )
>>                      , rowidx = seq.int( nrow( mydata ) )
>>                      )
>> objdta[ , c( "y", "t" ) ] <- mydata[ objdta$rowidx
>>                                    , c( "y", "x" ) ]
>> objdta$tf <- factor( objdta$t )
>> objdta$myfun <- with( objdta
>>                     , myfun( a = a, b = b, r = 2, t = t )
>>                     )
>> objdtass <- aggregate( ( objdta$myfun - objdta$y )^2
>>                      , objdta[ , c( "a", "b" ) ]
>>                      , FUN = function( x )
>>                               sum( x, na.rm=TRUE )
>>                      )
>> objdtassmin <- objdtass[ which.min( objdtass$x ), ]
>> 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 )
>>               )
>> a <- as.vector( coef( myfit )[ "a" ] )
>> b <- as.vector( coef( myfit )[ "b" ] )
>> brks <- c( 500, 1e7, 2e7, 3e7, 4e7 )
>> ggplot( objdtass, aes( x=a, y=b, z = x, fill=x ) ) +
>>     geom_tile() +
>>     geom_contour( breaks= brks ) +
>>     geom_point( x=a, y=b, colour="red" ) +
>>     geom_point( x=objdtassmin$a
>>               , y=objdtassmin$b
>>               , colour="green" ) +
>>     scale_fill_continuous( name="SumSq", breaks = brks )
>> # Green point is brute-force solution
>> # Red point is optimizer solution for myfun
>> ##############
>> myfun2 <- function( a, log1ab, r, t ) {
>>   ab <- 1000 - exp( log1ab )
>>   ab * ( 1 - exp( -ab/a * r * t ) )
>> }
>> objdta$log1ab <- with( objdta, log( 1000 - a * b ) )
>> objdta$myfun2 <- with( objdta
>>                      , myfun2( a = a
>>                              , log1ab = log1ab
>>                              , r = 2
>>                              , t = t
>>                              )
>>                      )
>> objdtass2 <- aggregate( ( objdta$myfun2 - objdta$y )^2
>>                       , objdta[ , c( "a", "b" ) ]
>>                       , FUN = function( x )
>>                                if ( all( is.na( x ) ) ) NA
>>                                else sum( x, na.rm=TRUE )
>>                       )
>> objdtass2min <- objdtass2[ which.min( objdtass2$x ), ]
>> myfit2 <- nlsLM( y ~ myfun2( a, log1ab, r = 2, t = x )
>>                , data = mydata
>>                , start = list( a = 2000
>>                              , log1ab = 4.60517
>>                              )
>>                , lower = c( 1000, 0 )
>>                , upper = c( 3000, 8.0063 )
>>                )
>> a2 <- as.vector( coef( myfit2 )[ "a" ] )
>> b2 <- ( 1000
>>       - exp( as.vector( coef( myfit2 )[ "log1ab" ] ) )
>>       ) / a
>> brks <- c( 500, 1e6, 2e6, 3e6, 4e6 )
>> ggplot( objdtass2, aes( x=a, y=b, z = x, fill=x ) ) +
>>     geom_tile() +
>>     geom_contour( breaks = brks ) +
>>     geom_point( x=a2, y=b2, colour="red" ) +
>>     geom_point( x=objdtass2min$a
>>               , y=objdtass2min$b
>>               , colour="green" ) +
>>     scale_fill_continuous( name="SumSq", breaks = brks )
>> # Green point is brute-force solution
>> # Red point is optimizer solution for myfun2
>> ##----------end
>> On Sun, 18 Jun 2017, J C Nash 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.
>>>
>>> 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.
>>>
>>> JN
>>>
>>> #  Here's the script #
>>> #
>>> # Manoranjan Muthusamy <[hidden email]>
>>> #
>>>
>>> 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 <[hidden email]> wrote:
>>>>>
>>>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>>>>>
>>>>>> ______________________________________________
>>>>>> [hidden email] 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
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] 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.
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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.
>>>>
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>>>
>> ---------------------------------------------------------------------------
>> Jeff Newmiller                        The     .....       .....  Go Live...
>> DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...
>>                                       Live:   OO#.. Dead: OO#..  Playing
>> Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
>> /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
>> ---------------------------------------------------------------------------

David Winsemius
Alameda, CA, USA

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: R_using non linear regression with constraints

J C Nash
In reply to this post by Jeff Newmiller
I took another look at the problem.

Essentially there isn't enough information in the data to estimate the parameters. What there is
doesn't have the right shape for the model.

You can see this by plotting the data -- essentially a straight line.

Get the slope as follows:

    rr <- mean(mydata$y/mydata$x, na.rm=TRUE)

    rr

    mydata$yx <- mydata$y - rr*mydata$x

    with(mydata, plot(x, yx))

Shows that the curve is bending the wrong way for the model.

More data is needed, especially at larger values of x.

JN

On 2017-06-18 06:23 PM, Jeff Newmiller wrote:

> I am not as expert as John, but I thought it worth pointing out that the variable substitution technique gives up one
> set of constraints for another (b=0 in this case). I also find that plots help me see what is going on, so here is my
> reproducible example (note inclusion of library calls for completeness). Note that NONE of the optimizers mentioned so far
> appear to be finding the true best fit. The fact that myfun() yields 0 always if t=0 and that condition is within the
> data given seems likely to be part of the problem. I don't know how to resolve this... perhaps John will look at it again.
>
> David: Your thinking makes fine sense if you are using a Monte Carlo or brute force solution, but the fact that it
> creates a discontinuity in the objective function will confuse any optimizer that uses analytic or numerically estimated
> slopes.
>
> ##----------begin
> library(minpack.lm)
> library(ggplot2)
>
> 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 ) {
>    a * b * ( 1 - exp( -b * r * t ) )
> }
>
> objdta <- expand.grid( a = seq( 1000, 3000, by=20 )
>                       , b = seq( -0.01, 1, 0.01 )
>                       , rowidx = seq.int( nrow( mydata ) )
>                       )
> objdta[ , c( "y", "t" ) ] <- mydata[ objdta$rowidx
>                                     , c( "y", "x" ) ]
> objdta$tf <- factor( objdta$t )
> objdta$myfun <- with( objdta
>                      , myfun( a = a, b = b, r = 2, t = t )
>                      )
> objdtass <- aggregate( ( objdta$myfun - objdta$y )^2
>                       , objdta[ , c( "a", "b" ) ]
>                       , FUN = function( x )
>                                sum( x, na.rm=TRUE )
>                       )
> objdtassmin <- objdtass[ which.min( objdtass$x ), ]
>
> 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 )
>                )
> a <- as.vector( coef( myfit )[ "a" ] )
> b <- as.vector( coef( myfit )[ "b" ] )
>
> brks <- c( 500, 1e7, 2e7, 3e7, 4e7 )
> ggplot( objdtass, aes( x=a, y=b, z = x, fill=x ) ) +
>      geom_tile() +
>      geom_contour( breaks= brks ) +
>      geom_point( x=a, y=b, colour="red" ) +
>      geom_point( x=objdtassmin$a
>                , y=objdtassmin$b
>                , colour="green" ) +
>      scale_fill_continuous( name="SumSq", breaks = brks )
> # Green point is brute-force solution
> # Red point is optimizer solution for myfun
>
> ##############
>
> myfun2 <- function( a, log1ab, r, t ) {
>    ab <- 1000 - exp( log1ab )
>    ab * ( 1 - exp( -ab/a * r * t ) )
> }
>
> objdta$log1ab <- with( objdta, log( 1000 - a * b ) )
> objdta$myfun2 <- with( objdta
>                       , myfun2( a = a
>                               , log1ab = log1ab
>                               , r = 2
>                               , t = t
>                               )
>                       )
> objdtass2 <- aggregate( ( objdta$myfun2 - objdta$y )^2
>                        , objdta[ , c( "a", "b" ) ]
>                        , FUN = function( x )
>                                 if ( all( is.na( x ) ) ) NA
>                                 else sum( x, na.rm=TRUE )
>                        )
> objdtass2min <- objdtass2[ which.min( objdtass2$x ), ]
>
> myfit2 <- nlsLM( y ~ myfun2( a, log1ab, r = 2, t = x )
>                 , data = mydata
>                 , start = list( a = 2000
>                               , log1ab = 4.60517
>                               )
>                 , lower = c( 1000, 0 )
>                 , upper = c( 3000, 8.0063 )
>                 )
> a2 <- as.vector( coef( myfit2 )[ "a" ] )
> b2 <- ( 1000
>        - exp( as.vector( coef( myfit2 )[ "log1ab" ] ) )
>        ) / a
>
> brks <- c( 500, 1e6, 2e6, 3e6, 4e6 )
> ggplot( objdtass2, aes( x=a, y=b, z = x, fill=x ) ) +
>      geom_tile() +
>      geom_contour( breaks = brks ) +
>      geom_point( x=a2, y=b2, colour="red" ) +
>      geom_point( x=objdtass2min$a
>                , y=objdtass2min$b
>                , colour="green" ) +
>      scale_fill_continuous( name="SumSq", breaks = brks )
> # Green point is brute-force solution
> # Red point is optimizer solution for myfun2
>
> ##----------end
>
> On Sun, 18 Jun 2017, J C Nash 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.
>>
>> 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.
>>
>> JN
>>
>> #  Here's the script #
>> #
>> # Manoranjan Muthusamy <[hidden email]>
>> #
>>
>> 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 <[hidden email]> wrote:
>>>>
>>>>> On Jun 18, 2017, at 6:24 AM, Manoranjan Muthusamy <[hidden email]> 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]]
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] 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
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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.
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>>>
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>
> ---------------------------------------------------------------------------
> Jeff Newmiller                        The     .....       .....  Go Live...
> DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...
>                                        Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
> /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
> ---------------------------------------------------------------------------

______________________________________________
[hidden email] 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.