nls 和有理函数的起始参数列表
Start list of parameters for nls and rational function
我想为我的数据拟合一个有理函数：
数据：
[1] 2.000000 3.000000 2.333333 1.750000 2.000000 1.833333 2.416667 1.916667
[9] 1.750000 2.166667 2.116667 1.916667 1.944444 1.611111 1.722222 1.777778
[17] 1.877778 1.944444 1.958333 1.833333 2.041667 2.020833 1.908333 1.916667
[25] 1.733333 1.833333 1.800000 1.933333 1.893333 1.866667 1.888889 1.805556
[33] 1.833333 1.847222 1.822222 1.805556 1.833333 1.904762 1.880952 1.833333
[41] 1.804762 1.809524 1.708333 1.708333 1.750000 1.708333 1.683333 1.687500
[49] 1.611111 1.666667 1.648148 1.611111 1.611111 1.611111 1.650000 1.600000
[57] 1.650000 1.625000 1.630000 1.616667 1.469697 1.560606 1.590909 1.651515
[65] 1.651515 1.651515 1.513889 1.555556 1.625000 1.638889 1.647222 1.652778
[73] 1.679487 1.717949 1.705128 1.698718.
我想要拟合的模型如下：
Model <- function(t, a, b, c, d) {
(a + b*t)/(1 + c*t + d*t^2)
}
I know that I firstly have to give a starting list of a, b, c... for nls but I really don't know how to set the parameters. Since I'm not an expert I found in this http://www.css.cornell.edu/faculty/dgr2/teach/R/R_rat.pdf document a useful guide. But at some point it says:
"Given a set of ordered pairs (ti,yi) where in general there are repeated measurements at each value of t, the parameters of a rational function can be fitted by non-linear least-squares estimation, for example with the nls method in R. One we have the four parameters, we can compute the value of t at which this is maximized, by computing the first derivative...."。
虽然我没有报告额外的数据，但我还有另一列代表时间（来自 1:76 的整数代表年份）。
谁能帮帮我？
最佳
E.
我相信你可以做类似的事情
nls(y ~ (a1 + b1*times) / (1 + c1*times + d1*times^2))
其中 y 是您在上面提供的数据，times=1:76。我在参数中添加了一个 1，因为 nls 没有将 c 识别为参数，而是将 c().
识别为函数
但是，当我 运行 这样做时，我收到 singular gradient 错误，并建议将起始值初始化为 1（默认值）以外的值。您可以使用参数 start = list("a1"=0.1, "b1"=0.1, "c1"=0.1, "d1"=0.1) 指定起始值，但这似乎没有帮助。也许您对起始值应该是多少有更好的了解？
问题中没有完全指定模型，但假设下面代码中的模型和下面注释 2 中可重复显示的数据如果我们设置 c = d = 0 那么它是一个线性模型，所以我们可以使用来自线性模型的系数作为起始值拟合：
fm1 <- lm(y ~ t)
st2 <- list(a = coef(fm1)[[1]], b = coef(fm1)[[2]], c = 0, d = 0)
fm2 <- nls(y ~ Model(t, a, b, c, d), start = st2)
给予：
> fm2
Nonlinear regression model
model: y ~ Model(t, a, b, c, d)
data: parent.frame()
a b c d
2.5097712 0.6038808 0.3205409 0.0008663
residual sum-of-squares: 1.684
Number of iterations to convergence: 16
Achieved convergence tolerance: 8.029e-06
以图形方式查看拟合：
# model is shown in red. See Note 1 for fm4 (blue) and fm5 (green) models.
plot(y ~ t)
lines(fitted(fm2) ~ t, col = "red")
lines(fitted(fm4) ~ t, col = "blue")
lines(fitted(fm5) ~ t, col = "green")
legend("topright", c("fm2", "fm4", "fm5"), col = c("red", "blue", "green"), lty = 1)
备注 1
以下是一个不同的模型，它几乎同样适用，但只使用了 3 个参数。见上图的蓝线。
fm3 <- lm(log(y) ~ log(t))
st4 <- list(a = coef(fm3)[[1]], b = 0, c = coef(fm3)[[2]])
fm4 <- nls(y ~ exp(a + b/t + c*log(t)), start = st4)
> fm4
Nonlinear regression model
model: y ~ exp(a + b/t + c * log(t))
data: parent.frame()
a b c
0.9845 -0.1767 -0.1157
residual sum-of-squares: 1.685
Number of iterations to convergence: 4
Achieved convergence tolerance: 2.625e-06
这个模型也不错。它仅使用两个参数，它们是线性的，残差平方和为 1.728837，而 fm2 模型为 1.684，fm4 模型为 1.685。请参见上图中的绿线。
fm5 <- lm(y ~ log(t))
> fm5
Call:
lm(formula = y ~ log(t))
Coefficients:
(Intercept) log(t)
2.4029 -0.1793
> deviance(fm5)
[1] 1.728837
注2
y <- c(2, 3, 2.333333, 1.75, 2, 1.833333, 2.416667, 1.916667, 1.75,
2.166667, 2.116667, 1.916667, 1.944444, 1.611111, 1.722222, 1.777778,
1.877778, 1.944444, 1.958333, 1.833333, 2.041667, 2.020833, 1.908333,
1.916667, 1.733333, 1.833333, 1.8, 1.933333, 1.893333, 1.866667,
1.888889, 1.805556, 1.833333, 1.847222, 1.822222, 1.805556, 1.833333,
1.904762, 1.880952, 1.833333, 1.804762, 1.809524, 1.708333, 1.708333,
1.75, 1.708333, 1.683333, 1.6875, 1.611111, 1.666667, 1.648148,
1.611111, 1.611111, 1.611111, 1.65, 1.6, 1.65, 1.625, 1.63, 1.616667,
1.469697, 1.560606, 1.590909, 1.651515, 1.651515, 1.651515, 1.513889,
1.555556, 1.625, 1.638889, 1.647222, 1.652778, 1.679487, 1.717949,
1.705128, 1.698718)
t <- seq_along(y)
我想为我的数据拟合一个有理函数：
数据：
[1] 2.000000 3.000000 2.333333 1.750000 2.000000 1.833333 2.416667 1.916667
[9] 1.750000 2.166667 2.116667 1.916667 1.944444 1.611111 1.722222 1.777778
[17] 1.877778 1.944444 1.958333 1.833333 2.041667 2.020833 1.908333 1.916667
[25] 1.733333 1.833333 1.800000 1.933333 1.893333 1.866667 1.888889 1.805556
[33] 1.833333 1.847222 1.822222 1.805556 1.833333 1.904762 1.880952 1.833333
[41] 1.804762 1.809524 1.708333 1.708333 1.750000 1.708333 1.683333 1.687500
[49] 1.611111 1.666667 1.648148 1.611111 1.611111 1.611111 1.650000 1.600000
[57] 1.650000 1.625000 1.630000 1.616667 1.469697 1.560606 1.590909 1.651515
[65] 1.651515 1.651515 1.513889 1.555556 1.625000 1.638889 1.647222 1.652778
[73] 1.679487 1.717949 1.705128 1.698718.
我想要拟合的模型如下：
Model <- function(t, a, b, c, d) { (a + b*t)/(1 + c*t + d*t^2) }I know that I firstly have to give a starting list of a, b, c... for nls but I really don't know how to set the parameters. Since I'm not an expert I found in this http://www.css.cornell.edu/faculty/dgr2/teach/R/R_rat.pdf document a useful guide. But at some point it says:
"Given a set of ordered pairs (ti,yi) where in general there are repeated measurements at each value of t, the parameters of a rational function can be fitted by non-linear least-squares estimation, for example with the nls method in R. One we have the four parameters, we can compute the value of t at which this is maximized, by computing the first derivative...."。
虽然我没有报告额外的数据，但我还有另一列代表时间（来自 1:76 的整数代表年份）。
谁能帮帮我？
最佳
E.
我相信你可以做类似的事情
nls(y ~ (a1 + b1*times) / (1 + c1*times + d1*times^2))
其中 y 是您在上面提供的数据，times=1:76。我在参数中添加了一个 1，因为 nls 没有将 c 识别为参数，而是将 c().
但是，当我 运行 这样做时，我收到 singular gradient 错误，并建议将起始值初始化为 1（默认值）以外的值。您可以使用参数 start = list("a1"=0.1, "b1"=0.1, "c1"=0.1, "d1"=0.1) 指定起始值，但这似乎没有帮助。也许您对起始值应该是多少有更好的了解？
问题中没有完全指定模型，但假设下面代码中的模型和下面注释 2 中可重复显示的数据如果我们设置 c = d = 0 那么它是一个线性模型，所以我们可以使用来自线性模型的系数作为起始值拟合：
fm1 <- lm(y ~ t)
st2 <- list(a = coef(fm1)[[1]], b = coef(fm1)[[2]], c = 0, d = 0)
fm2 <- nls(y ~ Model(t, a, b, c, d), start = st2)
给予：
> fm2
Nonlinear regression model
model: y ~ Model(t, a, b, c, d)
data: parent.frame()
a b c d
2.5097712 0.6038808 0.3205409 0.0008663
residual sum-of-squares: 1.684
Number of iterations to convergence: 16
Achieved convergence tolerance: 8.029e-06
以图形方式查看拟合：
# model is shown in red. See Note 1 for fm4 (blue) and fm5 (green) models.
plot(y ~ t)
lines(fitted(fm2) ~ t, col = "red")
lines(fitted(fm4) ~ t, col = "blue")
lines(fitted(fm5) ~ t, col = "green")
legend("topright", c("fm2", "fm4", "fm5"), col = c("red", "blue", "green"), lty = 1)
备注 1
以下是一个不同的模型，它几乎同样适用，但只使用了 3 个参数。见上图的蓝线。
fm3 <- lm(log(y) ~ log(t))
st4 <- list(a = coef(fm3)[[1]], b = 0, c = coef(fm3)[[2]])
fm4 <- nls(y ~ exp(a + b/t + c*log(t)), start = st4)
> fm4
Nonlinear regression model
model: y ~ exp(a + b/t + c * log(t))
data: parent.frame()
a b c
0.9845 -0.1767 -0.1157
residual sum-of-squares: 1.685
Number of iterations to convergence: 4
Achieved convergence tolerance: 2.625e-06
这个模型也不错。它仅使用两个参数，它们是线性的，残差平方和为 1.728837，而 fm2 模型为 1.684，fm4 模型为 1.685。请参见上图中的绿线。
fm5 <- lm(y ~ log(t))
> fm5
Call:
lm(formula = y ~ log(t))
Coefficients:
(Intercept) log(t)
2.4029 -0.1793
> deviance(fm5)
[1] 1.728837
注2
y <- c(2, 3, 2.333333, 1.75, 2, 1.833333, 2.416667, 1.916667, 1.75,
2.166667, 2.116667, 1.916667, 1.944444, 1.611111, 1.722222, 1.777778,
1.877778, 1.944444, 1.958333, 1.833333, 2.041667, 2.020833, 1.908333,
1.916667, 1.733333, 1.833333, 1.8, 1.933333, 1.893333, 1.866667,
1.888889, 1.805556, 1.833333, 1.847222, 1.822222, 1.805556, 1.833333,
1.904762, 1.880952, 1.833333, 1.804762, 1.809524, 1.708333, 1.708333,
1.75, 1.708333, 1.683333, 1.6875, 1.611111, 1.666667, 1.648148,
1.611111, 1.611111, 1.611111, 1.65, 1.6, 1.65, 1.625, 1.63, 1.616667,
1.469697, 1.560606, 1.590909, 1.651515, 1.651515, 1.651515, 1.513889,
1.555556, 1.625, 1.638889, 1.647222, 1.652778, 1.679487, 1.717949,
1.705128, 1.698718)
t <- seq_along(y)