从QR分解获得帽矩阵加权最小二乘回归

我想扩展包McSptial的lwr()函数，它适合作为非参数估计的加权回归。在lwr()函数的核心中，它使用solve()而不是QR分解来反转矩阵，导致数值不稳定。我想改变它，但不知道如何从QR分解中获得帽子矩阵（或其他衍生物）。从QR分解获得帽矩阵加权最小二乘回归

随着数据：

set.seed(0); xmat <- matrix(rnorm(500), nrow=50) ## model matrix 
y <- rowSums(rep(2:11,each=50)*xmat) ## arbitrary values to let `lm.wfit` work 
w <- runif(50, 1, 2) ## weights

的lwr()功能发生如下：

xmat2 <- w * xmat 
xx <- solve(crossprod(xmat, xmat2)) 
xmat1 <- tcrossprod(xx, xmat2) 
vmat <- tcrossprod(xmat1)

我需要的值，例如：

sum((xmat[1,] %*% xmat1)^2) 
sqrt(diag(vmat))

因为我使用reg <- lm.wfit(x=xmat, y=y, w=w)的那一刻但不能设法回到我看来是帽子矩阵（xmat1）。

来源

2013-12-13 Arthur

也许http://stackoverflow.com/questions/8065109/inverse-of-matrix-and-multiplication?rq=1是一个好的开始。 – Arthur

使用上面的链接，我可以做'xx < - qr（crossprod（xmat，k * xmat））; xxx < - qr.coef（xx，t（k * xmat [samp，]））; vmat < - tcrossprod（xxx）'给出了几乎相同的结果，但速度非常慢（约4倍）。有没有更有效的方法来做到这一点？有没有使用lm.wfit的方法？（对lm.wfit的一些处理是值得的，例如消除几乎冗余的列，以便微积分不会中止。） – Arthur

这个老问题是我刚刚回答的另一个老问题的延续：Compute projection/hat matrix via QR factorization, SVD (and Cholesky factorization?)。该答案讨论了用于计算普通最小二乘问题的帽子矩阵的3个选项，而这个问题是在加权最小二乘的情况下。但是，答案中的结果和方法将成为我在此答案的基础。具体来说，我只会演示QR方法。

OP提到的，我们可以用lm.wfit计算QR分解，但我们可以做到这一点使用qr.default自己，这是我将展示的方式。

我开始之前，我需要指出的是OP的代码没有做他的想法。xmat1不是帽子矩阵;相反，xmat %*% xmat1是。 vmat不是帽子矩阵，虽然我不知道它是什么。然后，我不明白这是什么意思：

sum((xmat[1,] %*% xmat1)^2) 
sqrt(diag(vmat))

第二个看起来像对角线帽子矩阵的，但正如我所说，vmat不是帽子矩阵。好吧，无论如何，我将继续进行帽子矩阵的正确计算，并展示如何获得其对角线和轨迹。

考虑一个玩具模型矩阵X和一些均匀，正权w：由行

set.seed(0); X <- matrix(rnorm(500), nrow = 50) 
w <- runif(50, 1, 2) ## weights must be positive 
rw <- sqrt(w) ## square root of weights

我们首先获得X1（X_tilde在胶乳段）重新缩放到X：

X1 <- rw * X

然后我们执行QR分解到X1。正如我在相关答案中所讨论的那样，我们可以使用或不使用列转轴来进行这种分解。 lm.fit或lm.wfit因此lm没有做旋转，但在这里我将使用旋转因子分解作为演示。

QR <- qr.default(X1, LAPACK = TRUE) 
Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))

报告中，我们并没有在计算tcrossprod(Q)走在链接的答案，因为这是普通最小二乘法。对于加权最小二乘，我们希望Q1和Q2：

Q1 <- (1/rw) * Q 
Q2 <- rw * Q

如果我们只希望对角线和帽子矩阵的痕迹，也没有必要做一个矩阵乘法先得到充分的帽子矩阵。我们可以使用

d <- rowSums(Q1 * Q2) ## diagonal 
# [1] 0.20597777 0.26700833 0.30503459 0.30633288 0.22246789 0.27171651 
# [7] 0.06649743 0.20170817 0.16522568 0.39758645 0.17464352 0.16496177 
#[13] 0.34872929 0.20523690 0.22071444 0.24328554 0.32374295 0.17190937 
#[19] 0.12124379 0.18590593 0.13227048 0.10935003 0.09495233 0.08295841 
#[25] 0.22041164 0.18057077 0.24191875 0.26059064 0.16263735 0.24078776 
#[31] 0.29575555 0.16053372 0.11833039 0.08597747 0.14431659 0.21979791 
#[37] 0.16392561 0.26856497 0.26675058 0.13254903 0.26514759 0.18343306 
#[43] 0.20267675 0.12329997 0.30294287 0.18650840 0.17514183 0.21875637 
#[49] 0.05702440 0.13218959 

edf <- sum(d) ## trace, sum of diagonals 
# [1] 10

在线性回归，d是每个数据的影响，并且它是用于（使用sqrt(1 - d)）产生的置信区间（使用sqrt(d)）和标准化残差是有用的。轨迹是模型的有效参数数量或有效自由度（因此我称之为edf）。我们看到edf = 10，因为我们已经使用了10个参数：X有10列，它不是等级不足的。

通常d和edf是我们所需要的。在极少数情况下，我们需要一个完整的帽子矩阵。为了得到它，我们需要一个昂贵的矩阵乘法：

H <- tcrossprod(Q1, Q2)

帽子矩阵是在帮助我们理解的模型是本地/稀疏的不是特别有用。我们绘制这个矩阵（关于如何在正确的方向绘制一个矩阵的细节和例子阅读?image）：

image(t(H)[ncol(H):1,])

我们看到，这个矩阵是完全致密。这意味着，每个数据的预测取决于所有的数据，即预测不是局部的。如果我们与其他非参数预测方法（如内核回归，黄土，P样条（惩罚B样条回归）和小波）进行比较，我们将观察稀疏的帽子矩阵。因此，这些方法被称为局部拟合。

来源

2016-09-04 05:12:43

从QR分解获得帽矩阵加权最小二乘回归

回答

相关问题