5

我正在尝试拟合非线性模型，但在网上找不到任何好的示例。

这个函数有名字吗？

可以线性化吗？

我试图估计参数a,b和c随机效应g（如在组中）作为 time 的函数t，如下所示。nls我可以在没有随机效应的情况下拟合模型，但无法让模型收敛。欢迎提出建议（最好在 R 中，但任何合适的包都可以）？

## 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")
4

1 回答 1

9

在物理学（和其他一些领域）中，我已经看到了这个或它的变体，称为 Hoerl 曲线或 Hoerl 函数，例如这里，尽管它有其他名称。如果c为负且a和b为正，则它是按比例缩放的伽马密度。

当您询问线性化它时，您必须小心；方程y = at ^ b。exp( ct ) 实际上并不是您的意思-观察结果y ( i ) 并不完全等于a。t (i)^ b。exp( ct ( i )) （否则几乎任何 3 个观察结果都会为您提供确切的参数值）。

所以噪音必须以某种方式进入你的模型。是添加剂吗？乘法还是别的？（同样重要，但出于其他原因：它的大小是否会随着t的变化而变化？不同观察的噪声项是否独立？）

如果您的实际模型是 y ( i ) = at ( i )^ b。exp( ct ( i ))+ε( i )，这不是线性化的。

如果您的实际模型是y ( i ) = at ( i )^ b。exp( ct ( i ))。ε(i)，对于一些（希望是零均值）η( i )，ε( i )=exp(η( i ) )，这是可线性化的。

采取第二种形式，

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

或者

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

这在参数a * = log( a )、b和c以及误差项 η( i ) 中是线性的；因此，如果您准备对错误做出这种假设，您应该能够使用适合此类线性模型的方法来拟合它；在这种情况下，您可能希望考虑有关上述错误术语的括号问题，这可能会影响您对其建模的方式。

于 2015-10-15T23:13:50.527 回答