QR 분해를 이용한 R의 다중 회귀 분석

Question

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) 

Answer 1

그것은 모두가이 문제를 해결하기 위해 사용하는 것이 허용되는 R에 내장 된 시설의 대부분에 따라 달라집니다. 나는 이미 lm이 허용되지 않는다는 것을 알고 있으므로 여기에 나머지 이야기가 있습니다. 당신이 그럼 당신은 간단하게 해결 QR 기반의 정규 방정식에 대한 lm.fit, .lm.fit 또는 lsfit을 사용할 수 있습니다 lm

보다 다른 루틴을 사용할 수 있도록 허용하는 경우

---

. 그 중

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

, .lm.fit는 가장 가벼운 무게, lm.fit 및 lsfit 꽤 비슷있다. 질문에

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

을 동료 급우에 의해 : 다음은 우리가 .lm.fit를 통해 무엇을 할 수 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이이를 제공합니다. 기능이 myqr를 사용

, 우리는 그것이 얇은 - Q 요인에도 불구하고, 우리는 명시 적으로 Q을 형성했다으로

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을 형성 할 필요가 없습니다.

---

기존 코드

이를 해결하려면

시간이 좀 걸릴 것 때문에 디버깅의 노력이 필요합니다.나는 나중에 이것에 대한 또 다른 대답을 할 것이다.