如何使用 R 中的 QR 分解计算最小二乘估计量的方差？
How to calculate variance of least squares estimator using QR decomposition in R?
我正在尝试学习 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)
感谢您的帮助！
设置（复制您的代码）
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] ## can be optimized; see Remark 1 below
R <- qr.R(qr.X) ## can be optimized; see Remark 2 below
beta <- as.vector(backsolve(R, b))
res <- as.vector(y - X %*% beta)
数学
计算
剩余自由度为n - p，因此估计方差为
se2 <- sum(res ^ 2) / (n - p)
因此，估计系数的方差协方差矩阵为
V <- chol2inv(R) * se2
# [,1] [,2]
#[1,] 0.22934170 -0.07352916
#[2,] -0.07352916 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（跳过形成 'Q' 因子）
代替b <- (t(qr.Q(qr.X)) %*% y)[1:p]，你可以使用函数qr.qty（避免形成'Q'矩阵）：
b <- qr.qty(qr.X, y)[1:p]
备注 2（跳过形成 'R' 因子）
您不必为 backsolve 提取 R <- qr.R(qr.X)；使用 qr.X$qr 就足够了：
beta <- as.vector(backsolve(qr.X$qr, b))
附录：估计函数
以上是最简单的演示。在实践中，需要处理列旋转和排名不足的问题。下面是一个实现。 X 是模型矩阵，y 是响应。结果应与 lm(y ~ X + 0).
进行比较
qr_estimation <- function (X, y) {
## QR factorization
QR <- qr(X)
r <- QR$rank
piv <- QR$pivot[1:r]
## estimate identifiable coefficients
b <- qr.qty(QR, y)[1:r]
beta <- backsolve(QR$qr, b, r)
## fitted values
yhat <- base::c(X[, piv] %*% beta)
## residuals
resi <- y - yhat
## error variance
se2 <- base::c(crossprod(resi)) / (nrow(X) - r)
## variance-covariance for coefficients
V <- chol2inv(QR$qr, r) * se2
## post-processing on pivoting and rank-deficiency
p <- ncol(X)
beta_full <- rep.int(NA_real_, p)
beta_full[piv] <- beta
V_full <- matrix(NA_real_, p, p)
V_full[piv, piv] <- V
## return
list(coefficients = beta_full, vcov = V_full,
fitted.values = yhat, residuals = resi, sig = sqrt(se2))
}
我正在尝试学习 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)
感谢您的帮助！
设置（复制您的代码）
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] ## can be optimized; see Remark 1 below
R <- qr.R(qr.X) ## can be optimized; see Remark 2 below
beta <- as.vector(backsolve(R, b))
res <- as.vector(y - X %*% beta)
数学
计算
剩余自由度为n - p，因此估计方差为
se2 <- sum(res ^ 2) / (n - p)
因此，估计系数的方差协方差矩阵为
V <- chol2inv(R) * se2
# [,1] [,2]
#[1,] 0.22934170 -0.07352916
#[2,] -0.07352916 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（跳过形成 'Q' 因子）
代替b <- (t(qr.Q(qr.X)) %*% y)[1:p]，你可以使用函数qr.qty（避免形成'Q'矩阵）：
b <- qr.qty(qr.X, y)[1:p]
备注 2（跳过形成 'R' 因子）
您不必为 backsolve 提取 R <- qr.R(qr.X)；使用 qr.X$qr 就足够了：
beta <- as.vector(backsolve(qr.X$qr, b))
附录：估计函数
以上是最简单的演示。在实践中，需要处理列旋转和排名不足的问题。下面是一个实现。 X 是模型矩阵，y 是响应。结果应与 lm(y ~ X + 0).
qr_estimation <- function (X, y) {
## QR factorization
QR <- qr(X)
r <- QR$rank
piv <- QR$pivot[1:r]
## estimate identifiable coefficients
b <- qr.qty(QR, y)[1:r]
beta <- backsolve(QR$qr, b, r)
## fitted values
yhat <- base::c(X[, piv] %*% beta)
## residuals
resi <- y - yhat
## error variance
se2 <- base::c(crossprod(resi)) / (nrow(X) - r)
## variance-covariance for coefficients
V <- chol2inv(QR$qr, r) * se2
## post-processing on pivoting and rank-deficiency
p <- ncol(X)
beta_full <- rep.int(NA_real_, p)
beta_full[piv] <- beta
V_full <- matrix(NA_real_, p, p)
V_full[piv, piv] <- V
## return
list(coefficients = beta_full, vcov = V_full,
fitted.values = yhat, residuals = resi, sig = sqrt(se2))
}