pershop <- read.csv("http://www.stat.ufl.edu/~winner/sta4702/personality_shopping.csv") attach(pershop) R <- as.matrix(pershop[,3:11]) R R11 <- R[6:9,6:9] R22 <- R[1:5,1:5] R12 <- R[6:9,1:5] R21 <- R[1:5,6:9] p <- nrow(R11); q <- nrow(R22) R.CC <- rbind(cbind(R11,R12),cbind(R21,R22)) R11.lam <- eigen(R11)$val R11.e <- eigen(R11)$vec R11.m12 <- matrix(rep(0,p^2),ncol=p) R11.m1 <- matrix(rep(0,p^2),ncol=p) for (i1 in 1:p) { R11.m12 <- R11.m12 + (1/sqrt(R11.lam[i1])) * (R11.e[,i1] %*% t(R11.e[,i1])) R11.m1 <- R11.m12 + (1/R11.lam[i1]) * (R11.e[,i1] %*% t(R11.e[,i1])) } ## solve(R11.m12 %*% R11.m12) ## check that R11^(-1/2) works R22.lam <- eigen(R22)$val R22.e <- eigen(R22)$vec R22.m12 <- matrix(rep(0,q^2),ncol=q) R22.m1 <- matrix(rep(0,q^2),ncol=q) for (i1 in 1:q) { R22.m12 <- R22.m12 + (1/sqrt(R22.lam[i1])) * (R22.e[,i1] %*% t(R22.e[,i1])) R22.m1 <- R22.m1 + (1/R22.lam[i1]) * (R22.e[,i1] %*% t(R22.e[,i1])) } CC.mat1 <- R11.m12 %*% R12 %*% R22.m1 %*% R21 %*% R11.m12 CC.mat2 <- R22.m12 %*% R21 %*% R11.m1 %*% R12 %*% R22.m12 CC.mat1.lam <- eigen(CC.mat1)$val; CC.mat1.e <- eigen(CC.mat1)$vec CC.mat2.lam <- eigen(CC.mat2)$val; CC.mat2.e <- eigen(CC.mat2)$vec (a1 <- t(R11.m12) %*% CC.mat1.e[,1]) (b1 <- t(R22.m12) %*% CC.mat2.e[,1]) t(a1) %*% R11 %*% a1 t(b1) %*% R22 %*% b1 CC.mat1.lam sqrt(CC.mat1.lam)