如何使用R中的QR分解計算最小二乘估計量的方差？

我試圖學習QR分解，但無法弄清楚如何獲得beta_hat的方差而不訴諸於傳統的矩陣計算。我與iris數據集練習，並在這裏是我到目前爲止有：如何使用R中的QR分解計算最小二乘估計量的方差？

y<-(iris$Sepal.Length) 
x<-(iris$Sepal.Width) 
X<-cbind(1,x) 
n<-nrow(X) 
p<-ncol(X) 
qr.X<-qr(X) 
b<-(t(qr.Q(qr.X)) %*% y)[1:p] 
R<-qr.R(qr.X) 
beta<-as.vector(backsolve(R,b)) 
res<-as.vector(y-X %*% beta)

感謝您的幫助！

來源

2016-09-19 Nudibranch

@ZheyuanLi接受的答案主要是關於表明一個答案工作，大約給人一種有點少獎勵回答者投入的努力。這不一定是關於職業行爲;當問題沒有（完全）解決時，不接受答案也可以歸類爲專業行爲。 – Jaap

@ZheyuanLi你好，我一直在路上，無法接受移動....謝謝！ – Nudibranch

自由的殘餘程度n - p（更準確地說，n - qr.X$rank。在下文中，我將使用qr.X$rank代替p。），所以估計方差是

因此，您估計的beta_hat的方差協方差矩陣只是

chol2inv(R) * se2 

#   [,1]   [,2] 
#[1,] 0.22934170 -0.07352916 
#[2,] -0.07352916 0.02405009

很多時候，我們只想要邊際方差，即這個完全方差 - 協方差矩陣的對角線。在這種情況下，不需要呼叫chol2inv首先明確地形成完全協方差矩陣。我們可以在低得多的計算成本，直接計算對角線：

Rinv <- backsolve(R, diag(qr.X$rank)) 
rowSums(Rinv^2) * se2 
# [1] 0.22934170 0.02405009

讓我們用lm比較檢查正確性：

fit <- lm(Sepal.Length ~ Sepal.Width, iris) 

vcov(fit) 

#   (Intercept) Sepal.Width 
#(Intercept) 0.22934170 -0.07352916 
#Sepal.Width -0.07352916 0.02405009

是的，正確的！

最後，在你的代碼貼一些意見：

而不是

b <- (t(qr.Q(qr.X)) %*% y)[1:p]

你可以使用

b <- qr.qty(qr.X, y)[1:qr.X$rank]

您不必提取R <- qr.R(qr.X) for backsolve;使用qr.X$qr充足：

beta <- as.vector(backsolve(qr.X$qr,b)) 
Rinv <- backsolve(qr.X$qr, diag(qr.X$rank))

總之，這裏是一個完整的會話：

## data 
y <- iris$Sepal.Length 
x <- iris$Sepal.Width 

## add intercept to model matrix 
X <- cbind(1,x) 

## QR factorization of model matrix 
qr.X <- qr(X) 

## get estimated coefficient 
b <- qr.qty(qr.X, y) 
beta <- as.vector(backsolve(qr.X$qr, b)) 

## compute residuals (I don't use `qr.resid` and `qr.fitted`, though I can) 
res <- as.vector(y - X %*% beta) 

## residual standard error 
se2 <- sum(res^2)/(nrow(X) - qr.X$rank) 

## full variance-covariance matrix 
chol2inv(qr.X$qr) * se2

來源

2016-09-25 14:29:27

如何使用R中的QR分解計算最小二乘估計量的方差？

回答

相關問題