# ??? is to nls() as abline() is to lm() ?

3 messages
Open this post in threaded view
|

## ??? is to nls() as abline() is to lm() ?

 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 ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.