Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Does the function y = at^b * exp(-ct) have a name? Can it be linearized? How can I estimate a, b, c?

Tags:

r

nlme

I am trying to fit a non-linear model, but can not find any good examples online.

Does this function have a name?

Can it be linearized?

I've attempted to estimate the parameters a, b, and c with a random effect g (as in group) as a function of time t, below. I can fit the model using nls without a random effect, but am having trouble getting the model to converge. Suggestions welcome (preferably within R, but any suitable package will do)?

## time, repeated 16 times for 4 replicates from each of 4 groups  
t <- rep(1:20, 16)
## g, group
g <- rep(1:4, each = 80)

## starting to create an example dataset, 
## to see if I can recover known parameters
a <- rep(c(3.5, 4, 4.1, 5), each = 80)
b <- rep(c(1.1, 1.4, 1.8, 2.5), each = 80)
c <- rep(c(0.125, 0.25), each = 160)

## error to add to above parameters
set.seed(1)
e_a <- runif(320, -0.5, 0.5)
e_b <- runif(320, -0.1, -0.1)
e_c <- runif(320, -0.02, 0.02)

## this is my function

f <- function(t, a, b, c) a * (t^b) * exp(-c * t)

## simulate y
y <- f(t = t, a + e_a, b + e_b, c + e_c)

mydata <- data.frame(t = t, y = y, g = g)

library(nlme)
## now fit the model to estimate a, b, c
fm1 <- nlme(y ~ a * (t^b) * exp(-c * t),
            data = mydata,
            fixed = a + b + c~1,
            random = a + b + c ~ 1|g,
            start = c(a = 4, b = 1, c = 0.25),
            method = "REML")
like image 482
Abe Avatar asked Oct 15 '15 20:10

Abe


1 Answers

In physics (and some other areas) I've seen this or variants of it called a Hoerl curve or Hoerl function e.g. here, though it has other names. If c is negative and a and b are positive it's a scaled gamma density.

When you ask about linearizing it, you have to be careful; the equation y = at^b . exp(ct) is not actually what you mean - the observations, y(i), are not exactly equal to a . t(i)^b . exp(ct(i)) (otherwise almost any 3 observations would give you the exact parameter values).

So the noise has to enter your model for y somehow. Is it additive? multiplicative, or something else? (Also important, but for other reasons: does its size change in some way as t changes, or not? Are the noise terms for different observations independent?)

If your actual model is y(i) = at(i)^b . exp(ct(i))+ε(i), that's not linearizable.

If your actual model is y(i) = at(i)^b . exp(ct(i)) . ε(i), and ε(i)=exp(η(i)) for some (hopefully zero-mean) η(i), that is linearizable.

Taking the second form,

log(y(i)) = log(a) + b log(t(i)) + c t(i) + log(ε(i))

or

y*(i) = a* + b.log(t(i)) + c.t(i) + η(i)

which is linear in the parameters a* = log(a), b and c, and the error term η(i); so if you're prepared to make that sort of an assumption about the error you should be able to fit it with methods suitable for such linear models; you may wish in that case to ponder the parenthetical questions about the error term above which may affect how you model it.

like image 171
Glen_b Avatar answered Nov 05 '22 03:11

Glen_b