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

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Boris Steipe
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: ??? is to nls() as abline() is to lm() ?

Duncan Murdoch-2
I haven't followed your example closely, but can't you use the predict()
method for this?  To draw a curve, the function that will be used in
curve() sets up a newdata dataframe and passes it to predict(fit,
newdata= ...) to get predictions at those locations.

Duncan Murdoch

On 17/10/2020 5:27 a.m., Boris Steipe wrote:

> 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-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: ??? is to nls() as abline() is to lm() ?

Boris Steipe
I had tried predict() but overlooked that it is picky about needing the data in a list-like object. When I passed a vector it just returned the predictions on the original values. Passing a list (as in the help example... of course) does the trick.

This works - and does exactly what I had hoped for:

curve(predict(myFit, list(x = x)), col = "#FF000055", lwd = 2, add = TRUE)


Very grateful!
Boris







> On 2020-10-17, at 20:26, Duncan Murdoch <[hidden email]> wrote:
>
>
> I haven't followed your example closely, but can't you use the predict()
> method for this?  To draw a curve, the function that will be used in
> curve() sets up a newdata dataframe and passes it to predict(fit,
> newdata= ...) to get predictions at those locations.
>
> Duncan Murdoch
>
> On 17/10/2020 5:27 a.m., Boris Steipe wrote:
>> 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-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>

______________________________________________
[hidden email] 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.