# 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
 # 
 # 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
 #
 # 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)
 #
 # compute coefficients and sse, etc.
   bhat <- backsolve(Xy[(1:4),(1:4)],Xy[(1:4),5])
   bhat
   sse <- sum(Xy[(5:22),5]^2)
   c(sse,sse/(12-4))                         # sse, variance est
   (sse/(12-4))*chol2inv(Xy[(1:4),(1:4)])    # est of Cov(bhat)
 #
 # 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
   sse <- sum(Xy[(5:22),5]^2)
   c(sse,sse/(22-4))                         # sse, variance est
   (sse/(22-4))*chol2inv(Xy[(1:4),(1:4)])    # est of Cov(bhat)
 # clean and close
 rm(list=ls())
 q()