从QR分解中获得帽子矩阵以进行加权最小二乘回归 [英] Get hat matrix from QR decomposition for weighted least square regression
问题描述
我正在尝试扩展程序包McSptial的lwr()函数，该函数适合将weigthed回归作为非参数估计.在lwr()函数的核心中，它使用solve()而不是QR分解来反转矩阵，从而导致数值不稳定.我想更改它，但无法弄清楚以后如何从QR分解中获得帽子矩阵(或其他导数).
I am trying to extend the lwr() function of the package McSptial, which fits weigthed regressions as non-parametric estimation. In the core of the lwr() function, it inverts a matrix using solve() instead of a QR decomposition, resulting in numerical instability. I would like to change it but can't figure out how to get the hat matrix (or other derivatives) from the QR decomposition afterward.
有数据:
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).
For the moment I use reg <- lm.wfit(x=xmat, y=y, w=w) but cannot manage to get back what seems to me to be the hat matrix (xmat1) out of it.
推荐答案
This old question is a continuation of another old question I have just answered: Compute projection / hat matrix via QR factorization, SVD (and Cholesky factorization?). That answer discusses 3 options for computing hat matrix for an ordinary least square problem, while this question is under the context of weighted least squares. But result and method in that answer will be the basis of my answer here. Specifically, I will only demonstrate the QR approach.
OP提到我们可以使用lm.wfit来计算QR因式分解，但是我们可以自己使用qr.default来进行计算，这就是我将演示的方式.
OP mentioned that we can use lm.wfit to compute QR factorization, but we could do so using qr.default ourselves, which is the way I will show.
在继续之前，我需要指出OP的代码没有按照他的想法进行操作. xmat1不是帽子矩阵；而是xmat %*% xmat1是. vmat不是帽子矩阵，尽管我不知道它到底是什么.然后我不明白这些是什么:
Before I proceed, I need point out that OP's code is not doing what he thinks. xmat1 is not hat matrix; instead, xmat %*% xmat1 is. vmat is not hat matrix, although I don't know what it is really. Then I don't understand what these are:
sum((xmat[1,] %*% xmat1)^2)
sqrt(diag(vmat))
第二个看起来像帽子矩阵的对角线，但是正如我说的，vmat不是帽子矩阵.好吧，我将继续对帽子矩阵进行正确的计算，并说明如何获取其对角线和轨迹.
The second one looks like the diagonal of hat matrix, but as I said, vmat is not hat matrix. Well, anyway, I will proceed with the correct computation for hat matrix, and show how to get its diagonal and trace.
考虑一个玩具模型矩阵X和一些均匀的正权重w:
Consider a toy model matrix X and some uniform, positive weights 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
我们首先通过将行缩放为X来获得X1(在Latex段落中为X_tilde):
We first obtain X1 (X_tilde in the latex paragraph) by row rescaling to X:
X1 <- rw * X
然后我们将QR因式分解为X1.正如我在链接的答案中讨论的那样，无论有无列透视，我们都可以进行分解. lm.fit或lm.wfit因此，lm没有进行透视，但是在这里我将以透视分解为例.
Then we perform QR factorization to X1. As discussed in my linked answer, we can do this factorization with or without column pivoting. lm.fit or lm.wfit hence lm is not doing pivoting, but here I will use pivoted factorization as a demonstration.
QR <- qr.default(X1, LAPACK = TRUE)
Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
请注意，我们没有像链接答案中那样继续计算tcrossprod(Q)，因为这是针对普通最小二乘法的.对于加权最小二乘，我们需要Q1和Q2:
Note we did not go on computing tcrossprod(Q) as in the linked answer, because that is for ordinary least squares. For weighted least squares, we want Q1 and Q2:
Q1 <- (1 / rw) * Q
Q2 <- rw * Q
如果我们只想要帽子矩阵的对角线和迹线，则无需进行矩阵乘法即可首先获得完整的帽子矩阵.我们可以使用
If we only want the diagonal and trace of the hat matrix, there is no need to do a matrix multiplication to first get the full hat matrix. We can use
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(d))和标准化残差(使用sqrt(1 - d))很有用.迹线是模型的有效参数数量或有效自由度(因此我将其称为edf).我们看到了edf = 10，因为我们使用了10个参数:X有10列，并且不是秩不足的.
In linear regression, d is the influence of each datum, and it is useful for producing confidence interval (using sqrt(d)) and standardised residuals (using sqrt(1 - d)). The trace, is the effective number of parameters or effective degree of freedom for the model (hence I call it edf). We see that edf = 10, because we have used 10 parameters: X has 10 columns and it is not rank-deficient.
通常，我们只需要d和edf.在极少数情况下，我们需要完整的帽子矩阵.要获得它，我们需要昂贵的矩阵乘法:
Usually d and edf are all we need. In rare cases do we want a full hat matrix. To get it, we need an expensive matrix multiplication:
H <- tcrossprod(Q1, Q2)
帽子矩阵在帮助我们了解模型是否为局部/稀疏模型方面特别有用.让我们绘制该矩阵(有关如何以正确方向绘制矩阵的详细信息和示例，请阅读?image):
image(t(H)[ncol(H):1,])
我们看到这个矩阵是完全密集的.这意味着，每个数据的预测取决于所有数据，即预测不是局部的.虽然如果我们与其他非参数预测方法进行比较，例如核回归，黄土，P样条曲线(惩罚性B样条曲线回归)和小波，我们将观察到稀疏的帽子矩阵.因此，这些方法被称为局部拟合.
We see that this matrix is completely dense. This means, prediction at each datum depends on all data, i.e., prediction is not local. While if we compare with other non-parametric prediction methods, like kernel regression, loess, P-spline (penalized B-spline regression) and wavelet, we will observe a sparse hat matrix. Therefore, these methods are known as local fitting.
这篇关于从QR分解中获得帽子矩阵以进行加权最小二乘回归的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
