Starting estimates for nls Exponential Fit

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Starting estimates for nls Exponential Fit

Anto
Hello everyone,

I have come across a bit of an odd problem: I am currently analysing data PCR reaction data of which part is behaving exponential. I would like to fit the exponential part to the following:

y0 + alpha * E^t

In which Y0 is the groundphase, alpha a fluorescence factor, E the efficiency of the reaction & t is time (in cycles)

I can get this to work for most of my reactions, but part fails the nls fit (Convergence failure: iteration limit reached without convergence). This mainly has to do with the starting estimates for the nls function. I have used various 'smart' ways of calculating starting estimates (grid searches, optim(), etc.) but however close the starting estimates are to the actual values, nls still fails for many datasets.

The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and E=1.85) which are totaly arbitray and these DO work for, say 99% of the datasets. I don't know why this is and I don't why my 'good estimates' do not work. I am desperatly seeking a way of calculating working starting estimates for my nls function. Can anybody give me a hand?

thanks,

Anto

R version 2.9.2

example dataset:

ExponValues
[1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47

ExponCycles
[1] 17 18 19 20 21 22 23 24 25


Example starting estimate calculation

Y0 <- ExponValues[1]
E <- weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3)
alpha <-  weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3)


Optional parameter optimisation:

Esp <- optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles))))})$par


nls function:

Emodel<-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp, algorithm="port"),silent=TRUE)
Reply | Threaded
Open this post in threaded view
|

Re: Starting estimates for nls Exponential Fit

Katharine Mullen
If you could reformulate your model as alpha * (y0 + E^t) then you could
use nls with alg="plinear" (alpha then would be eliminated from the
nonlinear param and treated as conditionally linear) and this would help
with convergence.

Else you can try package DEoptim to get the starting values; the advantage
over optim is that you need then lower and upper bounds on the parameters,
not more starting values; the bounds however should be appropriate and as
tight as possible.

Also: by default nls uses max. 50 iter.  Depending on where you start off
you may need more than this.

##

ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
                 2468.32,2778.47)

ExponCycles <- c(17,18,19,20,21,22,23,24,25)

library(DEoptim)

mod <- function(x) x[1] + x[2]*x[3]^ExponCycles
fun <- function(x) sum((ExponValues-mod(x))^2)

ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
              control=list(trace=FALSE))

pa <- ss$optim$bestmem

## now have par est that give OK fit
matplot(cbind(ExponValues, mod(pa)),type="l")

names(pa) <- c("Y0", "a", "E")

## fit w/port and GN
Emodel<- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)),
             start=pa, algorithm="port",
             control=list(maxiter=1000, warnOnly=TRUE))

Emodel1 <- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)), start=pa,
             control=list(maxiter=1000, warnOnly=TRUE))

## fit
matplot(cbind(ExponValues, fitted(Emodel), fitted(Emodel1)),type="l")

On Tue, 1 Dec 2009, Anto wrote:

>
> Hello everyone,
>
> I have come across a bit of an odd problem: I am currently analysing data
> PCR reaction data of which part is behaving exponential. I would like to fit
> the exponential part to the following:
>
> y0 + alpha * E^t
>
> In which Y0 is the groundphase, alpha a fluorescence factor, E the
> efficiency of the reaction & t is time (in cycles)
>
> I can get this to work for most of my reactions, but part fails the nls fit
> (Convergence failure: iteration limit reached without convergence). This
> mainly has to do with the starting estimates for the nls function. I have
> used various 'smart' ways of calculating starting estimates (grid searches,
> optim(), etc.) but however close the starting estimates are to the actual
> values, nls still fails for many datasets.
>
> The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and
> E=1.85) which are totaly arbitray and these DO work for, say 99% of the
> datasets. I don't know why this is and I don't why my 'good estimates' do
> not work. I am desperatly seeking a way of calculating working starting
> estimates for my nls function. Can anybody give me a hand?
>
> thanks,
>
> Anto
>
> R version 2.9.2
>
> example dataset:
>
> ExponValues
> [1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47
>
> ExponCycles
> [1] 17 18 19 20 21 22 23 24 25
>
>
> Example starting estimate calculation
>
> Y0 <- ExponValues[1]
> E <-
> weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3)
> alpha <-
> weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3)
>
>
> Optional parameter optimisation:
>
> Esp <-
> optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles))))})$par
>
>
> nls function:
>
> Emodel<-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp,
> algorithm="port"),silent=TRUE)
> --
> View this message in context: http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p932230.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> 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
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: Starting estimates for nls Exponential Fit

Prof J C Nash (U30A)
In reply to this post by Anto
Kate Mullen showed one approach to this problem by using DEOptim to get
some good starting values.

However, I believe that the real issue is scaling (Warning: well-ridden
hobby-horse!).

With appropriate scaling, as below, nls does fine. This isn't saying nls
is perfect -- I've been trying to figure out how to do a nice job of
helping folk to scale their problems. Ultimately, it would be nice to
has an nls version that will do the scaling and also watch for some
other situations that give trouble.

Cheers, JN


## JN test
rm(list=ls())

ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
                 2468.32,2778.47)

ExponCycles <- c(17,18,19,20,21,22,23,24,25)

mod <- function(x) x[1] + x[2]*x[3]^ExponCycles

modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles)

fun <- function(x) sum((ExponValues-mod(x))^2)



pa<-c(1,2,3)
lo<-c(0,0,0)
up<-c(20,20,20)
names(pa) <- c("Y0", "a", "E")

## fit w/port and GN
Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),
             start=pa, algorithm='port', lower=lo, upper=up,
             control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))

Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,
             control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))

## fit
matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l")

______________________________________________
[hidden email] mailing list
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: Starting estimates for nls Exponential Fit

Katharine Mullen
You used starting values:
   pa <- c(1,2,3)

but both algorithms (port and Gauss-Newton) fail if you use the slightly
different values:
   pa <- c(1,2,3.5)

Scaling does not fix the underlying sensitivity to starting values.
pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both
alg. also fail if there is insufficient spread (e.g., both fail for
pa <- c(1,1.25,1.5) ).

For the record, DE uses (by default at least) a random start and for bad
starts will sometimes fail to give good results; decrease the probability
this happens by increasing itermax from the default itermax=200, as in:

ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
              control=list(trace=FALSE, itermax=1000))

On Wed, 2 Dec 2009, Prof. John C Nash wrote:

> Kate Mullen showed one approach to this problem by using DEOptim to get
> some good starting values.
>
> However, I believe that the real issue is scaling (Warning: well-ridden
> hobby-horse!).
>
> With appropriate scaling, as below, nls does fine. This isn't saying nls
> is perfect -- I've been trying to figure out how to do a nice job of
> helping folk to scale their problems. Ultimately, it would be nice to
> has an nls version that will do the scaling and also watch for some
> other situations that give trouble.
>
> Cheers, JN
>
>
> ## JN test
> rm(list=ls())
>
> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
>                  2468.32,2778.47)
>
> ExponCycles <- c(17,18,19,20,21,22,23,24,25)
>
> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles
>
> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles)
>
> fun <- function(x) sum((ExponValues-mod(x))^2)
>
>
>
> pa<-c(1,2,3)
> lo<-c(0,0,0)
> up<-c(20,20,20)
> names(pa) <- c("Y0", "a", "E")
>
> ## fit w/port and GN
> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),
>              start=pa, algorithm='port', lower=lo, upper=up,
>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>
> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,
>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>
> ## fit
> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l")
>
> ______________________________________________
> [hidden email] mailing list
> 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
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: Starting estimates for nls Exponential Fit

Prof J C Nash (U30A)
Kate is correct. The parameter scaling helps quite a bit, but not enough
to render the problem "nice" so that many "reasonable" starting points
will give useful results. Indeed, a run using "all.methods=TRUE" in our
optimx package (on r-forge at
http://r-forge.r-project.org/R/?group_id=395) gives

                          par  fvalues      method   fns grs itns conv
4   2.194931, 1.000001, 1.000027 566407.6     nlm    NA  NA   3    0
1 2.1949335, -9.0413113, 0.7516435 566407.6   BFGS    37   6 NULL    0
2  2.1950009, 0.2548318, 1.1163498 566404.6 Nelder    82  NA NULL    0
3     1.974226, 1.829957, 1.681338 2409.771   SANN 10000  NA NULL    0
6  1.9904045, 0.6151421, 1.7559401 1696.497  ucminf    51  51 NULL    0
5  1.9906488, 0.6012996, 1.7575365 1696.248  nlminb    80 166   51    0
   KKT1  KKT2
4 FALSE FALSE
1  TRUE FALSE
2 FALSE FALSE
3 FALSE  TRUE
6 FALSE  TRUE
5 FALSE  TRUE
>

A bit of a dog's dinner.

On the positive side, this is a simple but nasty problem to illustrate
lots of the issues with nonlinear parameter estimation.

JN



Katharine Mullen wrote:

> You used starting values:
>    pa <- c(1,2,3)
>
> but both algorithms (port and Gauss-Newton) fail if you use the slightly
> different values:
>    pa <- c(1,2,3.5)
>
> Scaling does not fix the underlying sensitivity to starting values.
> pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both
> alg. also fail if there is insufficient spread (e.g., both fail for
> pa <- c(1,1.25,1.5) ).
>
> For the record, DE uses (by default at least) a random start and for bad
> starts will sometimes fail to give good results; decrease the probability
> this happens by increasing itermax from the default itermax=200, as in:
>
> ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
>               control=list(trace=FALSE, itermax=1000))
>
> On Wed, 2 Dec 2009, Prof. John C Nash wrote:
>
>> Kate Mullen showed one approach to this problem by using DEOptim to get
>> some good starting values.
>>
>> However, I believe that the real issue is scaling (Warning: well-ridden
>> hobby-horse!).
>>
>> With appropriate scaling, as below, nls does fine. This isn't saying nls
>> is perfect -- I've been trying to figure out how to do a nice job of
>> helping folk to scale their problems. Ultimately, it would be nice to
>> has an nls version that will do the scaling and also watch for some
>> other situations that give trouble.
>>
>> Cheers, JN
>>
>>
>> ## JN test
>> rm(list=ls())
>>
>> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
>>                  2468.32,2778.47)
>>
>> ExponCycles <- c(17,18,19,20,21,22,23,24,25)
>>
>> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles
>>
>> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles)
>>
>> fun <- function(x) sum((ExponValues-mod(x))^2)
>>
>>
>>
>> pa<-c(1,2,3)
>> lo<-c(0,0,0)
>> up<-c(20,20,20)
>> names(pa) <- c("Y0", "a", "E")
>>
>> ## fit w/port and GN
>> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),
>>              start=pa, algorithm='port', lower=lo, upper=up,
>>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>>
>> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,
>>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>>
>> ## fit
>> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l")
>>
>> ______________________________________________
>> [hidden email] mailing list
>> 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
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: Starting estimates for nls Exponential Fit

dave fournier
Figuring out the best parameterization for this kind of model
is a bit tricky until you get the hang of it.
Let the function be

        y_t = y_0 + alpha * E^t

where uppercase Y_t denotes an observed value
and lower case y_t is a predicted value.
Index the times by t_1 .... t_n
WLOG assume that t_1 is the smallest time and t_n is the largest time.

Now we already have a good idea what the predicted values  y_1 and y_n
should be as we have observations for them. We have the two equations

       y_1 = y_0 + alpha * E^t_1

       y_n = y_0 + alpha * E^t_n

we can solve these for alpha in terms of y_1,y_n,and E giving

     alpha = (y_n-y_1)/(E^t_n -E^t_1)       (1)

and solve for y_0  as

         y_0 = y_1 - alpha * E^t_1  using (1) for alpha

Now use the good estimates Y_1 and Y_n  as the starting values
for y_1 and y_n
and try some "reasonable value for E (say 0.1 < E < 100)

Do the minimization in two stages first holding y_1 and y_n fixed and
then estimate y_1,y_n,and E together. This converges in less than a second
using AD Model Builder and the starting values (large value for E.

2018.34 2778.47 exp(10.0)  where I deliberately
used exp(10) as a large initial value for E.
The parameter estimates together with the est std devs
are

  1   y_1     1.9994e+03 3.9177e-01
  2   y_n     2.7881e+03 6.7557e-01
  3   log_E  5.6391e-01 1.2876e-03
  4   alpha      6.0130e-04 1.9398e-05
  5   y_0     1.9906e+03 4.5907e-01
  6   E      1.7575e+00 4.3935e-02

There are problems for E near 1 which need to be dealt with if
you have to go there, but that is just a technicality.

This idea also works well for a logistic say 3 4 or 5 parameter.














--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

______________________________________________
[hidden email] mailing list
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: Starting estimates for nls Exponential Fit

Anto
In reply to this post by Prof J C Nash (U30A)
Thanks everybody,

This has been quite helpful, the problem remains tricky but at least now I've got a version of my script that handles all my reactions without error. The DEoptim solution produced good starting values for a lot of reactions, but sadly not for all. I now use scaled parameters and allow more iterations than the standard 50 (1000). I doubt it is completely stable and I think I will run into reactions that will fail the fit, but for now everything works fine so I won't be complaining.

Thanks a lot,

Antoon

 


Prof. John C Nash wrote
Kate is correct. The parameter scaling helps quite a bit, but not enough
to render the problem "nice" so that many "reasonable" starting points
will give useful results. Indeed, a run using "all.methods=TRUE" in our
optimx package (on r-forge at
http://r-forge.r-project.org/R/?group_id=395) gives

                          par  fvalues      method   fns grs itns conv
4   2.194931, 1.000001, 1.000027 566407.6     nlm    NA  NA   3    0
1 2.1949335, -9.0413113, 0.7516435 566407.6   BFGS    37   6 NULL    0
2  2.1950009, 0.2548318, 1.1163498 566404.6 Nelder    82  NA NULL    0
3     1.974226, 1.829957, 1.681338 2409.771   SANN 10000  NA NULL    0
6  1.9904045, 0.6151421, 1.7559401 1696.497  ucminf    51  51 NULL    0
5  1.9906488, 0.6012996, 1.7575365 1696.248  nlminb    80 166   51    0
   KKT1  KKT2
4 FALSE FALSE
1  TRUE FALSE
2 FALSE FALSE
3 FALSE  TRUE
6 FALSE  TRUE
5 FALSE  TRUE
>

A bit of a dog's dinner.

On the positive side, this is a simple but nasty problem to illustrate
lots of the issues with nonlinear parameter estimation.

JN



Katharine Mullen wrote:
> You used starting values:
>    pa <- c(1,2,3)
>
> but both algorithms (port and Gauss-Newton) fail if you use the slightly
> different values:
>    pa <- c(1,2,3.5)
>
> Scaling does not fix the underlying sensitivity to starting values.
> pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both
> alg. also fail if there is insufficient spread (e.g., both fail for
> pa <- c(1,1.25,1.5) ).
>
> For the record, DE uses (by default at least) a random start and for bad
> starts will sometimes fail to give good results; decrease the probability
> this happens by increasing itermax from the default itermax=200, as in:
>
> ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
>               control=list(trace=FALSE, itermax=1000))
>
> On Wed, 2 Dec 2009, Prof. John C Nash wrote:
>
>> Kate Mullen showed one approach to this problem by using DEOptim to get
>> some good starting values.
>>
>> However, I believe that the real issue is scaling (Warning: well-ridden
>> hobby-horse!).
>>
>> With appropriate scaling, as below, nls does fine. This isn't saying nls
>> is perfect -- I've been trying to figure out how to do a nice job of
>> helping folk to scale their problems. Ultimately, it would be nice to
>> has an nls version that will do the scaling and also watch for some
>> other situations that give trouble.
>>
>> Cheers, JN
>>
>>
>> ## JN test
>> rm(list=ls())
>>
>> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
>>                  2468.32,2778.47)
>>
>> ExponCycles <- c(17,18,19,20,21,22,23,24,25)
>>
>> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles
>>
>> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles)
>>
>> fun <- function(x) sum((ExponValues-mod(x))^2)
>>
>>
>>
>> pa<-c(1,2,3)
>> lo<-c(0,0,0)
>> up<-c(20,20,20)
>> names(pa) <- c("Y0", "a", "E")
>>
>> ## fit w/port and GN
>> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),
>>              start=pa, algorithm='port', lower=lo, upper=up,
>>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>>
>> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,
>>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>>
>> ## fit
>> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l")
>>
>> ______________________________________________
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>

______________________________________________
R-help@r-project.org mailing list
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: Starting estimates for nls Exponential Fit

dave fournier

I thought maybe my suggestion for reparameterizing this simple problem
was ignored because I didn't supply R code for the problem.
Here it is using optim for the optimization. It converges trivially with
an initial value for E of 1000.
As I stated before, there is nothing at all difficult about this
problem. You simply need to parameterize it properly. Of course that is
not to say that rescaling is not useful as well, but the important thing
is to parameterize the model properly.

ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9


#keep y1 and yn fixed and get initial value for E
Esp1 <- optim(c(E=E),method ="BFGS",
  function(x)
  {
    E=x[1]
    a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
    Y0=y1-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^ExponCycles)
    return(1000*sum(diff*diff))
  })$par

E=Esp1[1]

Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",
  function(x)
  {
    E=x[3]
    a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
    Y0=x[1]-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^ExponCycles)
    return(1000*sum(diff*diff))
  })$par

y1=Esp[1]
y2=Esp[2]
E=Esp[3]

a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];





--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

______________________________________________
[hidden email] mailing list
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: Starting estimates for nls Exponential Fit

dave fournier
Thanks to Dennis Murphy who pointed out that ExponCycles
is undefined. It is an R gotcha. I had shortened the name but R still
remembered it so the script worked but only on my computer.
This should fix that.

ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9


#keep y1 and yn fixed and get initial value for E
Esp1 <- optim(c(E=E),method ="BFGS",
  function(x)
  {
    E=x[1]
    a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
    Y0=y1-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

E=Esp1[1]

Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",
  function(x)
  {
    E=x[3]
    a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
    Y0=x[1]-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

y1=Esp[1]
y2=Esp[2]
E=Esp[3]

a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];




--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

______________________________________________
[hidden email] mailing list
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: Starting estimates for nls Exponential Fit

Anto
I have tried the method proposed by Dave, and I must say it works very well. Not to yield starting estimates for an nls-fit, but as an independent method for calculating E (which is, by the way, the only paramater that I am actually interested in).

The calculated values for E (Esp[3]) are on average 1.40e-06 lower than the values produced by the nls fit (calculated over 96 reactions). I am using the found E values to compenstate for differences between reaction efficiency in genetic quantification assays. An efficiency difference of a few %'s (+/- 0.05 in absolute value) can cause quantification differences of several dozen percentages. Since I am in the process of comparing different efficiency calculation methods (of which exopnential fit is one) I'll compare the nls-fit results with the results by Dave's method to see if they yield any significant differnces. I'll probably post the results by the end of next week (I have an alarming number of reports still due, for which I am recieving increasingly frequent angry looks by my coworkers. I had the revision of an article as an excuse for further procrastination, but now that it has been accepted I'll temporarly have to shift the focus of my work).

thanks a lot for all your time & effort

Antoon



dave fournier wrote
Thanks to Dennis Murphy who pointed out that ExponCycles
is undefined. It is an R gotcha. I had shortened the name but R still
remembered it so the script worked but only on my computer.
This should fix that.

ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9


#keep y1 and yn fixed and get initial value for E
Esp1 <- optim(c(E=E),method ="BFGS",
  function(x)
  {
    E=x[1]
    a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
    Y0=y1-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

E=Esp1[1]

Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",
  function(x)
  {
    E=x[3]
    a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
    Y0=x[1]-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

y1=Esp[1]
y2=Esp[2]
E=Esp[3]

a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];




--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

______________________________________________
R-help@r-project.org mailing list
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: Starting estimates for nls Exponential Fit

Christian Ritz-3
Hi Antoon,

now that you mention trying out different methods, maybe you should consider fitting a
sigmoidal curve to the entire dataset and not only the exponential part (which constitutes
a very small dataset) as seems to have been the endeavour that initiated the posting to
R-help.

One options would then be to use the package qpcR and simply fitting a sigmoidal model
using for instance drm() in the package drc or nls() combined with SSfpl().

Just a suggestion.


Christian

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