指数衰减适合 r
Exponential decay fit in r
我想将 R 中的指数衰减函数拟合到以下数据：
data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168,
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227,
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007,
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484,
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711,
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165,
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542,
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745,
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476,
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252,
0.109219168085624)), class = "data.frame", row.names = c(NA,
-39L), .Names = c("x", "y"))
我尝试过使用 nls 进行拟合，但生成的曲线与实际数据并不接近。
enter image description here
如果有人能解释如何处理此类非线性数据并找到最合适的函数，那将非常有帮助。
data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168,
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227,
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007,
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484,
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711,
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165,
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542,
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745,
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476,
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252,
0.109219168085624)), class = "data.frame", row.names = c(NA,
-39L), .Names = c("x", "y"))
# Do this because the log of 0 is not possible to calculate
data$x = data$x +1
fit = lm(log(y) ~ log(x), data = data)
plot(data$x, data$y)
lines(data$x, data$x ^ fit$coefficients[2], col = "red")
这比使用 nls forumla 好多了。并且在绘制拟合时似乎做得相当好。
双指数模型会更好，但仍不完美。这表明您可能有两个同时发生的衰变过程。
fit <- nls(y ~ SSbiexp(x, A1, lrc1, A2, lrc2), data = data)
#A1*exp(-exp(lrc1)*x)+A2*exp(-exp(lrc2)*x)
plot(y ~x, data = data)
curve(predict(fit, newdata = data.frame(x)), add = TRUE)
如果测量误差取决于量级，可以考虑用它来加权。
但是，您应该仔细考虑您期望从您的领域知识中得到什么样的模型。仅根据经验选择非线性模型通常不是一个好主意。非参数拟合可能是更好的选择。
尝试y ~ .lin / (b + x^c)。请注意，当使用 "plinear" 时，在将公式指定为 nls 时会省略 .lin 线性参数，并且还会省略它的起始值。
另请注意，.lin 和 b 参数在最佳状态下约为 1，因此我们也可以尝试单参数模型 y ~ 1 / (1 + x^c)。这是单参数对数逻辑生存曲线的形式。这个单参数模型的 AIC 比 3 参数模型差（比较 AIC(fm1) 和 AIC(fm3)），但单参数模型可能仍然更可取，因为它的简约性和视觉上的拟合与 3 参数模型无法区分。
opar <- par(mfcol = 2:1, mar = c(3, 3, 3, 1), family = "mono")
# data = data.frame with x & y col names; fm = model fit; main = string shown above plot
Plot <- function(data, fm, main) {
plot(y ~ x, data, pch = 20)
lines(fitted(fm) ~ x, data, col = "red")
legend("topright", bty = "n", cex = 0.7, legend = capture.output(fm))
title(main = paste(main, "- AIC:", round(AIC(fm), 2)))
}
# 3 parameter model
fo3 <- y ~ 1/(b + x^c) # omit .lin parameter; plinear will add it automatically
fm3 <- nls(fo3, data = data, start = list(b = 1, c = 1), alg = "plinear")
Plot(data, fm3, "3 parameters")
# one parameter model
fo1 <- y ~ 1 / (1 + x^c)
fm1 <- nls(fo1, data, start = list(c = 1))
Plot(data, fm1, "1 parameter")
par(read.only = opar)
AIC
在其他答案中添加解决方案我们可以比较 AIC 值。我们用它使用的参数数量标记了每个解决方案（自由度会比那个大一个）并且重新设计了对数对数解决方案以使用 nls 而不是 lm 并且有一个 LHS y 因为人们无法比较具有不同左侧或使用不同优化例程的模型的 AIC 值，因为使用的对数似然常数可能不同。
fo2 <- y ~ exp(a + b * log(x+1))
fm2 <- nls(fo2, data, start = list(a = 1, b = 1))
fo4 <- y ~ SSbiexp(x, A1, lrc1, A2, lrc2)
fm4 <- nls(fo4, data)
aic <- AIC(fm1, fm2, fm3, fm4)
aic[order(aic$AIC), ]
从最好的 AIC（即 fm3）到最差的 AIC（即 fm2）：
df AIC
fm3 4 -329.35
fm1 2 -307.69
fm4 5 -215.96
fm2 3 -167.33
我想将 R 中的指数衰减函数拟合到以下数据：
data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168,
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227,
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007,
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484,
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711,
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165,
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542,
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745,
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476,
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252,
0.109219168085624)), class = "data.frame", row.names = c(NA,
-39L), .Names = c("x", "y"))
我尝试过使用 nls 进行拟合，但生成的曲线与实际数据并不接近。
enter image description here
如果有人能解释如何处理此类非线性数据并找到最合适的函数，那将非常有帮助。
data <- structure(list(x = 0:38, y = c(0.991744340878828, 0.512512332368168,
0.41102449265681, 0.356621905557202, 0.320851602373477, 0.29499198506227,
0.275037747162642, 0.25938850981822, 0.245263623938863, 0.233655093612007,
0.224041426946405, 0.214152907133301, 0.207475138903635, 0.203270738895484,
0.194942528735632, 0.188107106969046, 0.180926819430008, 0.177028560207711,
0.172595416846822, 0.166729221891201, 0.163502461048814, 0.159286528409165,
0.156110097827889, 0.152655498715612, 0.148684858095915, 0.14733605355542,
0.144691873223729, 0.143118852619617, 0.139542186417186, 0.137730138713745,
0.134353615271572, 0.132197800438632, 0.128369567159113, 0.124971834736476,
0.120027536018095, 0.117678812415655, 0.115720611113327, 0.112491329844252,
0.109219168085624)), class = "data.frame", row.names = c(NA,
-39L), .Names = c("x", "y"))
# Do this because the log of 0 is not possible to calculate
data$x = data$x +1
fit = lm(log(y) ~ log(x), data = data)
plot(data$x, data$y)
lines(data$x, data$x ^ fit$coefficients[2], col = "red")
这比使用 nls forumla 好多了。并且在绘制拟合时似乎做得相当好。
双指数模型会更好，但仍不完美。这表明您可能有两个同时发生的衰变过程。
fit <- nls(y ~ SSbiexp(x, A1, lrc1, A2, lrc2), data = data)
#A1*exp(-exp(lrc1)*x)+A2*exp(-exp(lrc2)*x)
plot(y ~x, data = data)
curve(predict(fit, newdata = data.frame(x)), add = TRUE)
如果测量误差取决于量级，可以考虑用它来加权。
但是，您应该仔细考虑您期望从您的领域知识中得到什么样的模型。仅根据经验选择非线性模型通常不是一个好主意。非参数拟合可能是更好的选择。
尝试y ~ .lin / (b + x^c)。请注意，当使用 "plinear" 时，在将公式指定为 nls 时会省略 .lin 线性参数，并且还会省略它的起始值。
另请注意，.lin 和 b 参数在最佳状态下约为 1，因此我们也可以尝试单参数模型 y ~ 1 / (1 + x^c)。这是单参数对数逻辑生存曲线的形式。这个单参数模型的 AIC 比 3 参数模型差（比较 AIC(fm1) 和 AIC(fm3)），但单参数模型可能仍然更可取，因为它的简约性和视觉上的拟合与 3 参数模型无法区分。
opar <- par(mfcol = 2:1, mar = c(3, 3, 3, 1), family = "mono")
# data = data.frame with x & y col names; fm = model fit; main = string shown above plot
Plot <- function(data, fm, main) {
plot(y ~ x, data, pch = 20)
lines(fitted(fm) ~ x, data, col = "red")
legend("topright", bty = "n", cex = 0.7, legend = capture.output(fm))
title(main = paste(main, "- AIC:", round(AIC(fm), 2)))
}
# 3 parameter model
fo3 <- y ~ 1/(b + x^c) # omit .lin parameter; plinear will add it automatically
fm3 <- nls(fo3, data = data, start = list(b = 1, c = 1), alg = "plinear")
Plot(data, fm3, "3 parameters")
# one parameter model
fo1 <- y ~ 1 / (1 + x^c)
fm1 <- nls(fo1, data, start = list(c = 1))
Plot(data, fm1, "1 parameter")
par(read.only = opar)
AIC
在其他答案中添加解决方案我们可以比较 AIC 值。我们用它使用的参数数量标记了每个解决方案（自由度会比那个大一个）并且重新设计了对数对数解决方案以使用 nls 而不是 lm 并且有一个 LHS y 因为人们无法比较具有不同左侧或使用不同优化例程的模型的 AIC 值，因为使用的对数似然常数可能不同。
fo2 <- y ~ exp(a + b * log(x+1))
fm2 <- nls(fo2, data, start = list(a = 1, b = 1))
fo4 <- y ~ SSbiexp(x, A1, lrc1, A2, lrc2)
fm4 <- nls(fo4, data)
aic <- AIC(fm1, fm2, fm3, fm4)
aic[order(aic$AIC), ]
从最好的 AIC（即 fm3）到最差的 AIC（即 fm2）：
df AIC
fm3 4 -329.35
fm1 2 -307.69
fm4 5 -215.96
fm2 3 -167.33