R help - Dec 2001 - fit to spike with exponential decay : optim() question

I finally got (mostly) what I wanted.  In an attempt to figure out how to
get nls to deal with a non-differentiable function, I had (stupidly)
'simplified' the problem until it became singular.

Can I do something to make optim() less sensitive to my initial guess? For
this example, I get a lousy solution if I make the initial guess for t0 min(t) =
0.05.

Thanks again,
--
Robert Merithew
LASSP, Clark Hall
Cornell University, Ithaca NY


-----------

t <- c(0.05,0.9,1.4,2.38,3.42,5.4,8.31,12.4)
amp <- c(1.0,0.85,7.4, 6.1, 4.95, 3.5, 2.3, 1.5)

spike <- function (x, t) {
  b0 <- x[1]
  b1 <- x[2]
  tau <- x[3]
  t0 <- x[4]
  
  temp <- exp((-t+t0)/tau)
  (b0 + (b1 * temp)  * (t > t0))
}

spike.sos <- function (x) {
  sum((amp - spike(x, t))^2)
}

guess <- c(min(amp), max(amp)-min(amp), (max(t)-min(t))/3, 0)

decay.opt <- optim(guess, spike.sos, control=list(trace=T))

xr <- (0:140)/10
plot (xr, spike(decay.opt$par, xr), type="l", col="blue")
points (t, amp)

--------


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

On Fri, 14 Dec 2001, Robert D. Merithew wrote:
>
> I finally got (mostly) what I wanted.  In an attempt to figure out how to
> get nls to deal with a non-differentiable function, I had (stupidly)
> 'simplified' the problem until it became singular.
>
> Can I do something to make optim() less sensitive to my initial guess? For
> this example, I get a lousy solution if I make the initial guess for t0
> min(t) = 0.05.
1) Give it derivatives
2) Use a different method, such as BFGS.
3) Scale the problem as described in help(optim): not that yours seems
   badly scaled, but I would optimize the mean square in larger problems.

Recovering from lousy information is not what optimizers are designed for,
though.
>
> Thanks again,
> --
> Robert Merithew
> LASSP, Clark Hall
> Cornell University, Ithaca NY
>
>
> -----------
>
> t <- c(0.05,0.9,1.4,2.38,3.42,5.4,8.31,12.4)
> amp <- c(1.0,0.85,7.4, 6.1, 4.95, 3.5, 2.3, 1.5)
>
> spike <- function (x, t) {
>   b0 <- x[1]
>   b1 <- x[2]
>   tau <- x[3]
>   t0 <- x[4]
>
>   temp <- exp((-t+t0)/tau)
>   (b0 + (b1 * temp)  * (t > t0))
> }
>
> spike.sos <- function (x) {
>   sum((amp - spike(x, t))^2)
> }
>
> guess <- c(min(amp), max(amp)-min(amp), (max(t)-min(t))/3, 0)
>
> decay.opt <- optim(guess, spike.sos, control=list(trace=T))
>
> xr <- (0:140)/10
> plot (xr, spike(decay.opt$par, xr), type="l",
col="blue")
> points (t, amp)
>
> --------
>
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at
stat.math.ethz.ch
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._