從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分解獲得帽矩陣加權最小二乘迴歸

回答

相關問題