pdf("rpd2_13.pdf")


X1 <- matrix(
c(11.1,25.2,33.0,1576,12.6,26.0,48.2,1767,
9.3,27.3,27.3,1580,9.8,26.7,34.3,1578,
5.8,25.7,27.3,1222,12.1,28.3,20.4,1926),byrow=T,ncol=4)

Xmn <- colMeans(X1)

Xmeans <- rbind(t(Xmn),t(Xmn),t(Xmn),t(Xmn),t(Xmn),t(Xmn))

Xss <- c()

for (i in 1:ncol(X1)) {
Xss[i] <- (nrow(X1)-1)*var(X1[,i])
}

Xsse <- rbind(t(Xss),t(Xss),t(Xss),t(Xss),t(Xss),t(Xss))
X <- (X1-Xmeans)/sqrt(Xsse)

svdX <- svd(X)

U <- svdX$u
Z <- svdX$v
lambdasr <- svdX$d

Xchk <- U %*% diag(lambdasr) %*% t(Z)

L <- diag(lambdasr *lambdasr)

A1 <- lambdasr[1]*U[,1]%*%t(Z[,1])
gfA1X <- L[1,1]/sum(L)
A2 <- lambdasr[2]*U[,2]%*%t(Z[,2])
A1A2 <- A1+A2
gfA1A2X <- (L[1,1]+L[2,2])/sum(L)
A3 <- lambdasr[3]*U[,3]%*%t(Z[,3])
A1A2A3 <- A1+A2+A3
gfA1A2A3X <- (L[1,1]+L[2,2]+L[3,3])/sum(L)
A4 <- lambdasr[4]*U[,4]%*%t(Z[,4])
A1A2A3A4 <- A1+A2+A3+A4
gfA1A2A3A4X <- (L[1,1]+L[2,2]+L[3,3]+L[4,4])/sum(L)

X
U 
lambdasr
Z
Xchk
L
A1
gfA1X
A2
A1A2
gfA1A2X
A3
A1A2A3
gfA1A2A3X
A4
A1A2A3A4
gfA1A2A3A4X


dev.off()