如何找到非线性模型的起始值?

How do I find the starting values for a nonlinear model?

在 Windows 10.

上使用 RStudio v1.3.1093、R v.4.0.2

我正在处理作业问题,并给出了一个数据集:

x <- cbind(c(8, 8, 14, 14, 14, 16, 16, 16, 18, 18, 20, 20, 20, 22, 22, 22, 24, 24, 24, 26, 26, 26, 28, 28, 30, 30, 30, 32, 32, 34, 36, 36, 42))
y <- cbind(c(0.49, 0.49, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45, 0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 0.40, 0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 0.41, 0.40, 0.40, 0.41, 0.38, 0.39))

df <- data.frame('x' = x, 'y' = y)

我们想要拟合模型 Mitcherlich Law 模型:y = a - b*exp(-c*x) + e 然后讨论我们如何获得起始值。

我用过:

i <- getInitial(y ~ SSasymp(x, a, b, c), data = df)

获取我的起始值,但是当我拟合模型时:

fit <- nls(y ~ a - b*exp(-c*x), data = df, start = list(a = i[1], b = i[2], c = i[3]))

我得到:

Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates

我是不是使用了错误的函数,或者是否有其他方法可以计算起始值?

我想你可能用错了自启动功能。在这种情况下,可以绘制点并使用参数旋转,同时绘制生成的曲线以足够接近 nls 工作:

x <- cbind(c(8, 8, 14, 14, 14, 16, 16, 16, 18, 18, 20, 20, 20,
             22, 22, 22, 24, 24, 24, 26, 26, 26, 28, 28, 30, 30, 
             30, 32, 32, 34, 36, 36, 42))

y <- cbind(c(0.49, 0.49, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 
             0.45, 0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 
             0.40, 0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 
             0.41, 0.40, 0.40, 0.41, 0.38, 0.39))

df <- data.frame('x' = x, 'y' = y)

i <- c(a = -0.5, b = -1, c = 0.1)
fit <- nls(y ~ a - b*exp(-c*x), data = df, start = as.list(i))

fit
#> Nonlinear regression model
#>   model: y ~ a - b * exp(-c * x)
#>    data: df
#>        a        b        c 
#>  0.38621 -0.21016  0.09033 
#>  residual sum-of-squares: 0.003971
#> 
#> Number of iterations to convergence: 4 
#> Achieved convergence tolerance: 2.089e-08

plot(df)
lines(5:50, predict(fit, newdata = list(x = 5:50)), col = "red", lty = 2)

如果你想要一些能够自动接近起点的东西(并且不想写一个自启动),你可以做一些假设:

  1. 假设 c 为正(即该图显示像半衰期曲线一样的衰减,如您的数据所示),那么 exp(-c * x) 将趋于零且 x 较大,因此只要您有在您的数据合理范围内,y 的最小值可能接近 a
  2. 如果我们从估计的 a 中减去 y 数据,b 的值与通过结果点进行线性回归的截距不会太远。
  3. 通过取这些新 y 值的对数除以 b 创建的回归线的斜率将接近 -c

所以我们可以为这种类型的曲线创建一个粗略的估计器,如下所示:

roughstart <- function(x, y) {
    xy <- data.frame(x = x, y = y)
    z <- xy[["y"]]
    a <- min(z)
    xy$z <- a - z
    b <- coef(lm(z ~ x, xy))[1]
    xy$z <- log(xy$z/b)
    c <- -coef(lm(z ~ x, xy[is.finite(xy$z),]))[2]
    parms <- as.numeric(c(a, b, c))
    setNames(as.list(parms), c("a", "b", "c"))
}

所以我们可以这样做:

nls(y ~ a - b*exp(-c*x), data = df, start = roughstart(df$x, df$y))
Nonlinear regression model
  model: y ~ a - b * exp(-c * x)
   data: df
       a        b        c 
 0.38621 -0.21016  0.09033 
 residual sum-of-squares: 0.003971

Number of iterations to convergence: 4 
Achieved convergence tolerance: 6.616e-06

reprex package (v0.3.0)

于 2020-12-05 创建

1) 运行 nls with SSasymp 然后适当地转换参数给出如下 st 中的答案。如果你想要一个干净的 运行 使用你的参数化 运行 nls 第二次使用转换后的参数作为起始值也如下所示。第二个 运行 将有 0 次迭代,因为我们从最优开始。

fit0 <- nls(y ~ SSasymp(x, Asym, R0, lrc), df)
st <- with(as.list(coef(fit0)), c(a = Asym, b = Asym - R0, c = exp(lrc)))
nls(y ~ a - b * exp(-c*x), data = df, start = st)

给予:

Nonlinear regression model
  model: y ~ a - b * exp(-c * x)
   data: df
       a        b        c 
 0.38621 -0.21016  0.09033 
 residual sum-of-squares: 0.003971

Number of iterations to convergence: 0 
Achieved convergence tolerance: 6.531e-06

2) 正如所见,问题中的数据在 运行ning nls 中没有问题,但如果您对不同的数据确实有问题并想要第二个选项然后尝试 drc 包中的 AR.3 模型。下面的代码遵循上面的代码,除了我们进行更改以适应不同的参数化。

library(drc)
fit1 <- drm(y ~ x, data = df, fct = AR.3())
co <- setNames(coef(fit1), c("c", "d", "e"))
st <- with(as.list(co), c(a = d, b = d - c, c = 1/e))
nls(y ~ a - b * exp(-c*x), data = df, start = st)

更新

已修改。