[R] nls() fit to a lorentzian - can I specify partials?

[Douglas Bates]

"Robert D. Merithew" <merithew at ccmr.cornell.edu> writes:

> Now I'm trying to fit a lossy oscillator resonance to (the square root of)
> a lorentzian (testframe$y is oscillator amplitude, testframe$x is drive
> frequency):
> 
> lorentz <- nls( y ~ a0 / (Q*sqrt((1-(x/f0)^2)^2 + ((x/f0) * 1/Q)^2)),
>                data = testframe,
>                start = list (a0=1.4, f0=84321.6, Q=600000))
> 
> and am running into lots of convergence trouble (singular derivatives,
> step factors heading off to 1/infinity, etc).

The documentation for the nls function is not as complete as it should
be.  

> Does nls take symbolic partial derivatives of my expression?  (I suspect
> it's doing derivatives numerically) 

It uses numerical derivatives.

> Can I specify the partials in advance?

Yes.  Define a function Lorentzian with arguments x, a0, f0, and Q.
It should return the function value as a vector of the same length as
x with an attribute named "gradient" that is an length(x) by 3
matrix.  The column names of the matrix should be "a0", "f0", and "Q"

> I noticed the 'deriv' function, but am baffled by its output:
> 
> > deriv(~ a0 / (Q*sqrt((1-(x/f0)^2)^2 + ((x/f0) * 1/Q)^2)),"a0")
> expression({
>     .expr1 <- x/f0
>     .expr10 <- Q * sqrt(1 - .expr1^2^2 + .expr1 * 1/Q^2)
>     .value <- a0/.expr10
>     .grad <- array(0, c(length(.value), 1), list(NULL, c("a0")))
>     .grad[, "a0"] <- 1/.expr10
>     attr(.value, "gradient") <- .grad
>     .value
> })
> 
> 
> I think .expr10 should look more like:
> 
> 	Q * sqrt((1 - .expr1^2)^2 + (.expr1 * 1/Q)^2)
> 
> Is there simply a problem with the way the output of deriv() is displayed
> (missing 2 sets of parentheses)?

I imagine the actual expression tree for .expr10 corresponds to the
expression you wrote but it is being deparsed without the parentheses.

> Finally:  Is there a better way to do my lorentzian fits with R?

The parameter a0 is a conditionally linear parameter.  You can remove
that from the optimization by using the "plinear" algorithm.  I would
also recommend setting trace = TRUE so you can see the steps that are
being taken in the algorithm.

lorentz <- nls( y ~ 1 / (Q*sqrt((1-(x/f0)^2)^2 + ((x/f0) * 1/Q)^2)),
               data = testframe, alg = "plinear",
               start = list (f0=84321.6, Q=600000), trace = TRUE)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._