Как вычислить дисперсию оценки наименьших квадратов с использованием QR-декомпозиции в R?

Question

Я пытаюсь изучить QR-декомпозицию, но не могу понять, как получить дисперсию beta_hat, не прибегая к традиционным вычислениям матрицы. Я тренировался с набором iris данных, и вот что у меня до сих пор:

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) 

Спасибо за вашу помощь!

Answer 1

Остаточная степень свободы 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 

Да, все правильно!

---

Наконец, некоторые комментарии на ваш публикуемую код:

1. Вместо

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

   вы можете просто использовать

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

2. Вы не должны экстракт R <- qr.R(qr.X) для 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