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. > # flyreg.r Demonstration of Givens transformations > # for regression computations on the fly > # > all <- read.table("brown463.dat") > Xy <- cbind(rep(1,22),as.matrix(all[,c(2,3,4,5)])) # X and y > Xy V2 V3 V4 V5 [1,] 1 0.124 33 8 2.81 [2,] 1 0.049 31 6 2.55 [3,] 1 0.181 38 8 2.80 [4,] 1 0.004 17 2 2.24 [5,] 1 0.022 20 4 2.78 [6,] 1 0.152 39 6 2.52 [7,] 1 0.075 30 7 2.88 [8,] 1 0.054 29 7 2.45 [9,] 1 0.043 35 6 2.50 [10,] 1 0.041 31 5 2.69 [11,] 1 0.017 23 4 2.66 [12,] 1 0.022 21 3 2.45 [13,] 1 0.016 8 3 2.24 [14,] 1 0.010 23 3 2.43 [15,] 1 0.063 37 6 2.38 [16,] 1 0.170 40 8 2.72 [17,] 1 0.125 38 6 2.41 [18,] 1 0.015 25 4 2.38 [19,] 1 0.221 39 7 2.52 [20,] 1 0.171 33 7 2.52 [21,] 1 0.097 38 6 2.66 [22,] 1 0.254 39 8 2.89 > # > # function to do Givens on col i, rows i,j > Givens <- function(X,i,j) { + a <- X[i,i] ; b <- X[j,i] + r <- sqrt(a*a+b*b) # length + Uij <- matrix( c(a/r,-b/r,b/r,a/r),2,2) # essential + X[c(i,j),] <- Uij %*% X[c(i,j),] # main computation + Givens <- X } # return modified X > # > # go through a sequence of steps, putting zeros in X > Xy <- Givens(Xy,1,2) > Xy <- Givens(Xy,1,3) > Xy <- Givens(Xy,2,3) > Xy <- Givens(Xy,1,4) > Xy <- Givens(Xy,2,4) > Xy <- Givens(Xy,3,4) > Xy # U34*U24*U14*U23*U13*U12*Xy V2 V3 V4 V5 [1,] 2.000000e+00 1.790000e-01 59.500000 12.000000 5.2000000 [2,] 8.055775e-17 1.360625e-01 14.011945 4.365641 0.4288471 [3,] 7.423398e-17 1.348498e-17 6.812886 1.883031 0.1249131 [4,] 1.804765e-17 3.278448e-18 0.000000 1.181258 0.1291778 [5,] 1.000000e+00 2.200000e-02 20.000000 4.000000 2.7800000 [6,] 1.000000e+00 1.520000e-01 39.000000 6.000000 2.5200000 [7,] 1.000000e+00 7.500000e-02 30.000000 7.000000 2.8800000 [8,] 1.000000e+00 5.400000e-02 29.000000 7.000000 2.4500000 [9,] 1.000000e+00 4.300000e-02 35.000000 6.000000 2.5000000 [10,] 1.000000e+00 4.100000e-02 31.000000 5.000000 2.6900000 [11,] 1.000000e+00 1.700000e-02 23.000000 4.000000 2.6600000 [12,] 1.000000e+00 2.200000e-02 21.000000 3.000000 2.4500000 [13,] 1.000000e+00 1.600000e-02 8.000000 3.000000 2.2400000 [14,] 1.000000e+00 1.000000e-02 23.000000 3.000000 2.4300000 [15,] 1.000000e+00 6.300000e-02 37.000000 6.000000 2.3800000 [16,] 1.000000e+00 1.700000e-01 40.000000 8.000000 2.7200000 [17,] 1.000000e+00 1.250000e-01 38.000000 6.000000 2.4100000 [18,] 1.000000e+00 1.500000e-02 25.000000 4.000000 2.3800000 [19,] 1.000000e+00 2.210000e-01 39.000000 7.000000 2.5200000 [20,] 1.000000e+00 1.710000e-01 33.000000 7.000000 2.5200000 [21,] 1.000000e+00 9.700000e-02 38.000000 6.000000 2.6600000 [22,] 1.000000e+00 2.540000e-01 39.000000 8.000000 2.8900000 > # > # now recompute with each row/obs > for (i in(5:12)) { # stop partway & estimate + for ( j in(1:4)) { Xy <- Givens(Xy,j,i)} + print(Xy[i,]) } # zero and z2's (residuals) V2 V3 V4 V5 3.701933e-17 7.534715e-18 0.000000e+00 -1.110223e-16 2.252119e-01 V2 V3 V4 V5 1.041864e-17 5.302713e-18 0.000000e+00 0.000000e+00 1.187525e-01 V2 V3 V4 V5 -3.407617e-17 -6.883041e-18 1.985484e-16 0.000000e+00 8.383159e-02 V2 V3 V4 V5 -2.620147e-17 -1.874070e-18 -3.954203e-17 0.000000e+00 -2.987594e-01 V2 V3 V4 V5 -1.177732e-17 -5.868795e-18 -8.286824e-16 0.000000e+00 8.718340e-02 V2 V3 V4 V5 1.085066e-17 5.734126e-18 -3.537413e-16 0.000000e+00 2.186484e-01 V2 V3 V4 V5 -7.623241e-17 1.504843e-18 3.606311e-17 0.000000e+00 1.321907e-01 V2 V3 V4 V5 -4.920840e-17 3.567798e-18 1.118436e-16 0.000000e+00 -3.699531e-02 > # > # compute coefficients and sse, etc. > bhat <- backsolve(Xy[(1:4),(1:4)],Xy[(1:4),5]) > bhat [1] 2.52612573 0.77548062 -0.01527082 0.08647713 > sse <- sum(Xy[(5:22),5]^2) > c(sse,sse/(12-4)) # sse, variance est [1] 63.821659 7.977707 > (sse/(12-4))*chol2inv(Xy[(1:4),(1:4)]) # est of Cov(bhat) [,1] [,2] [,3] [,4] [1,] 22.8908629 85.605266 -0.94191140 -0.1058170 [2,] 85.6052661 740.763868 -3.48091502 -6.0627963 [3,] -0.9419114 -3.480915 0.06152643 -0.1108742 [4,] -0.1058170 -6.062796 -0.11087422 0.6741877 > # > # now add the remaining observations > for (i in(13:22)) { # to the end + for ( j in(1:4)) { Xy <- Givens(Xy,j,i)} } > # > # compute coefficients and sse, etc. > bhat <- backsolve(Xy[(1:4),(1:4)],Xy[(1:4),5]) > bhat [1] 2.239306072 -0.067849078 -0.004782399 0.084966161 > sse <- sum(Xy[(5:22),5]^2) > c(sse,sse/(22-4)) # sse, variance est [1] 0.46855053 0.02603059 > (sse/(22-4))*chol2inv(Xy[(1:4),(1:4)]) # est of Cov(bhat) [,1] [,2] [,3] [,4] [1,] 0.0254915147 0.077898727 -6.535723e-04 -0.0020071205 [2,] 0.0778987273 0.663263227 -1.726528e-03 -0.0148356655 [3,] -0.0006535723 -0.001726528 5.734347e-05 -0.0001656791 [4,] -0.0020071205 -0.014835665 -1.656791e-04 0.0014777267 > # clean and close > rm(list=ls()) > q()