#案例4.1
rm(list=ls())
X=matrix(c(2,0,1,2,1,0,0,1,2),3,3)
X=scale(X)
XX<-t(X)%*%X
eigen(XX)#求矩阵的特征值
X%*%eigen(XX)$vectors

#用R函数
C=princomp(X,cor=TRUE)
summary(C, loadings=TRUE)
#
princomp(covmat=XX)
loadings(princomp(covmat=XX))

Y<-matrix(c(1,2,2,5,2,0,3,1,3,1,0,1,0,3,1,2),4,4)
Y=scale(Y)
YY<-t(Y)%*%Y  #矩阵转置
eigen(YY)#求矩阵的特征值
Y%*%eigen(YY)$vectors

D0<-matrix(c(18.24,17.98,18.04,18.04,18.89,19.3,18.85,18.74,18.51,10.61,10.28,10.33,10.46,10.81,10.93,10.81,10.63,10.46,12.37,11.68,11.66,11.76,11.79,11.96,11.8,11.56,11.42,15.52,15.3,15.36,15.57,15.61,15.76,15.84,15.67,15.13,19.67,19.42,19.21,17.96,18.09,18.34,18.26,18.34,18.03),9,5)
D=scale(D0)
DD<-t(D)%*%D  
eigen(DD)#求矩阵的特征值
D%*%eigen(DD)$vectors

#另外直接用R函数求样本主成分的预测值
pr=princomp(D);pr
pre<-predict(pr);pre

pr=princomp(covmat=DD);pr#使用协方差矩阵
pre<-predict(pr);pre

#案例4.2:大样本数据矩阵
rm(list=ls())
read.csv("20140930.csv",header=T)->Z0
Z1<-Z0[,3:22]
Z=scale(Z1)
ZZ<-t(Z)%*%Z
eigen(ZZ)#求矩阵的特征值
ZZ%*%eigen(ZZ)$vectors

pr=princomp(Z);pr
pre<-predict(pr);pre
screeplot(pr,npcs=10,type=c("lines"))#主成分碎石图
biplot(pr)#缺省的是第一、二主成分散点图
biplot(pr,choices=3:4,scale=1)

#结果输出到文件
write.table(pre, file = "a1.txt", row.names = F, quote = F)
write.csv(pre, file = "a1.csv", row.names = F, quote = F)


pr=princomp(covmat=ZZ);pr#使用协方差验证标准差与直接特征值计算的一致性