#案例8.1 水泥数据
DA=data.frame(
X1=c( 7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10),
X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68),
X3=c( 6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8),
X4=c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12),
Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5,
93.1,115.9, 83.8, 113.3, 109.4))
lm.h=lm(Y~X1+X2+X3+X4,data=DA)
summary(lm.h)
#作主成分分析
D=scale(DA[,1:4])
R=cov(D)
eigen(R)
pr=princomp(D)
summary(pr,loadings=TRUE)
# 预测样本主成分
D[,1:4]%*%eigen(R)$vectors[,1:3]
pre<-predict(pr);pre#或用主成分预测值计算
#建立主成分回归方程
DA$z1<-pre[,1]
DA$z2<-pre[,2]
DA$z3<-pre[,3]
lm.z<-lm(DA$Y~DA$z1+DA$z2+DA$z3)
summary(lm.z)
# 作变换, 得到原坐标下的关系表达式
beta<-coef(lm.z); A<-loadings(pr)
x.bar<-colMeans(DA[,1:4]); x.sd<-sapply(DA[,1:4], sd)
coef<-(beta[2]*A[,1]+ beta[3]*A[,2]+ beta[4]*A[,3])/x.sd
beta0 <- beta[1]- sum(x.bar * coef)
c(beta0, coef)

#案例8.2 选择三个主成分的计算结果
rm(list=ls())
DA=read.csv("hs300.csv",header=T)
#### 主成分分析
D=scale(DA[1:300])
R=cov(D)
eigen(R)$values
pr=princomp(D)
summary(pr, loadings=F)
# 预测样本主成分, 并作主成分分析
pre<-predict(pr)
DA$z1<-pre[,1]
DA$z2<-pre[,2]
DA$z3<-pre[,3]
lm.z<-lm(DA$Y~DA$z1+DA$z2+DA$z3)
summary(lm.z)
#### 作变换, 得到原坐标下的关系表达式
beta<-coef(lm.z); A<-loadings(pr)
x.bar<-colMeans(DA[,1:300]); x.sd<-sapply(DA[,1:300], sd)
coef<-(beta[2]*A[,1]+ beta[3]*A[,2]+ beta[4]*A[,3])/x.sd
beta0 <- beta[1]- sum(x.bar * coef)
c(beta0, coef)
#计算残差
x<-as.matrix(DA[,1:300])
coef=as.vector(coef)
y<-DA[,301]
r<- y-beta0-x%*%coef
Res<-t(r)
Time<-1:480
plot(Time,Res,type="l")
lines(c(-15,515),c(0,0))
Res%*%r

#正确的主成分估计，舍弃10-4特征值对应的主成分，剩余252个主成分。
rm(list=ls())
DA=read.csv("hs300.csv",header=T)
# 主成分分析
D=scale(DA[1:300])
R=cov(D)
eigen(R)$values
pr=princomp(D)
summary(pr, loadings=F)
# 预测样本主成分, 并作主成分分析
pre<-predict(pr)
write.csv(pre,"pre.csv")
DAA=read.csv("pre.csv",header=T)
DAA=DAA[2:301]
lm.z<-lm(DA$Y~.,data=DAA)#不要常数项（去掉0+就是要常数项）
summary(lm.z)
#### 作变换, 得到原坐标下的关系表达式
beta<-coef(lm.z); A<-loadings(pr)
x.bar<-colMeans(DA[,1:300]); x.sd<-sapply(DA[,1:300], sd)
beta1=beta[2:301]
for (i in 253:300) beta1[i]=0#舍弃后面48个主成分
coef=A%*%beta1/x.sd
beta0 <- beta[1]- sum(x.bar * coef)
c(beta0, coef)
#计算残差
x<-as.matrix(DA[,1:300])
coef=as.vector(coef)
y<-DA[,301]
r<- y-beta0-x%*%coef
Res<-t(r)
Time<-1:480
plot(Time,Res,type="l")
lines(c(-15,515),c(0,0))
Res%*%r

#例8.3 岭估计 需安装MASS程序包
rm(list=ls())
library(MASS)
DA=read.csv("szzs.csv",header=T)
#岭迹图
plot(lm.ridge(cf~.,data=DA[,2:28],lambda=seq(0,0.04,0.0001)))
#岭参数选取
select(lm.ridge(cf~.,data=DA[,2:28],lambda=seq(0,0.04,0.0001)))
#岭系数
lm.r=lm.ridge(cf~.,data=DA[,2:28],lambda=0.0006)
lm.r

#例8.4
rm(list=ls())
library(MASS)
DA=read.csv("hs300.csv",header=T)
#岭迹图
plot(lm.ridge(Y~.,data=DA,lambda=seq(0,0.3,0.001)))
#岭参数选取
select(lm.ridge(Y~.,data=DA,lambda=seq(0,3,0.001)))
#岭系数
lm.r=lm.ridge(Y~.,data=DA,lambda=0.05)
lm.r

#计算残差
x<-as.matrix(DA[,1:300])
coef<- coef(lm.r)
coef=as.vector(coef[2:301])
y<-DA[,301]
r<- y-712.78-x%*%coef
Res<-t(r)
Time<-1:480
plot(Time,Res,type="l")
lines(c(-20,515),c(0,0))
Res%*%r

#因变量和预测值的拟合图
yv<-t(y)
pv<-t(712.78+x%*%coef)
Time<-1:480
plot(Time,yv,type="l",xlab="time",ylab="")
lines(Time,pv,xlab="time",ylab="")