# FUNCTIONS FOR LINEAR REGRESSION USING HOUSEHOLDER TRANSFORMATIONS. Written # by Radford Neal, 2000, using the method described by R. Thisted in Elements # of Statistical Computing, Chapter 3. # # The functions intended for outside use are lr.house and, for specialized # purposes, decomp.house and back.solve. These functions check their arguments # for validity, and call "stop" with an appropriate message if a problem is # found. # LINEAR REGRESSION USING HOUSEHOLDER TRANSFORMATIONS. This function finds # the least-squares regression coefficients for a linear regression model, # and their standard errors. The arguments are as follows: # # X - matrix of predictor variables, n x p, where n is the number # of cases and p is the number of predictor variables. # # y - vector of responses, of length n. # # The value returned is a list with three components: coef = the vector of # estimates for the p regression coefficients, res.std.dev = the estimated # residual standard deviation, and std.errs = the vector of standard errors # for the p regression coefficients. lr.house <- function(X,y) { # Check arguments for validity. if (!is.matrix(X) || !is.vector(y) || dim(X)[1]!=length(y) || dim(X)[1]1) { u[1:(i-1)] <- 0 } if (x[i]>0) { s <- -s } u[i] <- x[i] - s u <- u / sqrt(s^2 - x[i]*s) # Return the "u" value that defines the transformation, along with the # "s" value computed along the way. list (u = u, s = s) } # APPLY A HOUSEHOLDER TRANSFORMATION TO A VECTOR. The transformation is # represented by a list that contains a "u" component which is the vector on # which the transformation is based, as is returned by make.house. The # "u" vector is assumed to have been scaled to have a squared magnitude of 2. apply.house <- function(v,tr) { f <- as.vector(tr$u %*% v) v - f * tr$u } # SOLVE A SYSTEM OF EQUATIONS BY BACK-SUBSTITUTION. The value returned # is vector that is the solution of Xb=y, where X is a square upper-triangular # matrix, and y is a vector of matching length. # # Calls "stop" with an error message if X is singular. back.solve <- function(X,y) { if (!is.matrix(X) || !is.vector(y) || dim(X)[1]!=length(y) || dim(X)[2]!=length(y)) { stop("Bad arguments to back.solve") } p <- length(y) b <- rep(0,p) # Find the last component of the solution. if (X[p,p]==0) { stop("Singular matrix in back.solve") } b[p] <- y[p] / X[p,p] if (p>1) { # Find components of the solution before the last in reverse order. for (j in (p-1):1) { if (X[j,j]==0) { stop("Singular matrix in back.solve") } b[j] <- (y[j] - b[(j+1):p] %*% X[j,(j+1):p]) / X[j,j] } } b }