pdf("F:\\Rmisc\\rpd13_1r.pdf") x0 <-rep(1,20) x1 <- c(rep(seq(20,29),2)) x2 <- x1-25; x2[1] <- -4 x2[11] <- -4 x3 <- rep(c(5,4,3,2,1,2,3,4,5,6),2) x <- cbind(x0,x1,x2,x3) x <- as.matrix(x) xx <- t(x) %*% x xxi <- solve(xx) n <- nrow(x) xiv <- cbind(x1,x2,x3) xiv <- as.matrix(xiv) imat <- diag(n) jmat <- matrix(rep(1/n,n*n),ncol=n) z1 <- (imat-jmat) %*% xiv z2 <- diag(t(z1) %*% z1) z2 <- sqrt(z2) z2 <- diag(z2) z <- z1 %*% solve(z2) rho_hat <- t(z) %*% z rho_hat eigenrho <- eigen(rho_hat) eigval <- eigenrho$val eigvec <- eigenrho$vec print( cbind(eigval,eigvec)) w1 <- z %*% eigvec[,1] w2 <- z %*% eigvec[,2] w3 <- z %*% eigvec[,3] v <- eigvec l1_2 <- sqrt(diag(eigval)) l12_v <- l1_2 %*% t(v) l12_v l12_vpd2 <- t(l12_v[1:2,]) l12_vpd2 u <- z %*% v %*% solve(l1_2) u ud2 <- u[,1:2] theta12 <- 57.2958*acos(sum(l12_vpd2[1,] * t(l12_vpd2[2,]))) theta13 <- 57.2958*acos(sum(l12_vpd2[1,] * t(l12_vpd2[3,]))) theta23 <- 57.2958*acos(sum(l12_vpd2[2,] * t(l12_vpd2[3,]))) print(cbind(theta12,theta13,theta23)) l12_vpd2 <- cbind(l12_vpd2,matrix(rep(1,nrow(l12_vpd2)),ncol=1)) ud2 <- cbind(ud2,matrix(rep(2,nrow(ud2)),ncol=1)) lvu <- rbind(l12_vpd2,ud2) plot(ud2[,1],ud2[,2],xlim=c(-1,1),ylim=c(-1,1)) abline(h=0) abline(v=0) for (i in 1:nrow(l12_vpd2)) arrows(0,0,l12_vpd2[i,1],l12_vpd2[i,2]) dev.off()