如何找到非线性模型的起始值?
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)
如果你想要一些能够自动接近起点的东西(并且不想写一个自启动),你可以做一些假设:
- 假设 c 为正(即该图显示像半衰期曲线一样的衰减,如您的数据所示),那么
exp(-c * x)
将趋于零且 x 较大,因此只要您有在您的数据合理范围内,y
的最小值可能接近 a
- 如果我们从估计的
a
中减去 y 数据,b
的值与通过结果点进行线性回归的截距不会太远。
- 通过取这些新 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)
更新
已修改。
在 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)
如果你想要一些能够自动接近起点的东西(并且不想写一个自启动),你可以做一些假设:
- 假设 c 为正(即该图显示像半衰期曲线一样的衰减,如您的数据所示),那么
exp(-c * x)
将趋于零且 x 较大,因此只要您有在您的数据合理范围内,y
的最小值可能接近a
- 如果我们从估计的
a
中减去 y 数据,b
的值与通过结果点进行线性回归的截距不会太远。 - 通过取这些新 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)
更新
已修改。