#案例9.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))
X=as.matrix(DA[,1:4])
Y=as.matrix(DA[,5])
G=(t(Y)%*%Y-t(Y)%*%X%*%solve(t(X)%*%X)%*%t(X)%*%Y)/8#计算sigma
A1=sum(resid(lm(Y~X1,data=DA))^2)
A2=sum(resid(lm(Y~X2,data=DA))^2)
A3=sum(resid(lm(Y~X3,data=DA))^2)
A4=sum(resid(lm(Y~X4,data=DA))^2)
A5=sum(resid(lm(Y~X1+X2,data=DA))^2)
A6=sum(resid(lm(Y~X1+X3,data=DA))^2)
A7=sum(resid(lm(Y~X1+X4,data=DA))^2)
A8=sum(resid(lm(Y~X2+X3,data=DA))^2)
A9=sum(resid(lm(Y~X2+X4,data=DA))^2)
A10=sum(resid(lm(Y~X3+X4,data=DA))^2)
A11=sum(resid(lm(Y~X1+X2+X3,data=DA))^2)
A12=sum(resid(lm(Y~X1+X2+X4,data=DA))^2)
A13=sum(resid(lm(Y~X1+X3+X4,data=DA))^2)
A14=sum(resid(lm(Y~X2+X3+X4,data=DA))^2)
A15=sum(resid(lm(Y~X1+X2+X3+X4,data=DA))^2)
A=c(A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15)#RSS残差平方和
N=c(2,2,2,2,3,3,3,3,3,3,4,4,4,4,5)#变量个数（含常数项）
solve(diag(13-N))%*%A#RMS计算
A/G-13+2*N#Cp计算

#Cp图
q=c(2,2,2,2,3,3,3,3,3,3,4,4,4,4,5)
Cp=c(183.47,128.82,285.91,125.41,1.80,179.59,4.37,56.17,125.13,19.72,2.32,2.29,2.73,6.23,4.28)
plot(q,Cp,xlim=c(0,5),ylim=c(0,7),type="n")
text(q,Cp)
lines(c(0,5),c(0,5))

#案例9.2
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)
lm.step=step(lm.h)
summary(lm.step)
drop1(lm.step)

lm.h=lm(Y~X1+X2,data=DA)
summary(lm.h)
lm.step=step(lm.h)

#案例9.3
rm(list=ls())
read.csv("szzs.csv",header=T)->DA
lm.sz=lm(cf~zx+cy+bz+ll+cj+zc+sp+fz+mc+cz+sh+dz+js+jx+yy+sd+jz+ys+it+pl+jr+dc+fw+cb+hs+jj,data=DA)
summary(lm.sz)
lm.step=step(lm.sz)
summary(lm.step)
drop1(lm.step)#单项删除
lm.o=lm(cf~ll+cj+zc+sp+mc+cz+sh+dz+js+jx+sd+jz+it+pl+jr+dc+fw+cb+hs+jj,data=DA)
summary(lm.o)
drop1(lm.o)
lm.op=lm(cf~ll+cj+zc+sp+mc+cz+sh+dz+js+jx+sd+it+pl+jr+dc+fw+cb+hs+jj,data=DA)
summary(lm.op)

#案例9.4
library(lars)
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))
x=DA[,1:4]
y=DA[,5]
X<-scale(x,scale=F)
Y<-scale(y,scale=F)
lm.l<-lars(X,Y,type="lar");lm.l
plot.lars(lm.l)#作图
summary(lm.l)
coef.lars(lm.l)

#案例9.5
rm(list=ls())
library(lars)
DA=read.csv("hs300.csv",header=T)
x=DA[,1:300]
y=DA[,301]
X<-scale(x,scale=F)
Y<-scale(y,scale=F)
lm.l<-lars(X,Y,type="lar");lm.l
plot.lars(lm.l)#作图
summary(lm.l)
min(lm.l$Cp)#给出Cp的最小值
E=coef.lars(lm.l)[176,]#提取系数
a<-0
for(j in 1:300)
{if(E[j]==0)a=a+1}
m=300-a#变量子集含有的变量个数
#残差图
d.bar=colMeans(DA)
x.bar=colMeans(DA[,1:300])
y.bar=d.bar[301]
beta0=y.bar-E%*%x.bar#计算常数项
x<-as.matrix(x)
r<- y-beta0-x%*%E
Res<-t(r)
Time<-1:480
plot(Time,Res,type="l")
lines(c(-15,515),c(0,0))
Res%*%r

#图9.4 选择s使得BIC达到最小
s=seq(0.01,15,0.01)
b=c(1:1500)
for(i in 1:1500)
{coef<-predict(lm.l,s=0.01*i,type="coef",mode="lambda")
z<-coef$coefficients
z=as.matrix(z)
a<-0
for(j in 1:300)
{if(z[j,1]==0)a=a+1}
m=300-a
A=as.matrix(Y-X%*%z)
B=as.matrix(t(A))
BIC=log(B%*%A)+log(300)*m/300
b[i]<-BIC}
plot(s,b,type="l")
min(b)#最小BIC值
coef<-predict(lm.l,s=0.47,type="coef",mode="lambda")
z<-coef$coefficients
z=as.matrix(z)
a<-0
for(j in 1:300)
{if(z[j,1]==0)a=a+1}
m=300-a;m
A=as.matrix(Y-X%*%z)
B=as.matrix(t(A))
BIC=log(B%*%A)+log(300)*m/300

#残差图
d.bar=colMeans(DA)
x.bar=colMeans(DA[,1:300])
y.bar=d.bar[301]
E=coef.lars(lm.l)[123,]#提取系数
beta0=y.bar-E%*%x.bar#计算常数项
x<-as.matrix(x)
r<- y-beta0-x%*%E
Res<-t(r)
Time<-1:480
plot(Time,Res,type="l")
lines(c(-15,515),c(0,0))
Res%*%r