[R] nls() fit to a lorentzian - can I specify partials?
Douglas Bates
bates at stat.wisc.edu
Fri Oct 5 19:18:22 CEST 2001
"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list