指数曲线在R中的拟合 [英] Exponential curve fitting in R

问题描述

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

这是一条指数曲线。

1. 在不知道公式的情况下如何在此曲线上拟合线？我无法使用"NLS"，因为公式未知(仅提供数据点)。

2. 如何获取此曲线的方程式并确定方程式中的常量？
   我试过LOWAY，但它不能进行拦截。

推荐答案

您需要一个适合数据的模型。 在不了解模型的全部细节的情况下，我们假设这是一个 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()函数)
• 将我们的线性化&q；模型与对数转换后的数据相匹配(例如，使用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()，因为误差是恒定的 t。使用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)

这篇关于指数曲线在R中的拟合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！