I'm drawing a fitted normal distribution over a histogram. The use case is
trivial (fitting normal distributions on densities) but I want to extend it to
other fitting scenarios. What has stumped me so far is how to take the list
that is returned by nls() and use it for curve(). I realize that I can easily
do all of this with a few intermediate steps for any specific case. But I had
expected that it should be possible to get a parametrized(!) function that
computes predictions as one of the returned objects.
I.e. I want a function that works with the model of nls() like abline() works()
for lm().
I know that I can just pass parameters from coef(fit). But in the general case,
I don't know how many parameters there are. I thought since nls() is able to
put together the parametrized function internally, that function would be
passed into the results object. But that doesn't seem to be the case. And
though I could hack this together with parsing fit$m$formula to get() the
formula from the environment and then paste() the parameter list in there -
that sounds really, really awkward.
So I hope that there's an obvious way that I have overlooked.
Here's sample code to illustrate. :
fitNorm <- function(x, y) {
# fit a normal distribution
# Param: x domain
# y densities
# Value: the fit object
F <- function(x, a, mu, sig) { # some parametrized function
return( ( a / (sig*sqrt(2*pi)) ) * exp( (-1/2)*((x-mu)/sig)^2 ) )
}
mu <- weighted.mean(x, y) # estimate starting values
sig <- sd(sample(x, 1000, prob = y, replace = TRUE))
a <- 1
fit <- nls(y ~ F(x, a, mu, sig),
start = c(a = a, mu = mu, sig = sig)) # starting values
return(fit)
}
# Fit and plot ...:
# Values
x <- c(rnorm(5000, 3, 5), rnorm(2000, -5, 7)) # Two normal distributions ...
# Histogram
h <- hist(x, freq = FALSE,
breaks = seq(min(x)-2, max(x)+2, by = 1),
col = "#cfd7fa",
main = "", ylab = "density", xlab = "x")
# Fit
myFit <- fitNorm(h$mids, h$density)
# Now: what I can do is, patch together the model function ...
mF <- function(x,
a = coef(myFit)["a"],
mu = coef(myFit)["mu"],
sig = coef(myFit)["sig"]){
a / (sig*sqrt(2*pi)) * exp( (-1/2)*((x-mu)/sig)^2 )
}
# ... and add the curve:
curve(mF(x),
from = par("usr")[1], to = par("usr")[2],
col = "#FF000055", lwd = 2, add = TRUE)
# But what I would like to do is something much more general, like:
curve(myFit$func, from = myFit$from, to = myFit$to,
col = "#FF000055", lwd = 2, add = TRUE))
Thanks for any ideas and strategies.
Cheers,
Boris
--
Boris Steipe MD, PhD
University of Toronto
Associate Professor, Department of Biochemistry and
Department of Molecular Genetics
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
