用负 y 值拟合指数衰减

机器算法验证 r 模型非线性回归 nls

2022-03-09 10:49:49

我正在尝试将指数衰减函数拟合到在高 x 值时变为负数的 y 值，但无法nls正确配置我的函数。

目标

我对衰减函数的斜率（根据某些来源）感兴趣。我如何得到这个斜率并不重要，但模型应该尽可能适合我的数据（即，如果拟合良好，线性化问题是可以接受的；请参阅“线性化”）。然而，之前关于该主题的作品使用了以下指数衰减函数（Stedmon 等人的封闭访问文章，等式 3）： $\lambda$

$f(y) = a \times exp(-S \times x) + K$

S我感兴趣的斜率在哪里K，允许负值的校正因子和（即截距）a的初始值。x

我需要在 R 中执行此操作，因为我正在编写一个函数，它将发色溶解有机物 (CDOM) 的原始测量值转换为研究人员感兴趣的值。

示例数据

由于数据的性质，我不得不使用 PasteBin。示例数据可在此处获得。

dt <-将 PasteBin 中的代码编写并复制到您的 R 控制台。IE

dt <- structure(list(x = ...

数据如下所示：

library(ggplot2)
ggplot(dt, aes(x = x, y = y)) + geom_point()

时出现负 y 值。 $x > 540 nm$

试图找到解决方案`nls`

最初的尝试使用nls产生了一个奇点，看到我刚刚看到参数的起始值，这应该不足为奇：

nls(y ~ a * exp(-S * x) + K, data = dt, start = list(a = 0.5, S = 0.1, K = -0.1))

# Error in nlsModel(formula, mf, start, wts) : 
#  singular gradient matrix at initial parameter estimates

按照这个答案，我可以尝试制作更好的拟合启动参数来帮助该nls功能：

K0 <- min(dt$y)/2
mod0 <- lm(log(y - K0) ~ x, data = dt) # produces NaNs due to the negative values
start <- list(a = exp(coef(mod0)[1]), S = coef(mod0)[2], K = K0)
nls(y ~ a * exp(-S * x) + K, data = dt, start = start)

# Error in nls(y ~ a * exp(-S * x) + K, data = dt, start = start) : 
#  number of iterations exceeded maximum of 50

该函数似乎无法找到具有默认迭代次数的解决方案。让我们增加迭代次数：

nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000))

# Error in nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000)) : 
#  step factor 0.000488281 reduced below 'minFactor' of 0.000976562

更多错误。扔掉它！让我们强制函数给我们一个解决方案：

mod <- nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000, warnOnly = TRUE))
mod.dat <- data.frame(x = dt$x, y = predict(mod, list(wavelength = dt$x)))

ggplot(dt, aes(x = x, y = y)) + geom_point() + 
  geom_line(data = mod.dat, aes(x = x, y = y), color = "red")

好吧，这绝对不是一个好的解决方案......

线性化问题

许多人已经成功地将他们的指数衰减函数线性化（来源：1、2、3）。在这种情况下，我们需要确保没有 y 值是负数或 0。让我们在计算机的浮点限制内使最小 y 值尽可能接近 0 ：

K <- abs(min(dt$y)) 
dt$y <- dt$y + K*(1+10^-15)

fit <- lm(log(y) ~ x, data=dt)  
ggplot(dt, aes(x = x, y = y)) + geom_point() + 
geom_line(aes(x=x, y=exp(fit$fitted.values)), color = "red")

好多了，但是模型在低 x 值时不能完美地跟踪 y 值。

请注意，该nls函数仍然无法适应指数衰减：

K0 <- min(dt$y)/2
mod0 <- lm(log(y - K0) ~ x, data = dt) # produces NaNs due to the negative values
start <- list(a = exp(coef(mod0)[1]), S = coef(mod0)[2], K = K0)
nls(y ~ a * exp(-S * x) + K, data = dt, start = start)

# Error in nlsModel(formula, mf, start, wts) : 
#  singular gradient matrix at initial parameter estimates

负值重要吗？

负值显然是测量误差，因为吸收系数不能为负。那么，如果我将 y 值设为正数呢？是我感兴趣的坡度。如果加法不影响坡度，我应该解决：

dt$y <- dt$y + 0.1

fit <- lm(log(y) ~ x, data=dt)  
ggplot(dt, aes(x = x, y = y)) + geom_point() + geom_line(aes(x=x, y=exp(fit$fitted.values)), color = "red")

好吧，这并没有那么顺利......高 x 值显然应该尽可能接近零。

问题

我显然在这里做错了什么。使用 R 估计拟合在具有负 y 值的数据上的指数衰减函数的斜率的最准确方法是什么？

2个回答

使用自启动功能：

ggplot(dt, aes(x = x, y = y)) + 
  geom_point() +
  stat_smooth(method = "nls", formula = y ~ SSasymp(x, Asym, R0, lrc), se = FALSE)

fit <- nls(y ~ SSasymp(x, Asym, R0, lrc), data = dt)
summary(fit)
#Formula: y ~ SSasymp(x, Asym, R0, lrc)
#
#Parameters:
#       Estimate Std. Error  t value Pr(>|t|)    
#Asym -0.0001302  0.0004693   -0.277    0.782    
#R0   77.9103278  2.1432998   36.351   <2e-16 ***
#lrc  -4.0862443  0.0051816 -788.604   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 0.007307 on 698 degrees of freedom
#
#Number of iterations to convergence: 0 
#Achieved convergence tolerance: 9.189e-08

exp(coef(fit)[["lrc"]]) #lambda
#[1] 0.01680222

但是，如果您的领域知识不能证明将渐近线设置为零，我会认真考虑。我相信它确实如此，并且上述模型并不不一致（参见系数的标准误差/p 值）。

ggplot(dt, aes(x = x, y = y)) + 
  geom_point() +
  stat_smooth(method = "nls", formula = y ~ a * exp(-S * x), 
              method.args = list(start = list(a = 78, S = 0.02)), se = FALSE, #starting values obtained from fit above
              color = "dark red")

这个问题与其他几个问题有关系

关于这个问题中的一些观点，我还有三点补充意见。

1：为什么线性化模型不能很好地拟合 $y$

好多了，但是模型在低 x 值时不能完美地跟踪 y 值。

线性化拟合不会最小化相同的残差。在对数尺度下，较小值的残差将较大。下图通过在右图中的对数刻度上绘制 y 轴来显示比较：

必要时，您可以向最小二乘损失函数添加权重。

2：使用线性化拟合作为起始值

在您通过线性拟合获得估计后，您可以将这些作为非线性拟合的起点。*

# vectors x and y from data
x <- dat$x
y <- dat$y

# linearized fit with zero correction
K <- abs(min(y)) 
dty <- y + K*(1+10^-15)
fit <- lm(log(dty) ~x)  


# old fit that had a singluar gradient matrix error
#         nls(y ~ a * exp(-S * x) + K, 
#                 start = list(a = 0.5, 
#                              S = 0.1, 
#                              K = -0.1))
#

# new fit
fitnls <- nls(y ~ a * exp(-S * x) + K, 
                  start = list(a = exp(fit$coefficients[1]), 
                               S = -fit$coefficients[2], 
                               K = -0.1))
#

3：使用更通用的方法获取起点

如果你有足够的点，那么你也可以得到斜率，而不必担心渐近值和负值（不需要计算对数）。

您可以通过整合数据点来做到这一点。然后使用和描述为线性组合来获得的值向量、和截距：

y = a e^{s x} + k

$y = a e^{sx} + k$

Y = \frac{a}{s} e^{s x} + k x + C o n s t

$Y = \frac{a}{s} e^{sx} + kx + Const$

s

$s$

y

$y$

Y

$Y$

x

$x$

\begin{array}{rccccl} y & = & a e^{s x} + k & = & s (\frac{a}{s} e^{s x} + k x + C o n s t) & - s k x - s C o n s t \\ = & s Y & - s k x - s C o n s t \end{array}

$\begin{array}{rccccl}y &=& a e^{sx} + k &=& s(\frac{a}{s} e^{s x} + k x + Const) &- s k x - s Const \\ &&&=& sY &- sk x - s Const \end{array}$

这种方法的优点（参见 Tittelbach 和 Helmrich 1993 “一种用于分析多指数瞬态信号的积分方法”）是您可以将其扩展到多个指数衰减分量（添加更多积分）。

#
# using Tittelbach Helmrich
#

# integrating with trapezium rule assuming x variable is already ordered
ys <- c(0,cumsum(0.5*diff(x)*(y[-1]+y[-length(y)])))

# getting slope parameter
modth <- lm(y ~ ys + x)
slope <- modth$coefficients[2]

# getting other parameters 
modlm <- lm(y ~ 1 + I(exp(slope*x)))
K <- modlm$coefficients[1]
a <- modlm$coefficients[2]

# fitting with TH start

fitnls2 <- nls(y ~ a * exp(-S * x) + K, 
              start = list(a = a, 
                           S = -slope, 
                           K = K))

脚注： *在线性化问题中使用斜率正是自SSasymp启动函数的作用。它首先估计渐近线

> stats:::NLSstRtAsymptote.sortedXyData
function (xy) 
{
    in.range <- range(xy$y)
    last.dif <- abs(in.range - xy$y[nrow(xy)])
    if (match(min(last.dif), last.dif) == 2L) 
        in.range[2L] + diff(in.range)/8
    else in.range[1L] - diff(in.range)/8
}

然后是斜率（减去渐近线值并取对数值）

> stats:::NLSstAsymptotic.sortedXyData
function (xy) 
{
    xy$rt <- NLSstRtAsymptote(xy)
    setNames(coef(nls(y ~ cbind(1, 1 - exp(-exp(lrc) * x)), data = xy, 
        start = list(lrc = log(-coef(lm(log(abs(y - rt)) ~ x, 
            data = xy))[[2L]])), algorithm = "plinear"))[c(2, 
        3, 1)], c("b0", "b1", "lrc"))
}

注意线start = list(lrc = log(-coef(lm(log(abs(y - rt)) ~ x, data = xy))[[2L]]))

旁注：在的特殊情况下，您可以使用 $K=0$

plot(x,y)
mod <- glm(y~x, family = gaussian(link = log), start = c(2,-0.01))
lines(x,exp(predict(mod)),col=2)

它将观察到的参数建模为 $y$

y = e x p (X β) + ϵ = e x p (β_{0}) \cdot e x p (β_{1} \cdot x) + ϵ

$y = exp(X\beta) + \epsilon = exp(\beta_0) \cdot exp(\beta_1 \cdot x) + \epsilon$

其它你可能感兴趣的问题

上一篇一步一步的普通克里金例子？下一篇关于正则化方法的一些最重要的“早期论文”是什么？