8

QR 分解を使用して重回帰を解くための関数を作成しようとしています。入力: y ベクトルと X 行列。出力: b、e、R^2。これまでのところ、私はこれを手に入れましたが、ひどく立ち往生しています。私はすべてを複雑にしすぎたと思います：

QR.regression <- function(y, X) {
X <- as.matrix(X)
y <- as.vector(y)
p <- as.integer(ncol(X))
if (is.na(p)) stop("ncol(X) is invalid")
n <- as.integer(nrow(X))
if (is.na(n)) stop("nrow(X) is invalid")
nr <- length(y)
nc <- NCOL(X)

# Householder 
for (j in seq_len(nc)) {
id <- seq.int(j, nr)
sigma <- sum(X[id, j]^2)
s <- sqrt(sigma)
diag_ej <- X[j, j]
gamma <- 1.0 / (sigma + abs(s * diag_ej))
kappa <- if (diag_ej < 0) s else -s
X[j,j] <- X[j, j] - kappa
if (j < nc)
for (k in seq.int(j+1, nc)) {
yPrime <- sum(X[id,j] * X[id,k]) * gamma
X[id,k] <- X[id,k] - X[id,j] * yPrime
}
yPrime <- sum(X[id,j] * y[id]) * gamma
y[id] <- y[id] - X[id,j] * yPrime
X[j,j] <- kappa
} # end of Householder transformation

rss <- sum(y[seq.int(nc+1, nr)]^2)  # residuals sum of squares
e <- rss/nr
e <- mean(residuals(QR.regression)^2)
beta <- solve(t(X) %*% X, t(X) %*% y)
for (i in seq_len(ncol(X))) # set zeros in the lower triangular side of X
X[seq.int(i+1, nr),i] <- 0
Rsq <- (X[1:nc,1:nc])^2
return(list(Rsq=Rsq, y = y, beta = beta, e = e))
}


UPDATE:
my.QR <- function(y, X) {
X <- as.matrix(X)
y <- as.vector(y)
p <- as.integer(ncol(X))
if (is.na(p)) stop("ncol(X) is invalid")
n <- as.integer(nrow(X))
if (is.na(n)) stop("nrow(X) is invalid")
qr.X <- qr(X)
b <- solve(t(X) %*% X, t(X) %*% y)
e <- as.vector(y - X %*% beta) #e
R2 <- (X[1:p, 1:p])^2
return(list(b = b, e= e, R2 = R2 ))
}

X <- matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3)
y <- c(1,2,3,4)
my.QR(X, y)
4

1 に答える 1

12

それはすべて、この問題を解決するために使用できる R の組み込み機能の量に依存します。私はそれが許可されていないことをすでに知っているlmので、ここに残りの話があります.

以外のルーチンの使用が許可されている場合lm

次にlm.fit、.lm.fitまたはlsfitを使用して、QR ベースの通常の最小二乗法を解くことができます。

lm.fit(X, y)
.lm.fit(X, y)
lsfit(X, y, intercept = FALSE)

それらの中で、.lm.fitは最も軽量ですが、lm.fitとlsfitはかなり似ています。経由でできることは次のとおりです.lm.fit。

f1 <- function (X, y) {
  z <- .lm.fit(X, y)
  RSS <- crossprod(z$residuals)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  list(coefficients = z$coefficients, residuals = z$residuals, R2 = R2)
  }

同級生の質問:特異値分解による通常の最小二乗法を解くためのおもちゃの R 関数 で、SVD アプローチの正しさを確認するために既にこれを使用しました。

R の組み込み QR 因数分解ルーチンの使用が許可されていない場合qr.default

.lm.fitが許可されていない場合qr.defaultでも、それほど複雑ではありません。

f2 <- function (X, y) {
  ## QR factorization `X = QR`
  QR <- qr.default(X)
  ## After rotation of `X` and `y`, solve upper triangular system `Rb = Q'y` 
  b <- backsolve(QR$qr, qr.qty(QR, y))
  ## residuals
  e <- as.numeric(y - X %*% b)
  ## R-squared
  RSS <- crossprod(e)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  ## multiple return
  list(coefficients = b, residuals = e, R2 = R2)
  }

さらに推定係数の分散共分散が必要な場合は、R で QR 分解を使用して最小二乗推定量の分散を計算する方法に従ってください。.

使用すら許可されていない場合qr.default

次に、QR 分解を自分で作成する必要があります。RコードでHouseholder QR因数分解関数を書くと、これが得られます。

myqrそこの関数を使用して、次のように書くことができます

f3 <- function (X, y) {
  ## our own QR factorization
  ## complete Q factor is not required
  QR <- myqr(X, complete = FALSE)
  Q <- QR$Q
  R <- QR$R
  ## rotation of `y`
  Qty <- as.numeric(crossprod(Q, y))
  ## solving upper triangular system
  b <- backsolve(R, Qty)
  ## residuals
  e <- as.numeric(y - X %*% b)
  ## R-squared
  RSS <- crossprod(e)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  ## multiple return
  list(coefficients = b, residuals = e, R2 = R2)
  }

f3は、シンファクターですQが、明示的に形成したように、非常に効率的ではありません。Q原則として、yの QR 因数分解に沿って回転する必要があるためX、Q形成する必要はありません。

既存のコードを修正したい場合

これにはデバッグ作業が必要なため、時間がかかります。これについては、後で別の回答をします。

于 2016-10-25T23:39:10.337 に答える