r - R中的指数曲线拟合

Question

time = 1:100  
head(y)  
0.07841589 0.07686316 0.07534116 0.07384931 0.07238699 0.07095363   
plot(time,y)  

这是一条指数曲线。
1)如何在不知道公式的情况下在这条曲线上拟合线？我不能使用“nls”，因为公式未知(只给出了数据点)。
2) 我怎样才能得到这条曲线的方程并确定常数。在等式中？
我试过 loess 但它没有给出拦截

Answer 1

注:这个答案已经完全改写自 original ,
这在几个方面存在缺陷(感谢评论者强调这些)。我希望这个新答案是正确的。

您需要一个模型来拟合数据。
在不知道模型的全部细节的情况下，假设这是一个
exponential growth model ,
哪一个可以写成:y = a * e r*t

其中 y 是您的测量变量，t 是测量它的时间，
a 是 t = 0 时 y 的值，r 是增长常数。
我们要估计 a 和 r。

这是一个非线性问题，因为我们要估计指数 r。
但是，在这种情况下，我们可以使用一些代数并将其转换为线性方程，方法是对两边取对数并求解(记住
logarithmic rules )， 导致:
log(y) = log(a) + r * t

我们可以用一个例子来形象化，通过从我们的模型中生成一条曲线，假设 a 和 r 的一些值:

t <- 1:100      # these are your time points
a <- 10         # assume the size at t = 0 is 10
r <- 0.1        # assume a growth constant
y <- a*exp(r*t) # generate some y observations from our exponential model

# visualise
par(mfrow = c(1, 2))
plot(t, y)      # on the original scale
plot(t, log(y)) # taking the log(y)

因此，对于这种情况，我们可以探索两种可能性:

将我们的非线性模型拟合到原始数据(例如使用 nls() 函数)

将我们的“线性化”模型拟合到对数转换数据(例如使用 lm() 函数)

选择哪个选项(还有更多选项)，取决于我们的想法
(或假设)是我们数据背后的数据生成过程。

让我们用一些模拟来说明，其中包括添加的噪声(采样自
正态分布)，以模拟真实数据。请看这个
StackExchange post
对于此模拟背后的推理(由 Alejo Bernardin's comment 指出)。

set.seed(12) # for reproducible results

# errors constant across time - additive
y_add <- a*exp(r*t) + rnorm(length(t), sd = 5000) # or: rnorm(length(t), mean = a*exp(r*t), sd = 5000)

# errors grow as y grows - multiplicative (constant on the log-scale)
y_mult <- a*exp(r*t + rnorm(length(t), sd = 1))  # or: rlnorm(length(t), mean = log(a) + r*t, sd = 1)

# visualise
par(mfrow = c(1, 2))
plot(t, y_add, main = "additive error")
lines(t, a*exp(t*r), col = "red") 
plot(t, y_mult, main = "multiplicative error")
lines(t, a*exp(t*r), col = "red")

对于加法模型，我们可以使用 nls() ，因为误差是恒定的
吨。使用时 nls()我们需要为优化算法指定一些起始值(尝试“猜测”这些是什么，因为 nls() 经常难以收敛于解决方案)。

add_nls <- nls(y_add ~ a*exp(r*t), 
               start = list(a = 0.5, r = 0.2))
coef(add_nls)

#           a           r 
# 11.30876845  0.09867135 

使用 coef()函数我们可以得到两个参数的估计值。
这给了我们可以的估计，接近我们模拟的(a = 10 和 r = 0.1)。

通过绘制模型的残差，您可以看到误差方差在数据范围内是合理恒定的:

plot(t, resid(add_nls))
abline(h = 0, lty = 2)

对于乘法误差情况(我们的 y_mult 模拟值)，我们应该使用 lm()在对数转换的数据上，因为
相反，该错误在该范围内是恒定的。

mult_lm <- lm(log(y_mult) ~ t)
coef(mult_lm)

# (Intercept)           t 
#  2.39448488  0.09837215 

为了解释这个输出，再次记住我们的线性模型是 log(y) = log(a) + r*t，它等价于形式为 Y = β0 + β1 * X 的线性模型，其中 β0 是我们的截距， β1 我们的斜率。

因此，在此输出中 (Intercept)相当于我们模型的 log(a) 和 t是时间变量的系数，所以等价于我们的 r。
有意义地解释 (Intercept)我们可以取它的指数 ( exp(2.39448488) )，给我们 ~10.96，这非常接近我们的模拟值。

值得注意的是，如果我们拟合误差为乘法的数据会发生什么
使用 nls函数代替:

mult_nls <- nls(y_mult ~ a*exp(r*t), start = list(a = 0.5, r = 0.2))
coef(mult_nls)

#            a            r 
# 281.06913343   0.06955642 

现在我们高估了 a 并低估了 r
( Mario Reutter
在他的评论中强调了这一点)。我们可以想象使用错误的方法来拟合我们的模型的后果:

# get the model's coefficients
lm_coef <- coef(mult_lm)
nls_coef <- coef(mult_nls)

# make the plot
plot(t, y_mult)
lines(t, a*exp(r*t), col = "brown", lwd = 5)
lines(t, exp(lm_coef[1])*exp(lm_coef[2]*t), col = "dodgerblue", lwd = 2)
lines(t, nls_coef[1]*exp(nls_coef[2]*t), col = "orange2", lwd = 2)
legend("topleft", col = c("brown", "dodgerblue", "orange2"), 
       legend = c("known model", "nls fit", "lm fit"), lwd = 3)

我们可以看到 lm()对数转换数据的拟合明显优于 nls()拟合原始数据。

您可以再次绘制该模型的残差，以查看方差在数据范围内不是恒定的(我们也可以在上图中看到这一点，其中数据的分布随着 t 值的增加而增加):

plot(t, resid(mult_nls))
abline(h = 0, lty = 2)