############# HOMEWORK 5, STAT 541, DUE FRI, OCT 23, 2009 ###############

# This homework is an exercise in Least Squares computation.  We will
# build up R code that
# 1) computes a Q-R decomposition of a matrix X,
# 2) solves equations when the coefficients are upper triangular,
# 3) inverts an upper triangular matrix R, and finally
# 4) computes and returns most things we find useful in regression.
# We will be negligent in the sense that we will ignore precautions 
# that a numerical analyst would build into his code, and we will not
# consider measures to improve numerical conditioning such as pivoting
# (switching the order of the variables) either.

# BACKGROUND: Most LS solvers are based on decompositions of the
# design matrix of the form X = Q R, where R is square and triangular
# and Q is of the same size as X but has orthonormal columns.
# Note that here the columns of Q form an orthonormal basis
# of the column space of X, not residual space.  The matrix
#      P = Q Q^T
# is then the orthogonal projection onto the column space of X.
# With this knowledge write 'yhat' in two ways, where 'b' as usual
# is the vector of LS estimates:
#      yhat   =   X b = Q Q^T y
#               Q R b = Q Q^T y   (left-multiply both sides with Q^T)
#                 R b =   Q^T y
# The last equation is easily solved when  R  is triangular.
# Furthermore, for inference we want the matrix (X^T X)^{-1},
# but with a Q-R decomposition this is easily gotten by
# inverting R because
#   (X^T X)^{-1}  =  R^{-1} R^{-1}^T


#----------------------------------------------------------------

# PROBLEM 1: Write a function 'gs(X)' to implement the so-called
# modified Gram-Schmidt procedure for orthonormalizing the columns
# of a matrix X.  As this procedure orthonormalizes the columns
# of X to form the matrix Q of the same size, it stores the
# coefficients to reconstitute X from Q in an upper triangular
# matrix R.

# Here is the algorithm:
#   Initialize Q with a copy of X
#   Loop over the columns of Q (j=1,2,...,p-1) and do the following:
#     At stage j,
#       normalize Q[,j] in place
#       adjust (=orthogonalize) the columns of Q[,(j+1):p] w.r.t. Q[,j]
# As this algorithm proceeds, store suitable norms and inner products
# in the appropriate cells of R so as to reconstitute X from Q.
# If you do this right, it should hold at the end:  X = Q R.
# Finish the function gs() with  return(list(Q=Q, R=R))

# In this algorithm, use for-loops and sum() to calculate inner products
# and sqrt() for lengths, but not lm() or any other high-level function.

# Test your function on these data:
   X <- cbind(rep(1,10), 1:10, (1:10)^2)
   sol <- gs(X)
# Show the following and comment on what they mean and what they should show:
   sol$Q %*% sol$R
   t(sol$Q) %*% sol$Q

# YOUR SOLUTION:


#----------------------------------------------------------------

# PROBLEM 2: Write a function 'tri.solve(R,z)' that accepts an upper
# triangular matrix R and a column z and returns the solution b of 
# z = R b.  Show the results for these inputs:

   z <- 10:1
   R <- (row(diag(10)) <= col(diag(10)))

# You may use loops, but you may also use the sum() function inside 
# tri.solve() to spare yourself an inner loop.  You must not use 
# canned solvers such as solve() or ginv() (the latter in package MASS).

# YOUR SOLUTION:


#----------------------------------------------------------------

# PROBLEM 3: No programming, just thinking.
# Assume R and R1 are upper triangular square matrices.
# Questions:
# a) If the vector b has non-zero values only in the first k entries,
#    what can one say about z=Rb?
# b) When is R invertible?
# c) If the vector z has non-zero values only in the first k entries,
#    and if R is invertible, what can one say about the vector b?
# d) What can one say about R1%*%R?
# e) What can one say about R^{-1}?

# YOUR SOLUTION:


#----------------------------------------------------------------

# PROBLEM 4: For inference about the regression coefficients, we need
# the matrix (X^T X)^{-1}.  In preparation for it, we need the
# inversion of an upper triangular matrix R because from a Q-R
# decomposition X=QR we can easily obtain (X^T X)^{-1} if we have
# R^{-1}.  Package the inversion in a function inv(R).  Try it out on
# the R matrix of the above 'sol':

  inv(sol$R)
  inv(sol$R) %*% sol$R
# Show both and comment on both.

# YOUR SOLUTION:


#----------------------------------------------------------------

# PROBLEM 5: Using tri.solve(R,z), inv(R) and gs(X,y), write a
# function reg(X,y) that computes a list with the following named
# elements:

# a) "Coeffs": a matrix with one row per regression coefficient and columns
#    named "Coefficient", "Std.Err.Est", "t-Statistic", "P-Value"
# b) "Var": the estimated coeff variance/covariance matrix  s^2*(X^T X)^{-1}
# c) "RMSE": a number, s
# d) "R2": a number, R2
# e) "Fstat": the value of the overall F statistic
# f) "Fpval": the corresponding p-value
# g) "Residuals": the vector of residuals

# Assume that X does not have a column of 1's but an intercept is desired,
# hence first thing in the function do  X <- cbind(Icept=1,X) .

# You may use matrix multiplication if convenient, but no high-level
# functions other than things like sqrt() and pt().

# Try your solution on the following data and show the full result:
  X <- cbind(Pred1=1:10, Pred2=(1:10)^2)
  y <- rep(1,10) + 2*X[,1] + 3*X[,2] + resid(lm(rep(0:1,5)~X))
  reg(X,y)
# Compare with
  summary(lm(y~X))
# No need to comment, but you need to get the numbers right.

# YOUR SOLUTION:


#----------------------------------------------------------------