QR分解を用いたRの重回帰分析

Question

QR分解を用いた重回帰解を求める関数を書こうとしています。入力：yベクトルとX行列。出力：b、e、R^2。これまでのところ私はこれを持っており、ひどく詰まっています。それはすべてあなたがこの問題を解決するために使用することが許可されているどのくらいのR内蔵の施設のに依存

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) 

Answer 1

：私はあまりにも複雑なすべてのものを作っていると思います。私はすでにlmが許可されていないことを知っているので、ここには残りの話があります。あなたはその後、あなたは、単に解くQRベースの通常の最小二乗法のためlm.fit、.lm.fitまたはlsfitを使用することができますlm

よりも、他のルーチンを使用することを許可されている場合は

---

。これらの中でも

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

lm.fitとlsfitはかなり類似しているが、.lm.fitは、最も光秤量します。以下は、私たちが.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) 
    } 

あなたの仲間の同級生で：Toy R function for solving ordinary least squares by singular value decomposition、私はすでにSVDアプローチの正しさを確認するために、これを使用しています。あなたが使用することを許可されていない場合は

---

Rの内蔵.lm.fitが許可されていない場合qr.default

が、qr.defaultあるルーチンQR分解、それはまた、その複雑ではありません。

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) 
    } 

推定係数の分散共分散をさらに求めたい場合は、How to calculate variance of least squares estimator using QR decomposition in R?に従ってください。

---

あなたがさえその後、我々は、QRは、自分自身を分解記述する必要がqr.default

を使用することを許可されていない場合。 Writing a Householder QR factorization function in R codeがこれを与えています。我々はそれがシンQ要因であるにも関わらず、明示的Qを形成しているよう

が機能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を書くことができますが、非常に効率的ではありません。原理的にはXのQR分解と一緒にyを転記する必要があります。Qを形成する必要はありません。

---

、あなたの既存のコード

を修正したい場合これは、いくつかのデバッグの労力を必要とするので、いくつかの時間がかかるだろう。私はこれについて後でもう一度答えよう。