#  q309.r     Question 2 on Third Quiz, 2009
 #
 # index vectors
  i <- c(rep(1,4),rep(2,4))
  j <- rep(c(1,2),4)
  k <- c(2,2,3,3,2,2,3,3)
  t(cbind(i,j,k))             # print out indices
 #
 # mean vectors
  mu1 <- 2 + 3*j + 3*k 
  mu2 <- 2 + 3*i*i
  mu3 <- 2 + 3/i
  mu4 <- 4 + 3*j + 3*k
  mu <- cbind(mu1,mu2,mu3,mu4)
  mu                          # print out mean vectors
 #
 # design matrices
  one <- rep(1,8)
  XA <- matrix(c(rep(1,4),rep(0,8)),8,2) # recycle for last 4
  XB <- matrix(c(1,0,1,0,1,0,1,0,1,2,0,1,1,2,0,1,0,0,1,1,0,0,1,1),8,3)
  X <- cbind(one,XA,XB)       # make design matrix
  X
 #
 # mean vectors as Xb
  bs <- matrix(c(2,0,0,3,6,9,2,3,12,0,0,0,2,3,1.5,0,0,0,4,0,0,3,6,9),6,4)
  bs
  X %*% bs                    # same as mu?
 #
 # now projection matrices
  P14 <- matrix(.25,4,4)      #Pone for four dimensions
  P18 <- matrix(.125,8,8)     #Pone for eight dimensions
  PXA <- matrix(0,8,8)
  PXA[1:4,1:4] <- P14
  PXA[5:8,5:8] <- P14         # PXA is block diagonal with Pone's
 #
 # check
  sum(abs( PXA %*% PXA - PXA ))# idempotent?
  sum(abs( PXA %*% XA - XA )) # project onto right space?
 #
 # lambdaA will be zero if we get the same results
  P18 %*% mu
  PXA %*% mu
 #
 # lambdaB will be zero if we get the same results
  Xe <- cbind(XA[,1],XB)      # basis for Column space of X
  PX <- Xe %*% solve( t(Xe)%*%Xe, t(Xe) ) 
  sum(abs( PX %*% PX - PX ))  # idempotent?
  sum(abs( PX %*% X - X ))    # project onto right space?
  PXA %*% mu
  PX %*% mu
 #
 # lambdaE will be zero if this is zero
  sum(abs(mu - PX %*% mu))    # (I-Px)mu
 # 
  rm(list=ls())
 q()