R help - Oct 2001 - nls() fit to a lorentzian - can I specify partials?

First, thanks to all who helped me with my question about rescaling axes
on the fly.  Using unlist() and range() to set the axis ranges in advance
worked well.  I've since plotted about 300 datasets with relative ease.

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).

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

Can I specify the partials in advance?

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

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


thanks,

--
Robert Merithew
LASSP, Cornell University

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._