如何将对数回归拟合到 R 中的 "negative exponential" 散点图
How to fit logarithmic regression to a "negative exponential" scatterplot in R
我有每日降雨量 (x) 和观测值 (y) 的散点图,它看起来像 x^-2 图的 right/positive x 值的一半或以 1/ 为底的对数图2.基本上当 x 值非常低时 y 值非常高。 x 值越大,y 值越低。但是 y 值减小的速度变慢并且 y 永远不会为负。
这是一个代表性样本:
rain <- c(1, 1.2, 1.3, 2.5, 3.2, 4.2, 5, 7, 7.5, 10.3, 11.7, 12.9, 14.1, 15, 15.5, 17.5, 18.3, 20, 20.2, 20.3, 25, 28, 30, 34, 40)
obs <- c(42, 44, 43.9, 43.5, 35, 22, 18.4, 15.3, 10, 6.2, 5.7, 4, 3.7, 2.3, 2, 2.7, 3.5, 3, 2.9, 4, 1.6, 2.2, 1.6, 1.3, 0.8)
现在我想为这个散点图拟合一个回归模型。我已经尝试过 x^-4 之前的多项式回归,但我也想尝试对数回归,因为我认为它可能会成为更高质量的模型。
这是我到目前为止对多项式模型所做的:
y <- data$obs
x <- data$rain
xsq <- x^-2
xcub <- x^-3
xquar <- x^-4
fit4 <- lm(y~x+xsq+xcub+xquar) # I did the same for fit 1-3; until fit 4 it becomes more significant
xv <- seq(min(x), max(x), 0.01)
yv <- predict(fit5, list(x=xv, xsq=xv^-2, xcub=xv^-3, xquar=xv^-4))
lines(xv, yv)
这就是我对对数模型所做的尝试,但它只是 returns 与曲线不匹配的直线。感觉log()不是我真正需要的功能
xlog <- log(x)
fitlogx <- lm(y~xlog)
xv <- seq(min(xlog), max(xlog), 0.01)
yv <- predict(fitlogx, list(x=xv))
abline(fitlogx)
ylog <- log(y)
fitlogy <- lm(ylog~x)
xv <- seq(min(x), max(x), 0.01)
yv <- predict(fitlogy, list(x=xv))
abline(fitlogy)
现在我想知道如何拟合有意义的对数函数。如果您知道另一种可能有用的回归模型,我也非常感谢您提供任何建议。
您的 obs
变量非常适合 rain
的倒数。例如
dev.new(width=12, height=6)
oldp <- par(mfrow=c(1, 2))
plot(obs~rain)
lines(rain, 1/rain*40)
曲线需要高一点。我们可以反复猜测,例如尝试 rain*60,但使用 nls
函数更容易获得方程的最佳最小二乘拟合:
obs.nls <- nls(obs~1/rain*k, start=list(k=40))
summary(obs.nls)
#
# Formula: obs ~ 1/rain * k
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# k 57.145 4.182 13.66 8.12e-13 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 6.915 on 24 degrees of freedom
#
# Number of iterations to convergence: 1
# Achieved convergence tolerance: 2.379e-09
plot(obs~rain)
pred <- predict(obs.nls)
points(rain, pred, col="red", pch=18)
pred.rain <- seq(1, 40, length.out=100)
pred.obs <- predict(obs.nls, list(rain=pred.rain))
lines(pred.rain, pred.obs, col="blue", lty=2)
所以 k
的最佳估计值是 57.145。 nls
的主要缺点是您必须提供系数的起始值。它也可能无法收敛,但是对于我们在这里使用的简单函数,只要您可以估计合理的起始值,它就可以正常工作。
如果rain
可以有零值,您可以添加截距:
obs.nls <- nls(obs ~ k / (a + rain), start=list(a=1, k=40))
summary(obs.nls)
#
# Formula: obs ~ k/(a + rain)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a 1.4169 0.4245 3.337 0.00286 **
# k 117.5345 16.6878 7.043 3.55e-07 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 4.638 on 23 degrees of freedom
Number of iterations to convergence: 10
Achieved convergence tolerance: 6.763e-06
请注意,标准误差较小,但曲线高估了 rain > 10
的实际值。
我有每日降雨量 (x) 和观测值 (y) 的散点图,它看起来像 x^-2 图的 right/positive x 值的一半或以 1/ 为底的对数图2.基本上当 x 值非常低时 y 值非常高。 x 值越大,y 值越低。但是 y 值减小的速度变慢并且 y 永远不会为负。
这是一个代表性样本:
rain <- c(1, 1.2, 1.3, 2.5, 3.2, 4.2, 5, 7, 7.5, 10.3, 11.7, 12.9, 14.1, 15, 15.5, 17.5, 18.3, 20, 20.2, 20.3, 25, 28, 30, 34, 40)
obs <- c(42, 44, 43.9, 43.5, 35, 22, 18.4, 15.3, 10, 6.2, 5.7, 4, 3.7, 2.3, 2, 2.7, 3.5, 3, 2.9, 4, 1.6, 2.2, 1.6, 1.3, 0.8)
现在我想为这个散点图拟合一个回归模型。我已经尝试过 x^-4 之前的多项式回归,但我也想尝试对数回归,因为我认为它可能会成为更高质量的模型。
这是我到目前为止对多项式模型所做的:
y <- data$obs
x <- data$rain
xsq <- x^-2
xcub <- x^-3
xquar <- x^-4
fit4 <- lm(y~x+xsq+xcub+xquar) # I did the same for fit 1-3; until fit 4 it becomes more significant
xv <- seq(min(x), max(x), 0.01)
yv <- predict(fit5, list(x=xv, xsq=xv^-2, xcub=xv^-3, xquar=xv^-4))
lines(xv, yv)
这就是我对对数模型所做的尝试,但它只是 returns 与曲线不匹配的直线。感觉log()不是我真正需要的功能
xlog <- log(x)
fitlogx <- lm(y~xlog)
xv <- seq(min(xlog), max(xlog), 0.01)
yv <- predict(fitlogx, list(x=xv))
abline(fitlogx)
ylog <- log(y)
fitlogy <- lm(ylog~x)
xv <- seq(min(x), max(x), 0.01)
yv <- predict(fitlogy, list(x=xv))
abline(fitlogy)
现在我想知道如何拟合有意义的对数函数。如果您知道另一种可能有用的回归模型,我也非常感谢您提供任何建议。
您的 obs
变量非常适合 rain
的倒数。例如
dev.new(width=12, height=6)
oldp <- par(mfrow=c(1, 2))
plot(obs~rain)
lines(rain, 1/rain*40)
曲线需要高一点。我们可以反复猜测,例如尝试 rain*60,但使用 nls
函数更容易获得方程的最佳最小二乘拟合:
obs.nls <- nls(obs~1/rain*k, start=list(k=40))
summary(obs.nls)
#
# Formula: obs ~ 1/rain * k
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# k 57.145 4.182 13.66 8.12e-13 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 6.915 on 24 degrees of freedom
#
# Number of iterations to convergence: 1
# Achieved convergence tolerance: 2.379e-09
plot(obs~rain)
pred <- predict(obs.nls)
points(rain, pred, col="red", pch=18)
pred.rain <- seq(1, 40, length.out=100)
pred.obs <- predict(obs.nls, list(rain=pred.rain))
lines(pred.rain, pred.obs, col="blue", lty=2)
所以 k
的最佳估计值是 57.145。 nls
的主要缺点是您必须提供系数的起始值。它也可能无法收敛,但是对于我们在这里使用的简单函数,只要您可以估计合理的起始值,它就可以正常工作。
如果rain
可以有零值,您可以添加截距:
obs.nls <- nls(obs ~ k / (a + rain), start=list(a=1, k=40))
summary(obs.nls)
#
# Formula: obs ~ k/(a + rain)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a 1.4169 0.4245 3.337 0.00286 **
# k 117.5345 16.6878 7.043 3.55e-07 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 4.638 on 23 degrees of freedom
Number of iterations to convergence: 10
Achieved convergence tolerance: 6.763e-06
请注意,标准误差较小,但曲线高估了 rain > 10
的实际值。