# sweep.r           Demonstration of sweep operator 
 #                     for regression computations 
 #                   
 #
 #  define function
   regsweep <- function(A,k) {
              d <- A[k,k]               # sweep out row k, col k
              A[k,] <- A[k,]/d
              b <- A[,k]
              b[k] <- 0                 # don't change row k here
              A <- A - outer(b,A[k,])   # main operation
              A[,k] <- -b/d             # fix col k
              A[k,k] <- 1/d             # diagonal element & done
              regsweep <- A     }
 # 
 # read data from Brownlee, page 463
 #
   all <- read.table("brown463.dat")
   Xy <- cbind(rep(1,22),as.matrix(all[,c(2,3,4,5)])) # X and y
 # inner product matrix
   yXXy <- t(Xy) %*% Xy
   yXXy                       # print it out
 #
   for ( k in(c(1:4,1:4))){   # sweep in 1:4, then out 1:4
     print(k)
     yXXy <- regsweep(yXXy,k)
     print(yXXy)
                          }
 # last tableau should be (save rounding) same as first
 #
 # clean and close
 rm(list=ls())
 q()