##数据预处理
df<-read.csv(file="旅游板块和航空板块数据(20210210).csv")         
dim(df)
head(df)
X<-df[,-1]
y<-c(rep(1,25),rep(2,29))
n<-length(y)
rownames(X)<-df[,1]
group1<-1:20
group2<-26:50
X.train1<-X[group1,]
X.train2<-X[group2,]
y.train<-y[c(group1,group2)]
n1<-length(y.train)
X.test<-X[-c(group1,group2),]
y.test<-y[-c(group1,group2)]


##距离判别
#两总体判别分析
discriminiant.distance <- function(TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){
  if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
  if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
  if (is.matrix(TstX) != TRUE) TstX <- as.matrix(TstX)
  if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
  if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
  nx <- nrow(TstX)
  blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE,
                  dimnames=list("blong", 1:nx))
  mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
  if (var.equal == TRUE || var.equal == T){
    S <- var(rbind(TrnX1,TrnX2))
    w <- mahalanobis(TstX, mu2, S)- mahalanobis(TstX, mu1, S)
  }
  else{
    S1 <-var(TrnX1); S2 <- var(TrnX2)
    w <- mahalanobis(TstX, mu2, S2)- mahalanobis(TstX, mu1, S1)
  }
  for (i in 1:nx){
    if (w[i] > 0)
      blong[i] <- 1
    else
      blong[i] <- 2
  }
  blong
}

distance.nei11<-discriminiant.distance(X.train1,X.train2,var.equal=TRUE)
sum(diag(table(distance.nei11,y.train)))/length(y.train)
distance.wai11<-discriminiant.distance(X.train1,X.train2,TstX = X.test,var.equal=TRUE)
sum(diag(table(distance.wai11,y.test)))/length(y.test)
distance.wai11

distance.nei12<-discriminiant.distance(X.train1,X.train2)
sum(diag(table(distance.nei12,y.train)))/length(y.train)
distance.wai12<-discriminiant.distance(X.train1,X.train2,TstX = X.test)
sum(diag(table(distance.wai12,y.test)))/length(y.test)
distance.wai12

##Bayes判别
#两总体判别
discriminiant.bayes <- function(TrnX1, TrnX2, rate = 1, TstX = NULL, var.equal = FALSE){
  if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
  if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
  if (is.matrix(TstX) != TRUE)  TstX <- as.matrix(TstX)
  if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
  if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
  nx <- nrow(TstX)
  blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE,
                  dimnames=list("blong", 1:nx))
  mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
  if (var.equal == TRUE || var.equal == T){
    S <- var(rbind(TrnX1,TrnX2)); beta <- 2*log(rate)
    w <- mahalanobis(TstX, mu2, S)- mahalanobis(TstX, mu1, S)
  }
  else{
    S1 <- var(TrnX1); S2 <- var(TrnX2)
    beta <- 2*log(rate) + log(det(S1)/det(S2))
    w <- mahalanobis(TstX, mu2, S2)- mahalanobis(TstX, mu1, S2)
  }
  for (i in 1:nx){
    if (w[i] > beta)
      blong[i] <- 1
    else
      blong[i] <- 2
  }
  blong
}

distance.nei21<-discriminiant.bayes(X.train1,X.train2,rate=1,var.equal=TRUE)
sum(diag(table(distance.nei21,y.train)))/length(y.train)
distance.wai21<-discriminiant.bayes(X.train1,X.train2,TstX = X.test,rate=1,var.equal=TRUE)
sum(diag(table(distance.wai21,y.test)))/length(y.test)
distance.wai21

distance.nei22<-discriminiant.bayes(X.train1,X.train2,rate=1)
sum(diag(table(distance.nei22,y.train)))/length(y.train)
distance.wai22<-discriminiant.bayes(X.train1,X.train2,TstX = X.test,rate=1)
sum(diag(table(distance.wai22,y.test)))/length(y.test)
distance.wai22

##Fisher判别
#两总体判别分析
discriminiant.fisher <- function(TrnX1, TrnX2, TstX = NULL){
  if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
  if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
  if (is.matrix(TstX) != TRUE)  TstX <- as.matrix(TstX)
  if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
  if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
  nx <- nrow(TstX)
  blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE,
                  dimnames=list("blong", 1:nx))
  n1 <- nrow(TrnX1); n2 <- nrow(TrnX2)
  mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
  S <- (n1-1)*var(TrnX1) + (n2-1)*var(TrnX2)
  mu <- n1/(n1+n2)*mu1 + n2/(n1+n2)*mu2
  w <- (TstX-rep(1,nx) %o% mu) %*% solve(S, mu2-mu1);
  for (i in 1:nx){
    if (w[i] <= 0)
      blong[i] <- 1
    else
      blong[i] <- 2
  }
  blong
}

distance.nei31<-discriminiant.fisher(X.train1,X.train2)
sum(diag(table(distance.nei31,y.train)))/length(y.train)
distance.wai31<-discriminiant.fisher(X.train1,X.train2,TstX = X.test)
sum(diag(table(distance.wai31,y.test)))/length(y.test)
distance.wai31

#稳健的稀疏判别
#install.packages("glmnet")
#install.packages("quantreg")
library(glmnet)
library(quantreg)
library(MASS)
data_C1<-X.train1
data_C2<-X.train2
p<-ncol(data_C2)
x1<-data_C1
x2<-data_C2
mu1<-colMeans(x1)
mu2<-colMeans(x2)
X<-as.matrix(rbind(x1,x2))
nzc<-nrow(X)
y<-c(rep(-2,20),rep(2,25))
cvob1=cv.glmnet(X,y,alpha=1)
beta_lasso<-as.numeric(coef(cvob1,s=cvob1$lambda.min))
beta_lasso<-beta_lasso[-1]
beta_lasso[which(abs(beta_lasso)<10^-8)]=0
#一乘
lambda_q<-log(2*nzc)/sum(abs(beta_lasso))
cvob2=rq.fit.lasso(X,y,tau=0.5, lambda=lambda_q*1000)
beta_qlasso<-as.numeric(coef(cvob2))
beta_qlasso[which(abs(beta_qlasso)<10^-8)]=0

if(t(mu2-mu1)%*%beta_lasso<=0){
beta_lasso<--beta_lasso
print("调整系数")
}
if(t(mu2-mu1)%*%beta_qlasso<=0){
beta_qlasso<--beta_qlasso
print("调整系数")
}

beta0_hat<-sum(-(mu1+mu2)*beta_lasso/2)
beta0_qhat<-sum(-(mu1+mu2)*beta_qlasso/2)

##训练集判别结果
#旅行股判别
n_test=20
result<-0
resultq<-0
for(z in 1:n_test)
{
x_test<-X.train1[z,]
D<-sum(x_test*beta_lasso)+beta0_hat
Dq<-sum(x_test*beta_qlasso)+beta0_qhat
if(D>0){result<-result+1}
if(Dq>0){resultq<-resultq+1}
}
result/n_test
resultq/n_test

#航空股判别
n_test=25
result<-0
resultq<-0
for(z in 1:n_test)
{
x_test<-X.train2[z,]
D<-sum(x_test*beta_lasso)+beta0_hat
Dq<-sum(x_test*beta_qlasso)+beta0_qhat
if(D>0){result<-result+1}
if(Dq>0){resultq<-resultq+1}
}
result/n_test
resultq/n_test

##测试集判别结果
#旅行股判别
n_test=5
result<-0
resultq<-0
for(z in 1:n_test)
{
x_test<-X.test[z,]
D<-sum(x_test*beta_lasso)+beta0_hat
Dq<-sum(x_test*beta_qlasso)+beta0_qhat
if(D>0){result<-result+1}
if(Dq>0){resultq<-resultq+1}
}
result/n_test
resultq/n_test

#航空股判别
n_test=4
result<-0
resultq<-0
for(z in 1:n_test)
{
x_test<-X.test[5+z,]
D<-sum(x_test*beta_lasso)+beta0_hat
Dq<-sum(x_test*beta_qlasso)+beta0_qhat
if(D>0){result<-result+1}
if(Dq>0){resultq<-resultq+1}
}
result/n_test
resultq/n_test

#对应调整参数变化的弹性网拟合系数曲线（图6.1）
library(glmnet)
library(quantreg)
##随机数生成
x=matrix(rnorm(100*20),100,20)
y=rnorm(100)
fit1=glmnet(x,y)
predict(fit1,newx=x[1:5,],s=c(0.01,0.005))
predict(fit1,type="coef")
plot(fit1,xvar="lambda",xlab="log(λ)",ylab="系数")