Finding starting values for the parameters using nls() or nls2()

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

Finding starting values for the parameters using nls() or nls2()

Pinglei Gao
Hi,

I have some data that i'm trying to fit a double exponential model: data.
Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
exponential is: exp (b0*exp (b1*x^th)).

 

I failed to guess the initial parameter values and then I learned a measure
to find starting values from Nonlinear Regression with R (pp. 25-27):

 

> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

+ Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}

> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
c(0.02),b1 = seq(0.01, 4, by = 0.01)))

> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
= grid.Disperse, algorithm = "brute-force")

> Disperse.m2a

Nonlinear regression model

  model: Retention ~ expFct(Area, b0, th, b1)

   data: cl

b0   th   b1

3.82 0.02 0.01

residual sum-of-squares: 13596

Number of iterations to convergence: 160000

Achieved convergence tolerance: NA

 

I got no error then I use the output as starting values to nls2 ():

> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
list(b0 = 3.82, b1 = 0.02, th = 0.01))

Error in (function (formula, data = parent.frame(), start, control =
nls.control(),  :

Singular gradient

 

Why? Did I do something wrong or misunderstand something?

 

Later, I found another measure from Modern Applied Statistics with S (pp.
216-217):

 

> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
negexp.SSival, parameters = c("b0", "b1", "th"),

+ template = function(x, b0, b1, th) {})

> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
T)

         b0          b1          th

   4.208763  144.205455 1035.324595

Error in qr.default(.swts * attr(rhs, "gradient")) :

 NA/NaN/Inf (arg1) can not be called when the external function is called.

 

Error happened again. How can I fix it? I am desperate.

 

Best regards,

 

Pinglei Gao

 


        [[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: Finding starting values for the parameters using nls() or nls2()

Andrew Robinson-6
Here are some things to try.  Maybe divide Area by 1000 and retention
by 100.  Try plotting the data and superimposing the line that
corresponds to the 'fit' from nls2.  See if you can correct it with
some careful guesses.

Getting suitable starting parameters for non-linear modeling is one of
the black arts of statistical fitting. Good luck!  And don't forget to
check for sensitivity.

Andrew

On 9 October 2016 at 22:21, Pinglei Gao <[hidden email]> wrote:

> Hi,
>
> I have some data that i'm trying to fit a double exponential model: data.
> Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
> 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>
> Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
> exponential is: exp (b0*exp (b1*x^th)).
>
>
>
> I failed to guess the initial parameter values and then I learned a measure
> to find starting values from Nonlinear Regression with R (pp. 25-27):
>
>
>
>> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>
> + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )
>
>> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
>
>> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
> c(0.02),b1 = seq(0.01, 4, by = 0.01)))
>
>> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
> = grid.Disperse, algorithm = "brute-force")
>
>> Disperse.m2a
>
> Nonlinear regression model
>
>   model: Retention ~ expFct(Area, b0, th, b1)
>
>    data: cl
>
> b0   th   b1
>
> 3.82 0.02 0.01
>
> residual sum-of-squares: 13596
>
> Number of iterations to convergence: 160000
>
> Achieved convergence tolerance: NA
>
>
>
> I got no error then I use the output as starting values to nls2 ():
>
>> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
> list(b0 = 3.82, b1 = 0.02, th = 0.01))
>
> Error in (function (formula, data = parent.frame(), start, control =
> nls.control(),  :
>
> Singular gradient
>
>
>
> Why? Did I do something wrong or misunderstand something?
>
>
>
> Later, I found another measure from Modern Applied Statistics with S (pp.
> 216-217):
>
>
>
>> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
> negexp.SSival, parameters = c("b0", "b1", "th"),
>
> + template = function(x, b0, b1, th) {})
>
>> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
> T)
>
>          b0          b1          th
>
>    4.208763  144.205455 1035.324595
>
> Error in qr.default(.swts * attr(rhs, "gradient")) :
>
>  NA/NaN/Inf (arg1) can not be called when the external function is called.
>
>
>
> Error happened again. How can I fix it? I am desperate.
>
>
>
> Best regards,
>
>
>
> Pinglei Gao
>
>
>
>
>         [[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.



--
Andrew Robinson
Deputy Director, CEBRA, School of Biosciences
Reader & Associate Professor in Applied Statistics  Tel: (+61) 0403 138 955
School of Mathematics and Statistics                        Fax: +61-3-8344 4599
University of Melbourne, VIC 3010 Australia
Email: [hidden email]
Website: http://www.ms.unimelb.edu.au/~andrewpr

MSME: http://www.crcpress.com/product/isbn/9781439858028
FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
SPuR: http://www.ms.unimelb.edu.au/spuRs/

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

Bert Gunter-2
Well... (inline -- and I hope this isn't homework!)




On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
<[hidden email]> wrote:
> Here are some things to try.  Maybe divide Area by 1000 and retention
> by 100.  Try plotting the data and superimposing the line that
> corresponds to the 'fit' from nls2.  See if you can correct it with
> some careful guesses.
>
> Getting suitable starting parameters for non-linear modeling is one of
> the black arts of statistical fitting. ...
>
> Andrew

True. But it's usually worthwhile thinking about the math a bit before guessing.

Note that the model can be linearized to:

log(log(Retention)) = b0 + b1*Area^th

So a plot of log(log(Retention)) vs Area may be informative and useful
for finding starting values. e.g., for a grid of th's, do linear
regression fits .

However, when I look at that plot, it seems pretty linear with a
negative slope. This suggests that you may have an overparametrization
problem . i.e. fix th =1 and use the b0 and b1 from the above
regression for starting values.

Do note that this strategy isn't foolproof, as it ignores that the
error term is additive in the above transformed metric, rather than
the original. This can sometimes mislead. But this is just a
heuristic.

Cheers,
Bert







>
> On 9 October 2016 at 22:21, Pinglei Gao <[hidden email]> wrote:
>> Hi,
>>
>> I have some data that i'm trying to fit a double exponential model: data.
>> Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
>> 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>
>> Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
>> exponential is: exp (b0*exp (b1*x^th)).
>>
>>
>>
>> I failed to guess the initial parameter values and then I learned a measure
>> to find starting values from Nonlinear Regression with R (pp. 25-27):
>>
>>
>>
>>> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
>> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>
>> + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )
>>
>>> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
>>
>>> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
>> c(0.02),b1 = seq(0.01, 4, by = 0.01)))
>>
>>> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
>> = grid.Disperse, algorithm = "brute-force")
>>
>>> Disperse.m2a
>>
>> Nonlinear regression model
>>
>>   model: Retention ~ expFct(Area, b0, th, b1)
>>
>>    data: cl
>>
>> b0   th   b1
>>
>> 3.82 0.02 0.01
>>
>> residual sum-of-squares: 13596
>>
>> Number of iterations to convergence: 160000
>>
>> Achieved convergence tolerance: NA
>>
>>
>>
>> I got no error then I use the output as starting values to nls2 ():
>>
>>> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
>> list(b0 = 3.82, b1 = 0.02, th = 0.01))
>>
>> Error in (function (formula, data = parent.frame(), start, control =
>> nls.control(),  :
>>
>> Singular gradient
>>
>>
>>
>> Why? Did I do something wrong or misunderstand something?
>>
>>
>>
>> Later, I found another measure from Modern Applied Statistics with S (pp.
>> 216-217):
>>
>>
>>
>>> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
>> negexp.SSival, parameters = c("b0", "b1", "th"),
>>
>> + template = function(x, b0, b1, th) {})
>>
>>> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
>> T)
>>
>>          b0          b1          th
>>
>>    4.208763  144.205455 1035.324595
>>
>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>>
>>  NA/NaN/Inf (arg1) can not be called when the external function is called.
>>
>>
>>
>> Error happened again. How can I fix it? I am desperate.
>>
>>
>>
>> Best regards,
>>
>>
>>
>> Pinglei Gao
>>
>>
>>
>>
>>         [[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.
>
>
>
> --
> Andrew Robinson
> Deputy Director, CEBRA, School of Biosciences
> Reader & Associate Professor in Applied Statistics  Tel: (+61) 0403 138 955
> School of Mathematics and Statistics                        Fax: +61-3-8344 4599
> University of Melbourne, VIC 3010 Australia
> Email: [hidden email]
> Website: http://www.ms.unimelb.edu.au/~andrewpr
>
> MSME: http://www.crcpress.com/product/isbn/9781439858028
> FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
> SPuR: http://www.ms.unimelb.edu.au/spuRs/
>
> ______________________________________________
> [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: Finding starting values for the parameters using nls() or nls2()

Peter Dalgaard-2

> On 10 Oct 2016, at 00:40 , Bert Gunter <[hidden email]> wrote:
>
> Well... (inline -- and I hope this isn't homework!)
>

Pretty much same as I thought.

Fixing th=0.02 in the grid search looks wrong. Bert's plot is pretty linear, so th=1 is a good guesstimate. There's a slight curvature but to reduce it, you would increase th, not decrease it. Running the regression, as Bert suggests, indicates that b0=5.16 and b1= -0.00024 could work as reasonable starting values. Notice that the grid search had "b1 = seq(0.01, 4, by = 0.01)" which is wrong in both sign and scale.

Andrew's suggestion of dividing Retention by 100 is tempting, since it looks like a percentage, but that would make all Y values less than 1 and the double exponential function as written has values that are always bigger than 1. (It is conceivable that the model itself is wrong, though. E.g. it could be that Retention on a scale from 0 to 1 could be modeled as exp(-something), but we really have no idea of the context here.)

(If this was in fact homework, you should now go and write a proper SelfStart initializer routine for this model. Even if it isn't homework, you do need to study the text again, because you have clearly not understood how self-starting models work.)

-pd

>
>
>
> On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
> <[hidden email]> wrote:
>> Here are some things to try.  Maybe divide Area by 1000 and retention
>> by 100.  Try plotting the data and superimposing the line that
>> corresponds to the 'fit' from nls2.  See if you can correct it with
>> some careful guesses.
>>
>> Getting suitable starting parameters for non-linear modeling is one of
>> the black arts of statistical fitting. ...
>>
>> Andrew
>
> True. But it's usually worthwhile thinking about the math a bit before guessing.
>
> Note that the model can be linearized to:
>
> log(log(Retention)) = b0 + b1*Area^th
>
> So a plot of log(log(Retention)) vs Area may be informative and useful
> for finding starting values. e.g., for a grid of th's, do linear
> regression fits .
>
> However, when I look at that plot, it seems pretty linear with a
> negative slope. This suggests that you may have an overparametrization
> problem . i.e. fix th =1 and use the b0 and b1 from the above
> regression for starting values.
>
> Do note that this strategy isn't foolproof, as it ignores that the
> error term is additive in the above transformed metric, rather than
> the original. This can sometimes mislead. But this is just a
> heuristic.
>
> Cheers,
> Bert
>
>
>
>
>
>
>
>>
>> On 9 October 2016 at 22:21, Pinglei Gao <[hidden email]> wrote:
>>> Hi,
>>>
>>> I have some data that i'm trying to fit a double exponential model: data.
>>> Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
>>> 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>>
>>> Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
>>> exponential is: exp (b0*exp (b1*x^th)).
>>>
>>>
>>>
>>> I failed to guess the initial parameter values and then I learned a measure
>>> to find starting values from Nonlinear Regression with R (pp. 25-27):
>>>
>>>
>>>
>>>> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
>>> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>>
>>> + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )
>>>
>>>> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
>>>
>>>> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
>>> c(0.02),b1 = seq(0.01, 4, by = 0.01)))
>>>
>>>> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
>>> = grid.Disperse, algorithm = "brute-force")
>>>
>>>> Disperse.m2a
>>>
>>> Nonlinear regression model
>>>
>>>  model: Retention ~ expFct(Area, b0, th, b1)
>>>
>>>   data: cl
>>>
>>> b0   th   b1
>>>
>>> 3.82 0.02 0.01
>>>
>>> residual sum-of-squares: 13596
>>>
>>> Number of iterations to convergence: 160000
>>>
>>> Achieved convergence tolerance: NA
>>>
>>>
>>>
>>> I got no error then I use the output as starting values to nls2 ():
>>>
>>>> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
>>> list(b0 = 3.82, b1 = 0.02, th = 0.01))
>>>
>>> Error in (function (formula, data = parent.frame(), start, control =
>>> nls.control(),  :
>>>
>>> Singular gradient
>>>
>>>
>>>
>>> Why? Did I do something wrong or misunderstand something?
>>>
>>>
>>>
>>> Later, I found another measure from Modern Applied Statistics with S (pp.
>>> 216-217):
>>>
>>>
>>>
>>>> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
>>> negexp.SSival, parameters = c("b0", "b1", "th"),
>>>
>>> + template = function(x, b0, b1, th) {})
>>>
>>>> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
>>> T)
>>>
>>>         b0          b1          th
>>>
>>>   4.208763  144.205455 1035.324595
>>>
>>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>>>
>>> NA/NaN/Inf (arg1) can not be called when the external function is called.
>>>
>>>
>>>
>>> Error happened again. How can I fix it? I am desperate.
>>>
>>>
>>>
>>> Best regards,
>>>
>>>
>>>
>>> Pinglei Gao
>>>
>>>
>>>
>>>
>>>        [[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.
>>
>>
>>
>> --
>> Andrew Robinson
>> Deputy Director, CEBRA, School of Biosciences
>> Reader & Associate Professor in Applied Statistics  Tel: (+61) 0403 138 955
>> School of Mathematics and Statistics                        Fax: +61-3-8344 4599
>> University of Melbourne, VIC 3010 Australia
>> Email: [hidden email]
>> Website: http://www.ms.unimelb.edu.au/~andrewpr
>>
>> MSME: http://www.crcpress.com/product/isbn/9781439858028
>> FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
>> SPuR: http://www.ms.unimelb.edu.au/spuRs/
>>
>> ______________________________________________
>> [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.

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

J C Nash
I didn't try very hard, but got a solution from .1, 1, .1 with nlxb() from nlmrt. It took a lot
of iterations and looks to be pretty ill-conditioned. Note nlmrt uses analytic derivatives if it
can, and a Marquardt method. It is designed to be a pit bull -- tenacious, not fast.

I'm working on a replacement for this and nls(), but it will be a while. However, I welcome
short scripts like this as tests. My script below

Best, JN

cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
expf2 <- "Retention ~ exp(b0*exp(b1*Area^th))"

# grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
#c(0.02),b1 = seq(0.01, 4, by = 0.01)))

#> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
#= grid.Disperse, algorithm = "brute-force")

# Disperse.m2a

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)



On 16-10-09 07:40 PM, peter dalgaard wrote:

>
>> On 10 Oct 2016, at 00:40 , Bert Gunter <[hidden email]> wrote:
>>
>> Well... (inline -- and I hope this isn't homework!)
>>
>
> Pretty much same as I thought.
>
> Fixing th=0.02 in the grid search looks wrong. Bert's plot is pretty linear, so th=1 is a good guesstimate. There's a slight curvature but to reduce it, you would increase th, not decrease it. Running the regression, as Bert suggests, indicates that b0=5.16 and b1= -0.00024 could work as reasonable starting values. Notice that the grid search had "b1 = seq(0.01, 4, by = 0.01)" which is wrong in both sign and scale.
>
> Andrew's suggestion of dividing Retention by 100 is tempting, since it looks like a percentage, but that would make all Y values less than 1 and the double exponential function as written has values that are always bigger than 1. (It is conceivable that the model itself is wrong, though. E.g. it could be that Retention on a scale from 0 to 1 could be modeled as exp(-something), but we really have no idea of the context here.)
>
> (If this was in fact homework, you should now go and write a proper SelfStart initializer routine for this model. Even if it isn't homework, you do need to study the text again, because you have clearly not understood how self-starting models work.)
>
> -pd
>
>>
>>
>>
>> On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
>> <[hidden email]> wrote:
>>> Here are some things to try.  Maybe divide Area by 1000 and retention
>>> by 100.  Try plotting the data and superimposing the line that
>>> corresponds to the 'fit' from nls2.  See if you can correct it with
>>> some careful guesses.
>>>
>>> Getting suitable starting parameters for non-linear modeling is one of
>>> the black arts of statistical fitting. ...
>>>
>>> Andrew
>>
>> True. But it's usually worthwhile thinking about the math a bit before guessing.
>>
>> Note that the model can be linearized to:
>>
>> log(log(Retention)) = b0 + b1*Area^th
>>
>> So a plot of log(log(Retention)) vs Area may be informative and useful
>> for finding starting values. e.g., for a grid of th's, do linear
>> regression fits .
>>
>> However, when I look at that plot, it seems pretty linear with a
>> negative slope. This suggests that you may have an overparametrization
>> problem . i.e. fix th =1 and use the b0 and b1 from the above
>> regression for starting values.
>>
>> Do note that this strategy isn't foolproof, as it ignores that the
>> error term is additive in the above transformed metric, rather than
>> the original. This can sometimes mislead. But this is just a
>> heuristic.
>>
>> Cheers,
>> Bert
>>
>>
>>
>>
>>
>>
>>
>>>
>>> On 9 October 2016 at 22:21, Pinglei Gao <[hidden email]> wrote:
>>>> Hi,
>>>>
>>>> I have some data that i'm trying to fit a double exponential model: data.
>>>> Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
>>>> 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>>>
>>>> Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>>>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
>>>> exponential is: exp (b0*exp (b1*x^th)).
>>>>
>>>>
>>>>
>>>> I failed to guess the initial parameter values and then I learned a measure
>>>> to find starting values from Nonlinear Regression with R (pp. 25-27):
>>>>
>>>>
>>>>
>>>>> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
>>>> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
>>>>
>>>> + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
>>>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )
>>>>
>>>>> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
>>>>
>>>>> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
>>>> c(0.02),b1 = seq(0.01, 4, by = 0.01)))
>>>>
>>>>> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
>>>> = grid.Disperse, algorithm = "brute-force")
>>>>
>>>>> Disperse.m2a
>>>>
>>>> Nonlinear regression model
>>>>
>>>>  model: Retention ~ expFct(Area, b0, th, b1)
>>>>
>>>>   data: cl
>>>>
>>>> b0   th   b1
>>>>
>>>> 3.82 0.02 0.01
>>>>
>>>> residual sum-of-squares: 13596
>>>>
>>>> Number of iterations to convergence: 160000
>>>>
>>>> Achieved convergence tolerance: NA
>>>>
>>>>
>>>>
>>>> I got no error then I use the output as starting values to nls2 ():
>>>>
>>>>> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
>>>> list(b0 = 3.82, b1 = 0.02, th = 0.01))
>>>>
>>>> Error in (function (formula, data = parent.frame(), start, control =
>>>> nls.control(),  :
>>>>
>>>> Singular gradient
>>>>
>>>>
>>>>
>>>> Why? Did I do something wrong or misunderstand something?
>>>>
>>>>
>>>>
>>>> Later, I found another measure from Modern Applied Statistics with S (pp.
>>>> 216-217):
>>>>
>>>>
>>>>
>>>>> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
>>>> negexp.SSival, parameters = c("b0", "b1", "th"),
>>>>
>>>> + template = function(x, b0, b1, th) {})
>>>>
>>>>> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
>>>> T)
>>>>
>>>>         b0          b1          th
>>>>
>>>>   4.208763  144.205455 1035.324595
>>>>
>>>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>>>>
>>>> NA/NaN/Inf (arg1) can not be called when the external function is called.
>>>>
>>>>
>>>>
>>>> Error happened again. How can I fix it? I am desperate.
>>>>
>>>>
>>>>
>>>> Best regards,
>>>>
>>>>
>>>>
>>>> Pinglei Gao
>>>>
>>>>
>>>>
>>>>
>>>>        [[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.
>>>
>>>
>>>
>>> --
>>> Andrew Robinson
>>> Deputy Director, CEBRA, School of Biosciences
>>> Reader & Associate Professor in Applied Statistics  Tel: (+61) 0403 138 955
>>> School of Mathematics and Statistics                        Fax: +61-3-8344 4599
>>> University of Melbourne, VIC 3010 Australia
>>> Email: [hidden email]
>>> Website: http://www.ms.unimelb.edu.au/~andrewpr
>>>
>>> MSME: http://www.crcpress.com/product/isbn/9781439858028
>>> FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
>>> SPuR: http://www.ms.unimelb.edu.au/spuRs/
>>>
>>> ______________________________________________
>>> [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: Finding starting values for the parameters using nls() or nls2()

Gabor Grothendieck
In reply to this post by Pinglei Gao
If you are not tied to that model the SSasymp() model in R could be
considered and is easy to fit:

    # to plot points in order
    o <- order(cl$Area)
    cl.o <- cl[o, ]

    fm <- nls(Retention ~ SSasymp(Area, Asym, R0, lrc), cl.o)
    summary(fm)

    plot(Retention ~ Area, cl.o)
    lines(fitted(fm) ~ Area, cl.o, col = "red")

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

J C Nash
Despite nlmrt "solving" the OP's problem, I think Gabor's suggestion likely gives a more sensible approach
to the underlying modelling problem.

It is, of course, sometimes important to fit a particular model, in which case nls2 and nlmrt are set up to grind away.
And hopefully the follow-up to nlmrt I'm working on will have enough capability in getting analytic derivatives
to work for a wider class of models. Note that functional approaches in nlmrt and minpack.lm allow users to
provide derivatives. Too many users think numerical approximations are a panacea, but my experience is that
most problems benefit from very accurate derivatives, of which analytic expressions are generally the best.

JN

On 16-10-09 09:24 PM, Gabor Grothendieck wrote:

> If you are not tied to that model the SSasymp() model in R could be
> considered and is easy to fit:
>
>     # to plot points in order
>     o <- order(cl$Area)
>     cl.o <- cl[o, ]
>
>     fm <- nls(Retention ~ SSasymp(Area, Asym, R0, lrc), cl.o)
>     summary(fm)
>
>     plot(Retention ~ Area, cl.o)
>     lines(fitted(fm) ~ Area, cl.o, col = "red")
>
> ______________________________________________
> [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
|

答复: Finding starting values for the parameters using nls() or nls2()

Pinglei Gao
In reply to this post by Peter Dalgaard-2
Thanks very much for taking time on this. Your assistances are very much
appreciated. But, I am afraid that I still have a question to bother you.

I am working on a paper about weed seeds dispersal with harvest machine. I
found three general models for seed dispersal and retention after a review
of relevant literature. All models were optimized using nonlinear least
squares via the nls function in the statistical package R. The model that
best described the data will be determined by comparing Akaike Information
Criterion (AIC) values and the model with the lowest AIC score will be
selected.

The first general model incorporated simple exponential and power
exponential functions, its starting value was easily to be found. But, I am
stuck with model 2 (which was mentioned previously) and model 3 with the
form:  Retention = (b0*Area^th+1)^b1. The model 3 is totally different to
others. I tried the measures that you were mentioned. But I still can’t
find suitable starting values because of my limited knowledge. I hope you
can do me the favor again. I can send the draft to you when I finished the
paper, if it is necessary. Maybe you can give me some constructive
suggestion about statistic and model construction and I can name you as a
coauthor for your contributions.

Best,

Pinglei Gao

-----邮件原件-----
发件人: peter dalgaard [mailto:[hidden email]]
发送时间: 2016年10月10日 7:41
收件人: Pinglei Gao
抄送: Andrew Robinson; R help ([hidden email]); Bert Gunter
主题: Re: [R] Finding starting values for the parameters using nls() or
nls2()


> On 10 Oct 2016, at 00:40 , Bert Gunter <[hidden email]> wrote:
>
> Well... (inline -- and I hope this isn't homework!)
>

Pretty much same as I thought.

Fixing th=0.02 in the grid search looks wrong. Bert's plot is pretty linear,
so th=1 is a good guesstimate. There's a slight curvature but to reduce it,
you would increase th, not decrease it. Running the regression, as Bert
suggests, indicates that b0=5.16 and b1= -0.00024 could work as reasonable
starting values. Notice that the grid search had "b1 = seq(0.01, 4, by =
0.01)" which is wrong in both sign and scale.

Andrew's suggestion of dividing Retention by 100 is tempting, since it looks
like a percentage, but that would make all Y values less than 1 and the
double exponential function as written has values that are always bigger
than 1. (It is conceivable that the model itself is wrong, though. E.g. it
could be that Retention on a scale from 0 to 1 could be modeled as
exp(-something), but we really have no idea of the context here.)

(If this was in fact homework, you should now go and write a proper
SelfStart initializer routine for this model. Even if it isn't homework, you
do need to study the text again, because you have clearly not understood how
self-starting models work.)

-pd

>
>
>
> On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
> <[hidden email]> wrote:
>> Here are some things to try.  Maybe divide Area by 1000 and retention
>> by 100.  Try plotting the data and superimposing the line that
>> corresponds to the 'fit' from nls2.  See if you can correct it with
>> some careful guesses.
>>
>> Getting suitable starting parameters for non-linear modeling is one
>> of the black arts of statistical fitting. ...
>>
>> Andrew
>
> True. But it's usually worthwhile thinking about the math a bit before
guessing.

>
> Note that the model can be linearized to:
>
> log(log(Retention)) = b0 + b1*Area^th
>
> So a plot of log(log(Retention)) vs Area may be informative and useful
> for finding starting values. e.g., for a grid of th's, do linear
> regression fits .
>
> However, when I look at that plot, it seems pretty linear with a
> negative slope. This suggests that you may have an overparametrization
> problem . i.e. fix th =1 and use the b0 and b1 from the above
> regression for starting values.
>
> Do note that this strategy isn't foolproof, as it ignores that the
> error term is additive in the above transformed metric, rather than
> the original. This can sometimes mislead. But this is just a
> heuristic.
>
> Cheers,
> Bert
>
>
>
>
>
>
>
>>
>> On 9 October 2016 at 22:21, Pinglei Gao <[hidden email]> wrote:
>>> Hi,
>>>
>>> I have some data that i'm trying to fit a double exponential model:
data.

>>> Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
>>> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74,
>>> 2629.2),
>>>
>>> Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04,
>>> 23.46, 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of
>>> the double exponential is: exp (b0*exp (b1*x^th)).
>>>
>>>
>>>
>>> I failed to guess the initial parameter values and then I learned a
>>> measure to find starting values from Nonlinear Regression with R (pp.
25-27):

>>>
>>>
>>>
>>>> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46,
>>>> 524.91,
>>> 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74,
>>> 2629.2),
>>>
>>> + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24,
>>> + 33.04, 23.46,
>>> 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )
>>>
>>>> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
>>>
>>>> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
>>> c(0.02),b1 = seq(0.01, 4, by = 0.01)))
>>>
>>>> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl,
>>>> start
>>> = grid.Disperse, algorithm = "brute-force")
>>>
>>>> Disperse.m2a
>>>
>>> Nonlinear regression model
>>>
>>>  model: Retention ~ expFct(Area, b0, th, b1)
>>>
>>>   data: cl
>>>
>>> b0   th   b1
>>>
>>> 3.82 0.02 0.01
>>>
>>> residual sum-of-squares: 13596
>>>
>>> Number of iterations to convergence: 160000
>>>
>>> Achieved convergence tolerance: NA
>>>
>>>
>>>
>>> I got no error then I use the output as starting values to nls2 ():
>>>
>>>> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl,
>>>> start =
>>> list(b0 = 3.82, b1 = 0.02, th = 0.01))
>>>
>>> Error in (function (formula, data = parent.frame(), start, control =
>>> nls.control(),  :
>>>
>>> Singular gradient
>>>
>>>
>>>
>>> Why? Did I do something wrong or misunderstand something?
>>>
>>>
>>>
>>> Later, I found another measure from Modern Applied Statistics with S
(pp.

>>> 216-217):
>>>
>>>
>>>
>>>> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
>>> negexp.SSival, parameters = c("b0", "b1", "th"),
>>>
>>> + template = function(x, b0, b1, th) {})
>>>
>>>> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl,
>>>> trace =
>>> T)
>>>
>>>         b0          b1          th
>>>
>>>   4.208763  144.205455 1035.324595
>>>
>>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>>>
>>> NA/NaN/Inf (arg1) can not be called when the external function is
called.

>>>
>>>
>>>
>>> Error happened again. How can I fix it? I am desperate.
>>>
>>>
>>>
>>> Best regards,
>>>
>>>
>>>
>>> Pinglei Gao
>>>
>>>
>>>
>>>
>>>        [[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.
>>
>>
>>
>> --
>> Andrew Robinson
>> Deputy Director, CEBRA, School of Biosciences Reader & Associate
>> Professor in Applied Statistics  Tel: (+61) 0403 138 955
>> School of Mathematics and Statistics                        Fax:
+61-3-8344 4599

>> University of Melbourne, VIC 3010 Australia
>> Email: [hidden email]
>> Website: http://www.ms.unimelb.edu.au/~andrewpr
>>
>> MSME: http://www.crcpress.com/product/isbn/9781439858028
>> FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
>> SPuR: http://www.ms.unimelb.edu.au/spuRs/
>>
>> ______________________________________________
>> [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.

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000
Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

Pinglei Gao
In reply to this post by Pinglei Gao
Thanks very much for your kindness help. I run your script then came out
lots of outputs and I also studied the solution you posted. Forgive my
ignorance, I still can't find the suitable starting values. Did I
misunderstand something?

Best,

Pinglei Gao

-----邮件原件-----
发件人: ProfJCNash [mailto:[hidden email]]
发送时间: 2016年10月10日 10:41
收件人: Gabor Grothendieck; Pinglei Gao
主题: Re: [R] Finding starting values for the parameters using nls() or
nls2()

I forgot to post the "solution" found by nlmrt:

nlmrt class object: x
residual sumsquares =  1086.8  on  15 observations
    after  5001    Jacobian and  6991 function evaluations
  name            coeff          SE       tstat      pval      gradient
JSingval
b0            5.3274e-14            NA         NA         NA  -6.614e+13
1.735e+16
b1               33.5574            NA         NA         NA      -3.466
11518
th           -0.00721203            NA         NA         NA      -740.8
0.004635


Note the singular values -- this is the worst SV(max)/SV(min) ratio I've
observed!

JN

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

J C Nash
The key lines are

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)

Thus I started with .1 1 and .1. The "solution" from nlxb, which is using analytic derivatives
and a very aggressive Marquardt code to keep trying even in bad situations, was
as you included below. Note that the singular values of the Jacobian are given (they are
recorded on the same table as the parameters, but do NOT correspond to the parameters.
The placement was simply a tidy place to put these numbers.)

The ratio of these sv's is 1.735e+16/0.004635 or approx 4E+18, so the condition number
of the traditional Gauss Newton approach is about 1E+37. Not a nice problem!

You probably should reformulate.

JN




On 16-10-10 10:41 AM, Pinglei Gao wrote:

> Thanks very much for your kindness help. I run your script then came out
> lots of outputs and I also studied the solution you posted. Forgive my
> ignorance, I still can't find the suitable starting values. Did I
> misunderstand something?
>
> Best,
>
> Pinglei Gao
>
> -----邮件原件-----
> 发件人: ProfJCNash [mailto:[hidden email]]
> 发送时间: 2016年10月10日 10:41
> 收件人: Gabor Grothendieck; Pinglei Gao
> 主题: Re: [R] Finding starting values for the parameters using nls() or
> nls2()
>
> I forgot to post the "solution" found by nlmrt:
>
> nlmrt class object: x
> residual sumsquares =  1086.8  on  15 observations
>     after  5001    Jacobian and  6991 function evaluations
>   name            coeff          SE       tstat      pval      gradient
> JSingval
> b0            5.3274e-14            NA         NA         NA  -6.614e+13
> 1.735e+16
> b1               33.5574            NA         NA         NA      -3.466
> 11518
> th           -0.00721203            NA         NA         NA      -740.8
> 0.004635
>
>
> Note the singular values -- this is the worst SV(max)/SV(min) ratio I've
> observed!
>
> JN
>
>
>

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

dave fournier-2
In reply to this post by Pinglei Gao
Actually this converges very nicely if you use these starting values
that I obtained with
AD Model Builder

        th     9.1180e-01
        b0    5.2104e+00
        b1   -4.6725e-04

The R result looks like

nls.m2
Nonlinear regression model
   model: Retention ~ expFct(Area, b0, b1, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
         b0         b1         th
  5.2104466 -0.0004672  0.9118029
  residual sum-of-squares: 686.8

Number of iterations to convergence: 1
Achieved convergence tolerance: 1.75e-06

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

dave fournier-2
In reply to this post by Pinglei Gao
I've spent quite a bit of time trying to convince people on various
lists that the solution to these kinds of
problems lies in the stable parameterization of the model.  I write the
solutions in AD Model Builder because it
is easy.  But R people are generally stuck in R (or mired) so the
message somehow always gets lost.

  I thought I would bite the bullet and figure out how to do this in R.
It turns out that one can fit this model
  with only 6 function evaluations using nls2.   I apologize in advance
for my execrable R code. But it does do
the job.

The data were fit using 3 calls to nls2. For the first call only the
parameter th was estimated. This converged
in 4 function evaluations. For the second call all three of the
parameters were estimated
(but the other 2 parameters are different from b0, b1) .This converged
to the solution in four function evaluations.

The final call to nls2 uses the converged values to calculate b0,b1 and
theta and starts from there.
As you can see the model is already converged. This final call to nls2
is used to get the standard error
estimates for the parameters.

 >  nls.m1
Nonlinear regression model
   model: Retention ~ expFct1(Area, y1, yn, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
     th
0.9862
  residual sum-of-squares: 730

Number of iterations to convergence: 2
Achieved convergence tolerance: 1.616e-06

 >  nls.m2
Nonlinear regression model
   model: Retention ~ expFct2(Area, y1, yn, c, d, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
      c      d     th
0.9716 1.1010 0.9118
  residual sum-of-squares: 686.8

Number of iterations to convergence: 4
Achieved convergence tolerance: 2.398e-06

 >  nls.m
Nonlinear regression model
   model: Retention ~ expFct(Area, b0, b1, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
         b0         b1         th
  5.2104452 -0.0004672  0.9118040
  residual sum-of-squares: 686.8

Number of iterations to convergence: 0
Achieved convergence tolerance: 1.384e-06

Parameters:
      Estimate Std. Error t value Pr(>|t|)
b0  5.2104452  0.4999594  10.422 2.29e-07 ***
b1 -0.0004672  0.0013527  -0.345   0.7358
th  0.9118040  0.3583575   2.544   0.0257 *

So what is the stable parameterization for this model.  To simplify let
x be the independent variable and y be the
  dependent variable and write t instead of th  So the model is

             y = exp(b0*exp(b1*x^t) (1)

Now it would be nice to reorder the x's so that they are monotone
increasing or decreasing, but in this case the
first x is almost the largest and the last x is almost the smallest
(slight reach) so they will do. Call them x1 and xn.
The new parameters of the model are the predicted y's for x1 and xn.  
Call them y1 and yn.  Note that y1 and yn
are parameters, not observations.


The data were fit using 3 calls to nls2. For the first call only the
parameter th was estimated. this converged
in 4 function evaluestions. for the second call all three of the
parameters were estimated
(but the other 2 parameters are different from b0, b1) .This converged
to the solution in four function evaluations.

The final call to nls2 uses the converged values to calculate b0,b1 and
theta and starts from there.
as you can see the model is already converged. this final call to nls2
is used to get the standard error
estimates for the parameters.

 >  nls.m1
Nonlinear regression model
   model: Retention ~ expFct1(Area, y1, yn, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
     th
0.9862
  residual sum-of-squares: 730

Number of iterations to convergence: 2
Achieved convergence tolerance: 1.616e-06

 >  nls.m2
Nonlinear regression model
   model: Retention ~ expFct2(Area, y1, yn, c, d, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
      c      d     th
0.9716 1.1010 0.9118
  residual sum-of-squares: 686.8

Number of iterations to convergence: 4
Achieved convergence tolerance: 2.398e-06

 >  nls.m
Nonlinear regression model
   model: Retention ~ expFct(Area, b0, b1, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,
2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,
1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,
18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,
15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
         b0         b1         th
  5.2104452 -0.0004672  0.9118040
  residual sum-of-squares: 686.8

Number of iterations to convergence: 0
Achieved convergence tolerance: 1.384e-06

Parameters:
      Estimate Std. Error t value Pr(>|t|)
b0  5.2104452  0.4999594  10.422 2.29e-07 ***
b1 -0.0004672  0.0013527  -0.345   0.7358
th  0.9118040  0.3583575   2.544   0.0257 *

So what is the stable parameterization for this model.  To simplify letx
be the independent variable and y be the
  dependent variable and write t insted of th  So the model is

             y = exp(b0*exp(b1*x^t) (1)

Now it would be nice to reorder the x's so that they are monotone
increasing or decreasing, but in this case the
first x is almost the largest and the last x is almost the smallest
(slight reach) so they will do. Call them x1 and xn.
The new parameters of the model are the predicted y's for x1 and xn.  
Call them y1 and yn.  Note that y1 and yn
are parameters, not observations.

              y1 = exp(b0*exp(b1*x1^t)
(2)
              yn = exp(b0*exp(b1*xn^t)

One can solve for b1 and b0 in terms of y1 and yn.

             b1=log(log(y1)/log(yn))/(x1^t-xn^t);
(3)
             b0=log(y1)*exp(-b1*x1^t);


To use this  we do the fitting of the model in two stages. In the first
stage we use the obvious
  estimates of y[1] for y1 and y[n] for yn and fix these values and
estimate t using the obvious
starting value of 1 for t.

In the second stage we set

                y1=c*y[1]

                yn=d*y[n]

and estimate   c, d, and t using the obvious starting values of 1 for c
and d.
That's it. this converges as stated in 6 function evaluations.  This
technique works for
may curves such as 3,4,5 parameter logistic double exponential
vonBertalanffy
(where Schnute and I first discovered it in the early 1980's before
anyone was born).

One caveat  is that usually not all values of y1 and yn are
permissible.  For example in this case if b>0
then y1 and yn are both >1. To make this routine bulletproof one needs
to impose these kinds of
constraints.  But this was not needed here.

Here is my  R code for this model.

library("nls2")
cl=data.frame(Area=c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05,
1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

  y1 <<- cl$Retention[1];
  yn <<- cl$Retention[length(cl$Retention)]

  expFct1 <- function(x, y1, yn,th)
  {
   x1=x[1]
   xn=x[length(x)]
   b1=log(log(y1)/log(yn))/(x1^th-xn^th)
   b0=log(y1)*exp(-b1*x1^th)
   exp(b0*exp(b1*x^th));
  }

  expFct2 <- function(x, y_1,y_n,c,d,th)
  {
   x1=x[1]
   xn=x[length(x)]
   y1=c*y_1
   yn=d*y_n
   b1=log(log(y1)/log(yn))/(x1^th-xn^th)
   b0=log(y1)*exp(-b1*x1^th)
   exp(b0*exp(b1*x^th));
  }

  expFct <- function(x, b0, b1,th)
  {
    exp(b0*exp(b1*x^th))
  }


  nls.m1<- nls2(Retention ~ expFct1(Area, y1, yn, th), data = cl,
     start = list(th = 1))

  th_start=coef(nls.m1)


  nls.m2<- nls2(Retention ~ expFct2(Area, y1, yn, c, d,  th), data = cl,
     start = list(c=1, d=1, th = th_start))

  x=cl$Area
  x1=x[1]
  xn=x[length(x)]
  th_start=coef(nls.m2)[3]
  cy1=coef(nls.m2)[1]*y1

dyn=coef(nls.m2)[2]*yn
  b1_start=log(log(cy1)/log(dyn))/(x1^th_start-xn^th_start)
  b0_start=log(cy1)*exp(-b1_start*x1^th_start)

  th_start
  b1_start
  b0_start
  nls.m<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl,
     start = list(b0 = b0_start, b1 = b1_start, th = th_start))

  nls.m1
  nls.m2
  nls.m

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

Pinglei Gao
In reply to this post by Pinglei Gao
Dear Dave
Thanks for your kindness help. I am sorry, I am on a filed survey these days. I did not check my email for a long time. Please forgive me for the misunderstanding brought to you. Your answer works well for the model: Retention = exp (b0*exp (b1*Area^th)). Could you find the starting values for the another model: Retention = (b0*Area^(th+1))^b with the same data for me. I also asked the same question on R-help list before.

Bests,

PG

-----邮件原件-----
发件人: dave fournier [mailto:[hidden email]]
发送时间: 2016年10月20日 5:47
收件人: [hidden email]; [hidden email]
主题: Re: [R] Finding starting values for the parameters using nls() or nls2()

Actually this converges very nicely if you use these starting values that I obtained with AD Model Builder

        th     9.1180e-01
        b0    5.2104e+00
        b1   -4.6725e-04

The R result looks like

nls.m2
Nonlinear regression model
   model: Retention ~ expFct(Area, b0, b1, th)
    data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44, 15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class = "data.frame")
         b0         b1         th
  5.2104466 -0.0004672  0.9118029
  residual sum-of-squares: 686.8

Number of iterations to convergence: 1
Achieved convergence tolerance: 1.75e-06

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

dave fournier-2


  Unfortunately this problem does not appear to be well posed.

     Retention = (b0*Area^(th+1))^b

If b0, th, and b are the parameter only the product (th+1)*b is determined.

This comes from noting that powers satisfy



     (a^b)^c  =  a^(b*c)



So your model can be written as



           (b0^b) * Area^((th+1)*b)

______________________________________________
[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: Finding starting values for the parameters using nls() or nls2()

David Winsemius

> On Oct 25, 2016, at 9:29 AM, dave fournier <[hidden email]> wrote:
>
>
>
> Unfortunately this problem does not appear to be well posed.
>
>    Retention = (b0*Area^(th+1))^b
>
> If b0, th, and b are the parameter only the product (th+1)*b is determined.
>
> This comes from noting that powers satisfy
>
>
>
>    (a^b)^c  =  a^(b*c)
>
>
>
> So your model can be written as
>
>
>
>          (b0^b) * Area^((th+1)*b)
>

... which nicely completes the thread since one model had:

th1   =  9.1180e-01
 b01=    5.2104e+00
 b11 =  -4.6725e-04
(th1+1)*b11
[1] -0.0008932885


 b0  = 5.2104466   ;    b1 =   -0.0004672   ;  th =  0.9118029
((th+1)*b1)
[1] -0.0008931943

So both the R's nls2 and AD_MOdel_Builder results yield that same predictions for any given data point at least up to four decimal places.

So perhaps there was some other physical interpretation of those parameters and there exists an underlying theory that would support adding some extra constraints?

--
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: Finding starting values for the parameters using nls() or nls2()

dave fournier-2
 >
 >> On Oct 25, 2016, at 9:29 AM, dave fournier <davef at otter-rsch.com>
wrote:
 >>
 >>
 >>
 >> Unfortunately this problem does not appear to be well posed.
 >>
 >>    Retention = (b0*Area^(th+1))^b
 >>
 >> If b0, th, and b are the parameter only the product (th+1)*b is
determined.
 >>
 >> This comes from noting that powers satisfy
 >>
 >>
 >>
 >>    (a^b)^c  =  a^(b*c)
 >>
 >>
 >>
 >> So your model can be written as
 >>
 >>
 >>
 >>          (b0^b) * Area^((th+1)*b)
 >>
 >
 >... which nicely completes the thread since one model had:
 >
 >th1   =  9.1180e-01
 > b01=    5.2104e+00
 > b11 =  -4.6725e-04
 >(th1+1)*b11
 >[1] -0.0008932885
 >
 >
 > b0  = 5.2104466   ;    b1 =   -0.0004672   ;  th =  0.9118029
 >((th+1)*b1)
 >[1] -0.0008931943
 >
 >So both the R's nls2 and AD_MOdel_Builder results yield that same
predictions for any given data point at least up to four decimal places.
 >
 >So perhaps there was some other physical interpretation of those
parameters and there exists an underlying theory that would support
adding some extra constraints?
 >
 >--
 >>David Winsemius
 >Alameda, CA, USA

I'm not really sure what your point is. The OP has two models. One is well
posed and the other is not. I was discussing solution of the former model
which is well posed.  The form of the model, using a, b, and t and x,y to
simplify the expression is

         y  = exp( a * exp(b * x^t) )

My point is that there are many models like this where the obvious
parameterization (something like the parameters the user is interested
in) leads to parameter estimates with highly correlated errors.
This does not necessarily mean that
the model is badly determined. The model exists independent of the
particular parameterization. To fix  the ideas assume there are n
observations
(x_i,y_i)  and simplify by assuming that x_1<=x_2<=...<=x_n
(but not all equal lets hope)

A stable parameterization of the model can often be obtained by
picking as new parameters y1 and yn where y1 and yn are the predicted
values for y_1 and y_n, that is

            y1  = exp( a * exp(b * x_1^t) )
            yn  = exp( a * exp(b * x_n^t) )

Then one solves for some of the original model parameters in terms of y1
and yn.
The particular way this is done depends on the model. Often some has some
linear parameters and the procedure is easy. For this model as I showed
one can solve for a and b in terms of y1 and yn.
Then one can fit the model easily with AD Model Builder or
nls2 using this parameterization. nls2 provides the standard errors for
the parameter estimates.  The AD Model Builder solution provides the
estimated variance covariance matrix of the parameter estimates via the
standard maximum likelihood Hessian calculations. It also provides the
covariance for any desired dependent variable.  So one can fit the model
using y1,yn, and t and get the covariance matrix for a,b, and t
in one step. In nls2 one needs to fit the model using y1,yn and then
solve for a and b and run the model again from that point.  That is no big
deal, and I showed how to do it, but it is one more step for the user.

It is interesting to see the correlation matrix for the parameter estimates
and the dependent variables.
                     std dev          correlation
1  logt -9.233e-02 3.4026e-01  1.0000
2  c    9.7164e-01 3.7006e-02 -0.2690  1.0000
3  d    1.1010e+00 1.7852e-01 -0.7495  0.0325  1.000
4  t    9.1180e-01 3.1025e-01  1.0000 -0.2690 -0.749  1.0000
5  a    5.2104e+00 4.3404e-01 -0.9781  0.4191  0.621 -0.9781  1.000
6  b   -4.6725e-04 1.1714e-03  0.9994 -0.2836 -0.726  0.9994 -0.984 1.00

Here the independent variable are c, d, and logt
where y1=c*y_1  y2=d*y2
That is just an additional thing so that one can start with c=1 and d=1
Also logt is used where t=exp(logt) which makes sure t>0.
Notice that the correlation between c and d is 0.0325 which is small.


If  a and b were the parameters of the model

4   t    9.1180e-01 3.1025e-01   1.0000
5   a    5.2104e+00 4.3404e-01   -0.9781  1.000
6   b   -4.6725e-04 1.1714e-03    0.9994 -0.984  1.00

One can see how close to 1 and -1 the correlations are.

Notice that the parameter b is very badly determined.
So rather than saying the model is no good one sees that
the model is no good if one want to get information about
the parameter b. In contrast the parameter a is fairly well
determined and the parameter t is "kind of" determined.

   Because of this parameter confounding nls2 is extremely sensitive
   to the starting parameter values using this parameterization.
   If I change the parameters by about 15% or less as below

   #b0start=5.210445
   #b1_start=-0.0004672373
   #th_start=0.911804
   b0_start=5.7
   b1_start=-0.00055
   th_start=1.05

I get the dreaded

Error in (function (formula, data = parent.frame(), start, control =
nls.control(),  :
   singular gradient

So there is nothing wrong with either the data or the model. The model
just needs to be given a stable parameterization.

______________________________________________
[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.