[R] nls question
Michael A. Miller
mmiller3 at iupui.edu
Fri Aug 23 16:42:58 CEST 2002
I'm using the nls library to fit peaks in spectra. The model
that I'm using is a quadratic background plus a Gaussian peak.
The peak is well described by the Gaussian and the background is
mostly linear with a bit of quadratic. My goal is an estimate of
the integral (and the uncertainty of the integral) of the peak
without the back ground.
The procedure that I'm using is as follows - first, smooth the
data with modreg's smooth.spline. Use this smoothed curve to
make an initial estimate of the peak position, height and width.
Next, use the smoothed data below the peak to calculate initial
values for the quadratic background.
Finally, fit the spectrum like this:
back.ground.fnc <- function( E, a0, a1, a2 ) {
result <- a0 + a1*E + a2*E*E
return(result)
}
spectrum.fnc <- function(E, n0, x0, s, a0, a1, a2 ) {
back.ground <- back.ground.fnc( E, a0, a1, a2 )
## Make sure the background is non-negative:
back.ground <- back.ground * (back.ground > 0)
result <- n0 * exp( - ( ( E - x0 ) / s )^2 ) + back.ground
return(result)
}
##===================================================================================
weighted.spectrum.fnc <- function( counts, E, n0, x0, s, a0, a1, a2 ) {
epsilon <- 1e-20
pred <- spectrum.fnc( E, n0, x0, s, a0, a1, a2 )
w <- sqrt(abs(counts+epsilon))
(counts - pred) / w
}
spectrum.fit <- nls( ~ weighted.spectrum.fnc( counts, E, n0, x0, s, a0, a1, a2 ),
data = fit.df,
start = list( n0=n0, x0=x0, s=s, a0=a0, a1=a1, a2=a2 )
)
I've found that, for data sets that are similar sometimes the fit
works well and sometimes it fails when it hits step factor or
iteration limits. In trying to sort out what's going on, I've
managed to confuse myself pretty badly. Can someone help me to
understand how I can make nls give me a successful fit? Is there
an nls library tutorial kicking around anywhere?
Mike
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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