R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # 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 V2 V3 V4 V5 22.000 1.926000 667.000 124.000 56.48000 V2 1.926 0.288068 68.838 13.226 5.08791 V3 667.000 68.838000 21807.000 4042.000 1728.09000 V4 124.000 13.226000 4042.000 772.000 323.04000 V5 56.480 5.087910 1728.090 323.040 145.78240 > # > for ( k in(c(1:4,1:4))){ # sweep in 1:4, then out 1:4 + print(k) + yXXy <- regsweep(yXXy,k) + print(yXXy) + } [1] 1 V2 V3 V4 V5 0.04545455 0.08754545 30.31818 5.636364 2.5672727 V2 -0.08754545 0.11945545 10.44518 2.370364 0.1433427 V3 -30.31818182 10.44518182 1584.77273 282.545455 15.7190909 V4 -5.63636364 2.37036364 282.54545 73.090909 4.6981818 V5 -2.56727273 0.14334273 15.71909 4.698182 0.7828364 [1] 2 V2 V3 V4 V5 0.1096141 -0.7328711 22.663210 3.899193 2.4622210 V2 -0.7328711 8.3713214 87.439974 19.843076 1.1999680 V3 -22.6632096 -87.4399738 671.446303 75.280920 3.1852066 V4 -3.8991925 -19.8430758 75.280920 26.055604 1.8538212 V5 -2.4622210 -1.1999680 3.185207 1.853821 0.6108297 [1] 3 V2 V3 V4 V5 0.87456143 2.2184750 -0.033752825 1.3582488 2.354711259 V2 2.21847500 19.7583064 -0.130226309 10.0395194 0.785170339 V3 -0.03375283 -0.1302263 0.001489322 0.1121176 0.004743799 V4 -1.35824882 -10.0395194 -0.112117559 17.6152908 1.496703632 V5 -2.35471126 -0.7851703 -0.004743799 1.4967036 0.595719691 [1] 4 V2 V3 V4 V5 0.97929089 2.99258458 -0.025107862 -0.077106239 2.239306072 V2 2.99258458 25.48015062 -0.066326916 -0.569932082 -0.067849078 V3 -0.02510786 -0.06632692 0.002202927 -0.006364786 -0.004782399 V4 -0.07710624 -0.56993208 -0.006364786 0.056768861 0.084966161 V5 -2.23930607 0.06784908 0.004782399 -0.084966161 0.468550530 [1] 1 V2 V3 V4 V5 1.02114705 3.05586891 -0.025638819 -0.078736808 2.28666079 V2 -3.05586891 16.33520443 0.010399419 -0.334305524 -6.91087489 V3 0.02563882 0.01039942 0.001559191 -0.008341699 0.05263076 V4 0.07873681 -0.33430552 -0.008341699 0.050697762 0.26128197 V5 2.28666079 6.91087489 -0.052630764 -0.261281973 5.58908392 [1] 2 V2 V3 V4 V5 1.59281635 -0.1870725847 -0.0275842651 -0.016197410 3.57949602 V2 -0.18707258 0.0612174769 0.0006366262 -0.020465341 -0.42306632 V3 0.02758427 -0.0006366262 0.0015525702 -0.008128871 0.05703041 V4 0.01619741 0.0204653407 -0.0081288714 0.043856086 0.11984856 V5 3.57949602 -0.4230663241 -0.0570304079 -0.119848564 8.51284236 [1] 3 V2 V3 V4 V5 2.0829016 -0.19838342 17.7668394 -0.160621762 4.5927461 V2 -0.1983834 0.06147852 -0.4100466 -0.017132124 -0.4464515 V3 17.7668394 -0.41004663 644.0932642 -5.235751295 36.7329016 V4 0.1606218 0.01713212 5.2357513 0.001295337 0.4184456 V5 4.5927461 -0.44645145 36.7329016 -0.418445596 10.6077347 [1] 4 V2 V3 V4 V5 22.000 1.926000 667.000 124.000 56.48000 V2 1.926 0.288068 68.838 13.226 5.08791 V3 667.000 68.838000 21807.000 4042.000 1728.09000 V4 124.000 13.226000 4042.000 772.000 323.04000 V5 56.480 5.087910 1728.090 323.040 145.78240 > # last tableau should be (save rounding) same as first > # > # clean and close > rm(list=ls()) > q()