#creating the hypothetical data set of pairs of data from a multivariate
#normal distribution with N(100, 100), N(50, 50) and rho=0.707.
rm2 <- rmvnorm(100, mean=c(100, 50), cov=matrix(c(100, 50, 50, 50), 2)) 

#verifying that we used the correct cov matrix.
print( matrix(c(100, 50, 50, 50), 2) )

#Fig. 14.1, p. 332.
plot(rm2[, 1], rm2[, 2], xlab="Norm1", ylab="Norm2")

#Table 14.1, p. 333.
apply( rm2, 2, mean)
apply( rm2, 2, stdev)
apply( rm2, 2, var)
cor(rm2)

#Principal components, p. 335.
pr1 <- princomp(rm2)
loadings(pr1)
#eigenvalues:
pr1.eigen <- pr1$sdev^2
print(pr1.eigen)

#use a subset of the depress data set
attach(depress)
cesd.subset <- data.frame(c1, c2, c3, c4, c5, c6, c7, c8, c9, c10,
c11, c12, c13, c14, c15, c16, c17, c18, c19, c20)
detach(depress)

#need to standardize the variables
std.cesd <- cesd.subset
std.cesd$c1 <- cesd.subset$c1/stdev(cesd.subset$c1)
std.cesd$c2 <- cesd.subset$c2/stdev(cesd.subset$c2)
std.cesd$c3 <- cesd.subset$c3/stdev(cesd.subset$c3)
std.cesd$c4 <- cesd.subset$c4/stdev(cesd.subset$c4)
std.cesd$c5 <- cesd.subset$c5/stdev(cesd.subset$c5)
std.cesd$c6 <- cesd.subset$c6/stdev(cesd.subset$c6)
std.cesd$c7 <- cesd.subset$c7/stdev(cesd.subset$c7)
std.cesd$c8 <- cesd.subset$c8/stdev(cesd.subset$c8)
std.cesd$c9 <- cesd.subset$c9/stdev(cesd.subset$c9)
std.cesd$c10 <- cesd.subset$c10/stdev(cesd.subset$c10)

std.cesd$c11 <- cesd.subset$c11/stdev(cesd.subset$c11)
std.cesd$c12 <- cesd.subset$c12/stdev(cesd.subset$c12)
std.cesd$c13 <- cesd.subset$c13/stdev(cesd.subset$c13)
std.cesd$c14 <- cesd.subset$c14/stdev(cesd.subset$c14)
std.cesd$c15 <- cesd.subset$c15/stdev(cesd.subset$c15)
std.cesd$c16 <- cesd.subset$c16/stdev(cesd.subset$c16)
std.cesd$c17 <- cesd.subset$c17/stdev(cesd.subset$c17)
std.cesd$c18 <- cesd.subset$c18/stdev(cesd.subset$c18)
std.cesd$c19 <- cesd.subset$c19/stdev(cesd.subset$c19)
std.cesd$c20 <- cesd.subset$c20/stdev(cesd.subset$c20)

#table 14.2, p. 342.
pr.std <- princomp(std.cesd)
pr.std$coef
pr.std.eigen <- pr.std$sdev^2
print(pr.std.eigen)
r.pr <- 0.5/sqrt(pr.std.eigen)
print(r.pr)

#Fig. 14.5, p. 343.
plot(pr.std.eigen)