>  #  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
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
i    1    1    1    1    2    2    2    2
j    1    2    1    2    1    2    1    2
k    2    2    3    3    2    2    3    3
>  #
>  # 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
     mu1 mu2 mu3 mu4
[1,]  11   5 5.0  13
[2,]  14   5 5.0  16
[3,]  14   5 5.0  16
[4,]  17   5 5.0  19
[5,]  11  14 3.5  13
[6,]  14  14 3.5  16
[7,]  14  14 3.5  16
[8,]  17  14 3.5  19
>  #
>  # 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
     one          
[1,]   1 1 0 1 1 0
[2,]   1 1 0 0 2 0
[3,]   1 1 0 1 0 1
[4,]   1 1 0 0 1 1
[5,]   1 0 1 1 1 0
[6,]   1 0 1 0 2 0
[7,]   1 0 1 1 0 1
[8,]   1 0 1 0 1 1
>  #
>  # 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
     [,1] [,2] [,3] [,4]
[1,]    2    2  2.0    4
[2,]    0    3  3.0    0
[3,]    0   12  1.5    0
[4,]    3    0  0.0    3
[5,]    6    0  0.0    6
[6,]    9    0  0.0    9
>   X %*% bs                    # same as mu?
     [,1] [,2] [,3] [,4]
[1,]   11    5  5.0   13
[2,]   14    5  5.0   16
[3,]   14    5  5.0   16
[4,]   17    5  5.0   19
[5,]   11   14  3.5   13
[6,]   14   14  3.5   16
[7,]   14   14  3.5   16
[8,]   17   14  3.5   19
>  #
>  # 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?
[1] 0
>   sum(abs( PXA %*% XA - XA )) # project onto right space?
[1] 0
>  #
>  # lambdaA will be zero if we get the same results
>   P18 %*% mu
     mu1 mu2  mu3 mu4
[1,]  14 9.5 4.25  16
[2,]  14 9.5 4.25  16
[3,]  14 9.5 4.25  16
[4,]  14 9.5 4.25  16
[5,]  14 9.5 4.25  16
[6,]  14 9.5 4.25  16
[7,]  14 9.5 4.25  16
[8,]  14 9.5 4.25  16
>   PXA %*% mu
     mu1 mu2 mu3 mu4
[1,]  14   5 5.0  16
[2,]  14   5 5.0  16
[3,]  14   5 5.0  16
[4,]  14   5 5.0  16
[5,]  14  14 3.5  16
[6,]  14  14 3.5  16
[7,]  14  14 3.5  16
[8,]  14  14 3.5  16
>  #
>  # 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?
[1] 2.775558e-16
>   sum(abs( PX %*% X - X ))    # project onto right space?
[1] 3.330669e-16
>   PXA %*% mu
     mu1 mu2 mu3 mu4
[1,]  14   5 5.0  16
[2,]  14   5 5.0  16
[3,]  14   5 5.0  16
[4,]  14   5 5.0  16
[5,]  14  14 3.5  16
[6,]  14  14 3.5  16
[7,]  14  14 3.5  16
[8,]  14  14 3.5  16
>   PX %*% mu
     mu1 mu2 mu3 mu4
[1,]  11   5 5.0  13
[2,]  14   5 5.0  16
[3,]  14   5 5.0  16
[4,]  17   5 5.0  19
[5,]  11  14 3.5  13
[6,]  14  14 3.5  16
[7,]  14  14 3.5  16
[8,]  17  14 3.5  19
>  #
>  # lambdaE will be zero if this is zero
>   sum(abs(mu - PX %*% mu))    # (I-Px)mu
[1] 1.776357e-15
>  # 
>   rm(list=ls())
>  q()