6

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)

ご協力いただきありがとうございます！

4

1 に答える 1

10

セットアップ (コードにコピー)

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」係数の形成をスキップ)

R <- qr.R(qr.X)forを抽出する必要はありませんbacksolve。を使用する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))
  }
于 2016-09-25T14:29:27.350 に答える