编码非线性回归、指数衰减

Question

还是 R 的新手。

背景

我有 Li 等人的数据。 2003年论文“加拿大林业部门碳预算模型中的地下生物量动态”。(https:/ /cdnsciencepub.com/doi/10.1139/x02-165）我正在尝试在 R 中重新创建方程 6，以便我可以生成相同的参数估计。 （等式 6），然而，Pf 的最大值为 0.426，因为更大的值将产生 0。

为了给函数施加最大值，我们可以对参数施加限制。最简单的是我是否可以根据最大值重新参数化函数。在 RB>0 且指数项衰减的情况下，我们对函数 Pf(RB)=a+bexp(cRB) 有最大值

m = max(Pf) = a+b

所以这可以重写为

a = m − b

并拟合方程

Pf(RB) = m − b + b·exp(c·RB) = m(1 − b·exp(c·RB))

我想使用 r 中的 nls 对其进行编码，但我不知道如何进行继续以上限来做这件事。我正在尝试在 R 中生成模型，以便我可以获得表 3 方程 6 中生成的相同参数。任何帮助将不胜感激。

Answer 1

我们可以按如下方式使用nls。

fit = nls(fine_prop ~ m * (1 - b * exp(c*total_root)),
    data = fr,
    start = list(m=0.1, b = 2, c=-0.05))

coefficients(fit)
#          m           b           c 
# 0.06851410 -5.06265496 -0.05369425 

plot(fr)
x = seq(0, max(fr$total_root), length.out=100)
lines(x, predict(fit, newdata = data.frame(total_root = x)))

同样，我们也可以将方程采用与已发表文章中相同的格式，而不是OP中的形式。我们有：

fit = nls(fine_prop ~ a + b * exp(c*total_root),
    data = fr,
    start = list(a=1, b = 1, c=-0.1))

coefficients(fit)
#          a          b          c 
#  0.0685124  0.3468601 -0.0536925

数据：

fr = structure(list(total_root = c(28, 47.5, 104.5, 62, 59.8304, 50.9335, 
40.1412, 43.3121, 131.0057, 24.74, 137.71, 25.004, 88.1, 88.1, 
57.6, 57.6, 54.16, 82.82, 88.6, 134.68, 20.92, 25.004, 141.31, 
143.8, 204.04, 104.88, 172.8, 122.62, 157.7, 18.184, 13.538, 
10.4, 27, 36.7, 44.2, 53.9, 78.6, 47.5, 6.6, 10.1, 14.6, 16.8, 
18.3, 8, 9.5, 6.2, 14.1, 15.8, 21.6, 29.1, 33.2, 41, 45, 46, 
59, 37.6, 9.3, 15, 22.7, 43.2, 37.1, 51, 28.6, 60, 75.2, 12.2, 
43.8, 43.3, 6.8, 6, 7.1, 8.9, 94.3, 100.5, 44.9, 9.98, 9.17, 
7.38, 20.9, 9.6, 11.48, 22.68, 27.63, 27.63, 15.19, 27.78, 26.63, 
12.44, 16.7, 0.5, 3.13), fine_prop = c(0.033796, 0.037486, 0.033818, 
0.033774, 0.10147, 0.032788, 0.068508, 0.045253, 0.005412, 0.373484, 
0.092876, 0.196049, 0.030647, 0.051078, 0.144097, 0.182292, 0.182422, 
0.051195, 0.059594, 0.113306, 0.291587, 0.141337, 0.115562, 0.0758, 
0.053911, 0.075324, 0.075231, 0.08563, 0.061509, 0.125825, 0.130743, 
0.203846, 0.37037, 0.174387, 0.144796, 0.079777, 0.075064, 0.042316, 
0.090909, 0.059406, 0.047945, 0.039881, 0.040984, 0.0875, 0.063158, 
0.091935, 0.061702, 0.063924, 0.057407, 0.052234, 0.052711, 0.047317, 
0.042667, 0.043913, 0.033898, 0.035106, 0.26129, 0.28, 0.129956, 
0.306019, 0.190566, 0.117647, 0.158042, 0.106667, 0.142553, 0.142623, 
0.025571, 0.010624, 0.633824, 0.616667, 0.56338, 0.278652, 0.059597, 
0.013433, 0.088864, 0.305611, 0.314068, 0.298103, 0.398086, 0.491667, 
0.408537, 0.260582, 0.299312, 0.103511, 0.357472, 0.236501, 0.241833, 
0.276527, 0.353293, 0.4, 0.389776)), row.names = c(NA, -91L), class = "data.frame")