r - 如何在 R 中使用共享参数执行非线性最小二乘?
问题描述
我想在 R 中执行非线性最小二乘回归,同时最小化三个模型的平方残差(见下文)。现在,这三个模型共享一些参数,在我的示例中,参数b
和d
.
nls()
有没有办法用包minpack.lm
或包来做到这一点nlsr
?
因此,理想情况下,我想生成目标函数(所有模型的最小二乘之和)并一次回归所有参数:a1
、a2
、a3
、b
、c1
、c2
和。c3
d
(我试图避免运行三个独立的回归,然后对b
和进行一些平均d
。)
my_model <- function(x, a, b, c, d) {
a * b ^ (x - c) + d
}
# x values
x <- seq(0, 10, 0.2)
# Shared parameters
b <- 2
d <- 10
a1 <- 1
c1 <- 1
y1 <- my_model(x,
a = a1,
b = b,
c = c1,
d = d) + rnorm(length(x))
a2 <- 2
c2 <- 5
y2 <- my_model(x,
a = a2,
b = b,
c = c2,
d = d) + rnorm(length(x))
a3 <- -2
c3 <- 3
y3 <- my_model(x,
a = a3,
b = b,
c = c3,
d = d) + rnorm(length(x))
plot(
y1 ~ x,
xlim = range(x),
ylim = d + c(-50, 50),
type = 'b',
col = 'red',
ylab = 'y'
)
lines(y2 ~ x, type = 'b', col = 'green')
lines(y3 ~ x, type = 'b', col = 'blue')
解决方案
下面我们运行nls
(使用稍微修改的模型)和nlxb
(来自 nlsr)但nlxb
在收敛之前停止。尽管存在这些问题,但这两个问题确实给出了在视觉上很好地拟合数据的结果。这些问题表明模型本身存在问题,因此在“其他”部分中,根据nlxb
输出,我们展示了如何修复模型,给出原始模型的子模型,该子模型可以轻松地拟合数据nls
,nlxb
并且可以很好地拟合. 在注释部分的最后,我们以可重现的形式提供数据。
nls
假设在最后的注释中显示的设置可重现,通过定义一个右侧矩阵,其列分别与每个线性参数 a1、a2、a3 和 d 相乘,重新表述 nls plinear 算法的问题。plinear 不需要那些简化设置的起始值。它将分别报告为 .lin1、.lin2、.lin3 和 .lin4。
为了获得起始值,我们使用了一个没有分组的更简单的模型,并nls2
在同名包中使用从 1 到 10 的 b 和从 1 到 10 的 c 进行网格搜索。我们还发现nls
仍然会产生错误,但通过abs
在公式中使用,如图所示,它运行完成。
该模型的问题表明它存在根本问题,在“其他”部分中,我们将讨论如何修复它。
xx <- c(x, x, x)
yy <- c(y1, y2, y3)
# startingi values using nls2
library(nls2)
fo0 <- yy ~ cbind(b ^ abs(xx - c), 1)
st0 <- data.frame(b = c(1, 10), c = c(1, 10))
fm0 <- nls2(fo0, start = st0, alg = "plinear-brute")
# run nls using starting values from above
g <- rep(1:3, each = length(x))
fo <- yy ~ cbind((g==1) * b ^ abs(xx - c[g]),
(g==2) * b ^ abs(xx - c[g]),
(g==3) * b ^ abs(xx - c[g]),
1)
st <- with(as.list(coef(fm0)), list(b = b, c = c(c, c, c)))
fm <- nls(fo, start = st, alg = "plinear")
plot(yy ~ xx, col = g)
for(i in unique(g)) lines(predict(fm) ~ xx, col = i, subset = g == i)
fm
给予:
Nonlinear regression model
model: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx - c[g]), (g == 3) * b^abs(xx - c[g]), 1)
data: parent.frame()
b c1 c2 c3 .lin1 .lin2 .lin3 .lin4
1.997 0.424 1.622 1.074 0.680 0.196 -0.532 9.922
residual sum-of-squares: 133
Number of iterations to convergence: 5
Achieved convergence tolerance: 5.47e-06
(剧情后续)
nlsr
使用 nlsr 可以这样做。不需要网格搜索起始值,abs
也不需要添加。b 和 d 值似乎与 nls 解决方案相似,但其他系数不同。从视觉上看,这两种解决方案似乎都适合数据。
另一方面,从 JSingval 列中,我们看到 jacobian 秩不足,这导致它停止并且不产生 SE 值,并且收敛性存在疑问(尽管从视觉上看,未显示的情节似乎就足够了,看起来像一个很合适)。我们将在“其他”部分讨论如何解决此问题。
g1 <- g == 1; g2 <- g == 2; g3 <- g == 3
fo2 <- yy ~ g1 * (a1 * b ^ (xx - c1) + d) +
g2 * (a2 * b ^ (xx - c2) + d) +
g3 * (a3 * b ^ (xx - c3) + d)
st2 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, c1 = 1, c2 = 1, c3 = 1, d = 1)
fm2 <- nlxb(fo2, start = st2)
fm2
给予:
vn: [1] "yy" "g1" "a1" "b" "xx" "c1" "d" "g2" "a2" "c2" "g3" "a3" "c3"
no weights
nlsr object: x
residual sumsquares = 133.45 on 153 observations
after 16 Jacobian and 22 function evaluations
name coeff SE tstat pval gradient JSingval
a1 3.19575 NA NA NA 9.68e-10 4097
a2 0.64157 NA NA NA 8.914e-11 662.5
a3 -1.03096 NA NA NA -1.002e-09 234.9
b 1.99713 NA NA NA -2.28e-08 72.57
c1 2.66146 NA NA NA -2.14e-09 10.25
c2 3.33564 NA NA NA -3.955e-11 1.585e-13
c3 2.0297 NA NA NA -7.144e-10 1.292e-13
d 9.92363 NA NA NA -2.603e-12 3.271e-14
我们可以使用 nls2 作为第二阶段来计算 SE,但这仍然不能解决奇异值所暗示的整个问题。
summary(nls2(fo2, start = coef(fm2), algorithm = "brute-force"))
给予:
Formula: yy ~ g1 * (a1 * b^(xx - c1) + d) + g2 * (a2 * b^(xx - c2) + d) +
g3 * (a3 * b^(xx - c3) + d)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a1 3.20e+00 5.38e+05 0.0 1
a2 6.42e-01 3.55e+05 0.0 1
a3 -1.03e+00 3.16e+05 0.0 1
b 2.00e+00 2.49e-03 803.4 <2e-16 ***
c1 2.66e+00 9.42e-02 28.2 <2e-16 ***
c2 3.34e+00 2.43e+05 0.0 1
c3 2.03e+00 8.00e+05 0.0 1
d 9.92e+00 4.42e+05 0.0 1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.959 on 145 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: NA
其他
当nls
无法拟合模型时,它通常表明模型本身存在问题。在上面 nlsr 输出中的 JSingval 列的指导下,我们发现如果我们将所有c
参数值固定为 0,那么模型很容易拟合给定足够好的起始值和它仍然给出较低的残差平方和。d
c
library(nls2)
fo3 <- yy ~ cbind((g==1) * b ^ xx, (g==2) * b ^ xx, (g==3) * b ^ xx, 1)
st3 <- coef(fm0)["b"]
fm3 <- nls(fo3, start = st3, alg = "plinear")
给予:
Nonlinear regression model
model: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx, 1)
data: parent.frame()
b .lin1 .lin2 .lin3 .lin4
1.9971 0.5071 0.0639 -0.2532 9.9236
residual sum-of-squares: 133
Number of iterations to convergence: 4
Achieved convergence tolerance: 1.67e-09
fm
尽管参数少了 3 个,但以下 anova 表示的结果与上面的结果相当:
anova(fm3, fm)
给予:
Analysis of Variance Table
Model 1: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx, 1)
Model 2: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx - c[g]), (g == 3) * b^abs(xx - c[g]), 1)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 148 134
2 145 133 3 0.385 0.14 0.94
我们可以fm3
像nlxb
这样重做:
fo4 <- yy ~ g1 * (a1 * b ^ xx + d) +
g2 * (a2 * b ^ xx + d) +
g3 * (a3 * b ^ xx + d)
st4 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, d = 1)
fm4 <- nlxb(fo4, start = st4)
fm4
给予:
nlsr object: x
residual sumsquares = 133.45 on 153 observations
after 24 Jacobian and 33 function evaluations
name coeff SE tstat pval gradient JSingval
a1 0.507053 0.005515 91.94 1.83e-132 8.274e-08 5880
a2 0.0638554 0.0008735 73.11 4.774e-118 1.26e-08 2053
a3 -0.253225 0.002737 -92.54 7.154e-133 -4.181e-08 2053
b 1.99713 0.002294 870.6 2.073e-276 -2.55e-07 147.5
d 9.92363 0.09256 107.2 3.367e-142 -1.219e-11 10.26
笔记
下面假设的输入与问题中的相同,只是我们另外设置了种子以使其可重现。
set.seed(123)
my_model <- function(x, a, b, c, d) a * b ^ (x - c) + d
x <- seq(0, 10, 0.2)
b <- 2; d <- 10 # shared
a1 <- 1; c1 <- 1
y1 <- my_model(x, a = a1, b = b, c = c1, d = d) + rnorm(length(x))
a2 <- 2; c2 <- 5
y2 <- my_model(x, a = a2, b = b, c = c2, d = d) + rnorm(length(x))
a3 <- -2; c3 <- 3
y3 <- my_model(x, a = a3, b = b, c = c3, d = d) + rnorm(length(x))
推荐阅读
- javascript - JS - 将对象数组转换为对象,使新对象包含唯一的单词
- java - Java 线程被破坏
- python - Python FutureWarning:使用多个键进行索引(隐式转换为键元组)
- java - 使用“默认”作为变量
- java - 春季升级后@Value默认空值不起作用
- angular - 具有不同拦截器的新 HttpClient 实例
- javascript - 使用滑动滑块垂直滑块模式无法正确计算可见幻灯片计数
- android - 上传 2 张图片时无法访问的语句
- postgresql - 如何在postgres(sql文件)中的do块内设置变量
- apache-spark - 想要使用 spark 连接到 smb 服务器并在 spark 中从该服务器加载文件.. 让我们说