library(MASS)#mvrnorm
library(car)#scatterplotMatrix；；；；vip(model,....)计算线性模型、广义线性模型和其他模型的方差-通货膨胀和广义方差-通货膨胀因子。
library(glmnet)#用惩罚极大似然拟合广义线性模型。在正则化参数lambda的数值网格上计算套索或elasticnet惩罚的正则化路径。可以处理各种形状的数据，包括非常大的稀疏数据矩阵。拟合线性，logistic和多项式，泊松和Cox回归模型。

library(psych)#corr.tes
library(ncvreg)#SCAD 拟合MCP-或scad -惩罚回归模型的系数路径在一个网格值的正则化参数lambda。拟合线性和逻辑回归模型，选择额外的L2惩罚。
library(msgps)#Adaptive Lasso基于广义寻径算法计算了套索、弹性网、广义弹性网和自适应套索的自由度。采用Mallows’Cp、偏置校正AIC (AICc)、广义交叉验证(GCV)和BIC等模型选择准则，可以筛选出最优模型。

#=================================数值模拟环境=============
#样本  n=80 #（训练集+测试集），预测变量维数 p=50，稀疏矩阵：beta=c(1,1.2,-3.2,4.1,-5.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
 #m=300 试验次数
#===================================X:多元正态数据集==========================-
options(digits=3)
p=50
n=80
mean=rnorm(p,0,1)
cname=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10"
,"x11","x12","x13","x14","x15","x16","x17","x18","x19","x20"
,"x21","x22","x23","x24","x25","x26","x27","x28","x29","x30"
,"x31","x32","x33","x34","x35","x36","x37","x38","x39","x40"
,"x41","x42","x43","x44","x45","x46","x47","x48","x49","x50")
#sigma08:生成相关系数为0.8的协方差矩阵（正定）50*50
sigma08=matrix(0,nrow=50,ncol=50,dimnames=list(cname,cname))
for (i in 1:50){
    for (j in 1:50){
    if (i==j) sigma08[i,j]=1
    else sigma08[i,j]=0.8}}
XD08=mvrnorm(n, mean, sigma08)

#------------------相关系数及相关性检验--------------------

corr.test(XD08,use="complete")#use=c(pairwise,complete)
#-------------------------------------------------------------------

beta=c(1,1.2,-3.2,4.1,-5.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
 
beta1=beta/(1+1.2^2+(-3.2)^2+4.1^2+(-5.)^2)^0.5
#0.135  0.163 -0.434  0.555 -0.677  0.000  0.000  0.000  0.000  0.000

train=sample(c(TRUE,FALSE),nrow(XD08),rep=TRUE)
test=(!train)

#-------------训练集--------------------
XD08Tr=XD08[train,]     # X矩阵
esp=rnorm(nrow(XD08Tr),0,1)
YD08Tr=XD08Tr%*%beta+esp   #Y矩阵
dim(XD08Tr)                
#-------------测试集---------------------------

XD08Te=XD08[test,]  #X矩阵
esp=rnorm(nrow(XD08Te),0,1)
YD08Te=XD08Te%*%beta+esp #Y矩阵
dim(XD08Te)
#------------散点图---------------------------------
YXD08Tr=cbind(YD08Tr,XD08Tr) #[Y,X]矩阵
dfYXD08Tr=as.data.frame(YXD08Tr)#数据框
names(dfYXD08Tr)=c("ymo","x1","x2","x3","x4","x5","x6","x7","x8","x9","x10"
,"x11","x12","x13","x14","x15","x16","x17","x18","x19","x20"
,"x21","x22","x23","x24","x25","x26","x27","x28","x29","x30"
,"x31","x32","x33","x34","x35","x36","x37","x38","x39","x40"
,"x41","x42","x43","x44","x45","x46","x47","x48","x49","x50")
attach(dfYXD08Tr)#[Y,X]数据框
scatterplotMatrix(~ymo+x1+x2+x50, data=dfYXD08Tr, spread=FALSE,
    lty.smooth=2, main="ScatterplotmydatadfYXD08Tr")
#----------------------------------------------------------------------------
#---------------elastic-net(alpha=0.5)--------------------
YD08Trs=as.vector(YD08Tr)
YD08Tes=as.vector(YD08Te)

fEND08=msgps(XD08Tr,YD08Trs,penalty="enet",alpha=0.5)#(penalty="enet",alpha=0---Lasso);(penalty="enet"/gent",0<alpha<1);(penalty="alasso",lambda(控制w）,gamma=“The value of γ on "alasso".)
plot.msgps(fEND08)
coef(fEND08)
summary.msgps(fEND08)


#err箱线图,角度~试验次数,运行时间
system.time({
errsumEND08=0
agsumEND08=0
m=300
errEND08=numeric(m)
agEND08=numeric(m)
for (i in 1:m){
mean=rnorm(p,0,1)
XD08=mvrnorm(n, mean, sigma08)#生成X多元正态矩阵数据集，以下为X的行命名

train=sample(c(TRUE,FALSE),nrow(XD08),rep=TRUE)
test=(!train)

#-------------训练集-------------------
XD08Tr=XD08[train,]     # X矩阵
esp=rnorm(nrow(XD08Tr),0,1)
YD08Tr=XD08Tr%*%beta+esp   #Y矩阵
dim(XD08Tr)                
#-------------测试集---------------------------

XD08Te=XD08[test,]  #X矩阵
esp=rnorm(nrow(XD08Te),0,1)
YD08Te=XD08Te%*%beta+esp #Y矩阵
dim(XD08Te)
YD08Trs=as.vector(YD08Tr)
YD08Tes=as.vector(YD08Te)

# 
cvEND08=msgps(XD08Tr,YD08Trs,penalty="enet",alpha=0.5)
predmi=predict(cvEND08,XD08Te,s="lambda.min")
errEND08[i]=mean((predmi-YD08Tes)^2)
errsumEND08=errsumEND08+sqrt(errEND08[i])
btLAi=as.vector(coef(cvEND08,s = "lambda.min"))
btj=btLAi[104:153]
btk=btj/(sum(btj^2))^0.5
agEND08[i]=(180/pi)*acos(beta1%*%btk)
agsumEND08=agsumEND08+agEND08[i]
     i=i+1
}
})

x=1:m

avagEND08=agsumEND08/m
averrEND08=errsumEND08/m
plot(x,agEND08,pch=20,lty=2,type="b",col="blue",main="ElaNET",xlab="试验序号",ylab="Angle")
abline(h=avagEND08,lty=1,col="red",lwd=2)

yEND08=sqrt(errEND08)
boxplot(yEND08,main="ElaNET-ERR-D08",ylab="err")

boxplot(agEND08,main="ElaNET-Angle-D08",ylab="Angle")

boxplot.stats(yEND08)
boxplot.stats(agEND08)

#-------------LASSO---------alpha=1--------------------------
YD08Trs=as.vector(YD08Tr)
YD08Tes=as.vector(YD08Te)


fLAD08=msgps(XD08Tr,YD08Trs)#(penalty="enet",alpha=0---Lasso);(penalty="enet"/gent",0<alpha<1);(penalty="alasso",lambda(控制w）,gamma=“The value of γ on "alasso".)
plot.msgps(fLAD08)
coef(fLAD08)
summary.msgps(fLAD08)
system.time({
errsumLAD08=0
agsumLAD08=0
m=300
errLAD08=numeric(m)
agLAD08=numeric(m)
for (i in 1:m){
mean=rnorm(p,0,1)
XD08=mvrnorm(n, mean, sigma08)#生成X多元正态矩阵数据集，以下为X的行命名

train=sample(c(TRUE,FALSE),nrow(XD08),rep=TRUE)
test=(!train)

#-------------训练集-------------------
XD08Tr=XD08[train,]     # X矩阵
esp=rnorm(nrow(XD08Tr),0,1)
YD08Tr=XD08Tr%*%beta+esp   #Y矩阵
dim(XD08Tr)                
#-------------测试集---------------------------

XD08Te=XD08[test,]  #X矩阵
esp=rnorm(nrow(XD08Te),0,1)
YD08Te=XD08Te%*%beta+esp #Y矩阵
dim(XD08Te)
YD08Trs=as.vector(YD08Tr)
YD08Tes=as.vector(YD08Te)


#-----------------------  
cvLAD08=msgps(XD08Tr,YD08Trs)
predmi=predict(cvLAD08,XD08Te,s="lambda.min")
errLAD08[i]=mean((predmi-YD08Tes)^2)
errsumLAD08=errsumLAD08+sqrt(errLAD08[i])
btLAi=as.vector(coef(cvLAD08,s = "lambda.min"))
btj=btLAi[104:153]
btk=btj/(sum(btj^2))^0.5
agLAD08[i]=(180/pi)*acos(beta1%*%btk)
agsumLAD08=agsumLAD08+agLAD08[i]
     i=i+1
}
})

x=1:m

avagLAD08=agsumLAD08/m
averrLAD08=errsumLAD08/m
plot(x,agLAD08,pch=20,lty=2,type="b",col="blue",main="ElaNET",xlab="试验序号",ylab="Angle")
abline(h=avagLAD08,lty=1,col="red",lwd=2)

yLAD08=sqrt(errLAD08)
boxplot(yLAD08,main="Lasso-ERR-D08",ylab="err")

boxplot(agLAD08,main="Lasso-Angle-D08",ylab="Angle")

boxplot.stats(yLAD08)
boxplot.stats(agLAD08)


#---------------------------#SCAD拟合--------------------------------------

fSDD08=ncvreg(XD08Tr,YD08Tr,family="gaussian",penalty="SCAD",gamma=3.7,lambda.min=0.05,alpha=1)
plot(fSDD08)

cvSDD08=cv.ncvreg(XD08Tr,YD08Tr,family="gaussian",penalty="SCAD",gamma=3.7,lambda.min=0.05,alpha=1)
plot(cvSDD08)
plot(cvSDD08$fit)
plot(cvSDD08$lambda,cvSDD08$cve,lty=1,pch=20,col="blue")

opar=par(no.readonly=TRUE)
par(mfrow=c(2,2))
plot(cvSDD08, type='cve')
plot(cvSDD08, type='scale')  # sigma hat
plot(cvSDD08, type='snr')
plot(cvSDD08, type='rsq')
par(opar)
fit=cvSDD08$fit
cvSDD08$lambda.min#具有最小交叉验证错误的lambda值coef(cvSDH08)
betasd <- fit$beta[,cvSDD08$min]#cvSDH08$min=0.151(对应于lambda。min的索引)
summary(cvSDD08)
#summary(cvSDH08)$min/lambda/r.squared#/对应于lambda.min的索引)/lambda序列/模型解释的方差比例，通过交叉验证估计

system.time({
errsumSDD08=0
agsumSDD08=0
m=300
errSDD08=numeric(m)
agSDD08=numeric(m)
for (i in 1:m){
mean=rnorm(p,0,1)
XD08=mvrnorm(n, mean, sigma08)#生成X多元正态矩阵数据集，以下为X的行命名

train=sample(c(TRUE,FALSE),nrow(XD08),rep=TRUE)
test=(!train)

#-------------训练集-------------------
XD08Tr=XD08[train,]     # X矩阵
esp=rnorm(nrow(XD08Tr),0,1)
YD08Tr=XD08Tr%*%beta+esp   #Y矩阵
dim(XD08Tr)                
#-------------测试集---------------------------

XD08Te=XD08[test,]  #X矩阵
esp=rnorm(nrow(XD08Te),0,1)
YD08Te=XD08Te%*%beta+esp #Y矩阵
dim(XD08Te)


cvSDD08=cv.ncvreg(XD08Tr,YD08Tr,family="gaussian",penalty="SCAD",gamma=3.7,lambda.min=0.05,alpha=1)
predmi=predict(cvSDD08,XD08Te,s="lambda.min")
errSDD08[i]=mean((predmi-YD08Te)^2)
errsumSDD08=errsumSDD08+sqrt(errSDD08[i])
btLAi=as.vector(coef(cvSDD08,s = "lambda.min"))
btj=btLAi[2:51]
btk=btj/(sum(btj^2))^0.5
agSDD08[i]=(180/pi)*acos(beta1%*%btk)
agsumSDD08=agsumSDD08+agSDD08[i]
     i=i+1
}
})

x=1:m

avagSDD08=agsumSDD08/m
averrSDD08=errsumSDD08/m
plot(x,agSDD08,pch=20,lty=2,type="b",col="blue",main="ElaNET",xlab="试验序号",ylab="Angle")
abline(h=avagSDD08,lty=1,col="red",lwd=2)

ySDD08=sqrt(errSDD08)
boxplot(ySDD08,main="SCAD-ERR-D08",ylab="err")

boxplot(agSDD08,main="SCAD-Angle-D08",ylab="Angle")

boxplot.stats(ySDD08)
boxplot.stats(agSDD08)

#===================Adaptive  Lasso---===============================
YD08Trs=as.vector(YD08Tr)
YD08Tes=as.vector(YD08Te)

ALD08=msgps(XD08Tr,YD08Trs,penalty="alasso",gamma=1,lambda=0.001)#(penalty="enet",alpha=0---Lasso);(penalty="enet"/gent",0<alpha<1);(penalty="alasso",lambda(控制w）,gamma=“The value of γ on "alasso".)
plot.msgps(ALD08)
coef(ALD08)
summary.msgps(ALD08)


system.time({
errsumALD08=0
agsumALD08=0
m=300
errALD08=numeric(m)
agALD08=numeric(m)
for (i in 1:m){
mean=rnorm(p,0,1)
XD08=mvrnorm(n, mean, sigma08)#生成X多元正态矩阵数据集，以下为X的行命名

train=sample(c(TRUE,FALSE),nrow(XD08),rep=TRUE)
test=(!train)

#-------------训练集-------------------
XD08Tr=XD08[train,]     # X矩阵
esp=rnorm(nrow(XD08Tr),0,1)
YD08Tr=XD08Tr%*%beta+esp   #Y矩阵
YD08Trs=as.vector(YD08Tr)
dim(XD08Tr)                
#-------------测试集---------------------------

XD08Te=XD08[test,]  #X矩阵
esp=rnorm(nrow(XD08Te),0,1)
YD08Te=XD08Te%*%beta+esp #Y矩阵
YD08Tes=as.vector(YD08Te)
dim(XD08Te)

ALD08=msgps(XD08Tr,YD08Trs,penalty="alasso",gamma=1,lambda=0.001)#(penalty="enet",alpha=0---Lasso);(penalty="enet"/gent",0<alpha<1);(penalty="alasso",lambda(控制w）,gamma=“The value of γ on "alasso".)
predmi=predict(ALD08,XD08Te,s="lambda.min")
errALD08[i]=mean((predmi-YD08Tes)^2)
errsumALD08=errsumALD08+sqrt(errALD08[i])
btLAi=as.vector(coef(ALD08,s = "lambda.min"))
btj=btLAi[104:153]
btk=btj/(sum(btj^2))^0.5
agALD08[i]=(180/pi)*acos(beta1%*%btk)
agsumALD08=agsumALD08+agALD08[i]
     i=i+1
}
})

x=1:m

avagALD08=agsumALD08/m
averrALD08=errsumALD08/m
plot(x,agALD08,pch=20,lty=2,type="b",col="blue",main="Alasso",xlab="试验序号",ylab="Angle")
abline(h=avagALD08,lty=1,col="red",lwd=2)

yALD08=sqrt(errALD08)
boxplot(yALD08,main="Alasso-ERR-D08",ylab="err")

boxplot(agALD08,main="Alasso-Angle-D08",ylab="Angle")

boxplot.stats(yALD08)
boxplot.stats(agALD08)


#========================ERR箱线图==============================================
par(opar)

par(mfrow=c(2,2))
boxplot(yEND08,main="ElaNet-ERR-D08",ylab="err")
boxplot(yLAD08,main="Lasso-ERR-D08",ylab="err")
boxplot(ySDD08,main="SCAD-ERR-D08",ylab="err")
boxplot(yALD08,main="Alasso-ERR-D08",ylab="err")
par(opar)

boxplot.stats(yEND08)

boxplot.stats(yLAD08)

boxplot.stats(ySDD08)

boxplot.stats(yALD08)



#========================Angle箱线图==============================================

par(mfrow=c(2,2))
boxplot(agEND08,main="ElaNet-Angle-D08",ylab="Angle")
boxplot(agLAD08,main="Lasso-Angle-D08",ylab="Angle")
boxplot(agSDD08,main="SCAD-Angle-D08",ylab="Angle")
boxplot(agALD08,main="Alasso-Angle-D08",ylab="Angle")
par(opar)

boxplot.stats(agEND08)

boxplot.stats(agLAD08)

boxplot.stats(agSDD08)

boxplot.stats(agALD08)