Mehrere Regressionsanalyse in R mit QR-Zerlegung

Question

Ich versuche, eine Funktion zum Lösen multipler Regression mit QR-Dekomposition zu schreiben. Eingabe: y Vektor und X Matrix; Ausgabe: b, e, R^2. Bis jetzt habe ich das und bin furchtbar fest; Ich glaube, ich habe alles gemacht habe, viel zu kompliziert:

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

Es hängt alles davon ab, wie viel von R-interner Einrichtung, die Sie dieses Problem lösen benutzen dürfen. Ich weiß bereits, dass lm nicht erlaubt ist, also hier ist der Rest der Geschichte.

---

Wenn Sie irgendwelche anderen Routinen verwenden dürfen als lm

Dann können Sie einfach lm.fit, .lm.fit oder lsfit für QR-basierte gewöhnlichen kleinsten Quadrate verwenden zu lösen.

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

Unter denen, .lm.fit ist das Licht-gewogen, während lm.fit und lsfit ziemlich ähnlich sind. Hier finden Sie, was wir können über .lm.fit tun:

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

In der Frage, indem Sie Ihre Mitschüler: Toy R function for solving ordinary least squares by singular value decomposition, ich habe dies bereits verwendet, um die Richtigkeit der SVD-Ansatz zu überprüfen.

---

Wenn Sie nicht eingebaut, um zu verwenden R erlaubt QR-Faktorisierung Routine qr.default

Wenn .lm.fit nicht erlaubt ist, aber qr.default ist, dann ist es auch nicht so kompliziert.

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

Wenn Sie weitere Varianz-Kovarianz der geschätzten Koeffizienten möchten, How to calculate variance of least squares estimator using QR decomposition in R?.

---

Wenn Sie nicht einmal qr.default

Dann müssen wir uns QR Zersetzung schreiben benutzen. Writing a Householder QR factorization function in R code gibt dies.

die Funktion myqr dort verwenden, können wir

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 nicht sehr effizient ist, schreiben wir Q ausdrücklich gebildet haben, auch wenn es der Dünn- Q Faktor ist. Im Prinzip sollten wir y zusammen mit der QR-Faktorisierung von X drehen, also Q braucht nicht gebildet werden.

---

Wenn Sie Ihre vorhandenen Code

Dies erfordert einige Debug-Aufwand zu beheben, so einige Zeit dauern würde.Ich werde später noch eine Antwort darauf geben.