4

我正在尝试扩展lwr()package 的功能McSptial,它将加权回归拟合为非参数估计。在函数的核心中,它使用QR 分解而不是 QR 分解来lwr()反转矩阵,从而导致数值不稳定。solve()我想更改它,但无法弄清楚如何从 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) 的东西。

4

1 回答 1

4

这个老问题是我刚刚回答的另一个老问题的延续: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)像链接的答案那样继续计算,因为那是针对普通最小二乘法的。对于加权最小二乘,我们想要Q1Q2

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(d))和标准化残差(使用sqrt(1 - d))很有用。迹线是模型的有效参数数量或有效自由度(因此我称之为edf)。我们看到edf = 10,因为我们使用了 10 个参数:X有 10 列并且它不是秩不足的。

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

H <- tcrossprod(Q1, Q2)

Hat 矩阵在帮助我们理解模型是否是局部的/稀疏的方面特别有用。让我们绘制这个矩阵(阅读?image有关如何以正确方向绘制矩阵的详细信息和示例):

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

在此处输入图像描述

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

于 2016-09-04T05:12:43.847 回答