##数据预处理
df<-read.csv(file="银行板块和航空板块数据(20210210).csv")         
dim(df)
head(df)
X<-df[,-1]
X1<-X
rownames(X1)<-df[,1]
##层次聚类
DisMatrix1<-dist(X1,method="manhattan")     #绝对距离
DisMatrix2<-dist(X1,method="euclidean")     #欧式距离
DisMatrix3<-dist(X1,method="maximum")       #切比雪夫距离
DisMatrix4<-dist(X1,method="canberra")      #兰氏距离
pingjia<-function(clustering,realclass){
  AA<-table(clustering,realclass)
  max(sum(diag(AA)),sum(AA)-sum(diag(AA)))/sum(AA)
}
clust<-function(DisMatrix,realclass){
  meth<-c("single","complete","centroid","average","ward.D")
  result<-0
  for (i in 1:5) {
    fit.hc<-hclust(d=DisMatrix,method=meth[i])    #类最短
    result[i]<-pingjia(cutree(fit.hc,k=2),realclass)
  }
  return(result)
}

clust(DisMatrix1,c(rep(1,25),rep(2,29)))
clust(DisMatrix2,c(rep(1,25),rep(2,29)))
clust(DisMatrix3,c(rep(1,25),rep(2,29)))
clust(DisMatrix4,c(rep(1,25),rep(2,29)))

fit.hc1<-hclust(d=DisMatrix4,method="complete")    
cutree(fit.hc1,k=2)
pingjia(cutree(fit.hc1,k=2),c(rep(1,25),rep(2,29)))
plot(fit.hc1,ylab="高度")
rect.hclust(fit.hc1,k=2,border = 2)

fit.hc2<-hclust(d=DisMatrix4,method="average")    
cutree(fit.hc2,k=2)
pingjia(cutree(fit.hc2,k=2),c(rep(1,25),rep(2,29)))
plot(fit.hc2,ylab="高度")
rect.hclust(fit.hc2,k=2,border = 3)

fit.hc3<-hclust(d=DisMatrix1,method="ward.D")    
cutree(fit.hc3,k=2)
pingjia(cutree(fit.hc3,k=2),c(rep(1,25),rep(2,29)))
plot(fit.hc3,ylab="高度")
rect.hclust(fit.hc3,k=2,border = 4)

##变量聚类法
DisMatrix5<-as.dist(cor(X1))
par(mfrow=c(1,2))
fit.hc5<-hclust(d=DisMatrix5,method="single")    
cutree(fit.hc5,k=2)
plot(fit.hc5,ylab="高度")

fit.hc6<-hclust(d=DisMatrix5,method="complete")    
cutree(fit.hc6,k=2)
plot(fit.hc6,ylab="高度")


#############################################################################K-means聚类
fit.km<-kmeans(x=X1,centers=2,iter.max = 50,nstart = 3)
fit.km$cluster
pingjia(fit.km$cluster,c(rep(1,25),rep(2,29)))

###############################################################################PAM聚类
library(cluster)
fit.pam<-pam(x=X1,k=2,do.swap=TRUE,stand=FALSE)
fit.pam$clustering
pingjia(fit.pam$clustering,c(rep(1,25),rep(2,29)))

###############################################################################EM聚类
#install.packages("mclust")
library("mclust") 
fit.em<-Mclust(data=X1,2)  
pingjia(fit.em$classification,c(rep(1,25),rep(2,29)))
summary(fit.em)
summary(fit.em,parameters=TRUE)   
plot(fit.em$classification,pch=fit.em$classification+14,col=fit.em$classification,ylab="类别编号",xlab="股票名称",main="传统数据的EM聚类",axes=FALSE)
par(las=2)
axis(1,at=1:dim(X)[1],labels=rownames(X1),cex.axis=0.6)
axis(2,at=1:3,labels=1:3,cex.axis=0.6)
box()
legend("topright",c("第一类","第二类","第三类"),pch=15:17,col=1:3,cex=0.6)

############################################################################主成分聚类
rm(list=ls())  
read.csv("新冠检测概念数据(20210219).csv",header=T)->DA
X<-DA[1:53,2:10]
d<-dist(X)
d <- dist(scale(X))
X<-scale(X) #标准化
########检验PCA是否适用
#install.packages("psych")
#install.packages("MASS")
library(psych)
cortest.bartlett(cor(X), n = 120)
# KMO Kaiser-Meyer-Olkin Measure of Sampling Adequacy
kmo = function( data ){
  library(MASS)
  X <- cor(as.matrix(data))
  iX <- ginv(X)
  S2 <- diag(diag((iX^-1)))
  AIS <- S2%*%iX%*%S2                      # anti-image covariance matrix
  IS <- X+AIS-2*S2                         # image covariance matrix
  Dai <- sqrt(diag(diag(AIS)))
  IR <- ginv(Dai)%*%IS%*%ginv(Dai)         # image correlation matrix
  AIR <- ginv(Dai)%*%AIS%*%ginv(Dai)       # anti-image correlation matrix
  a <- apply((AIR - diag(diag(AIR)))^2, 2, sum)
  AA <- sum(a)
  b <- apply((X - diag(nrow(X)))^2, 2, sum)
  BB <- sum(b)
  MSA <- b/(b+a)                        # indiv. measures of sampling adequacy
  
  AIR <- AIR-diag(nrow(AIR))+diag(MSA)  # Examine the anti-image of the
  # correlation matrix. That is the
  # negative of the partial correlations,
  # partialling out all other variables.
  
  kmo <- BB/(AA+BB)                     # overall KMO statistic
  
  # Reporting the conclusion
  if (kmo >= 0.00 && kmo < 0.50){
    test <- 'The KMO test yields a degree of common variance
    unacceptable for FA.'
    } else if (kmo >= 0.50 && kmo < 0.60){
      test <- 'The KMO test yields a degree of common variance miserable.'
    } else if (kmo >= 0.60 && kmo < 0.70){
      test <- 'The KMO test yields a degree of common variance mediocre.'
    } else if (kmo >= 0.70 && kmo < 0.80){
      test <- 'The KMO test yields a degree of common variance middling.'
    } else if (kmo >= 0.80 && kmo < 0.90){
      test <- 'The KMO test yields a degree of common variance meritorious.'
    } else {
      test <- 'The KMO test yields a degree of common variance marvelous.'
    }

    ans <- list(  overall = kmo,
                  report = test,
                  individual = MSA,
                  AIS = AIS,
                  AIR = AIR )
    return(ans)

}    # end of kmo()
kmo(X)  
#########################################

pr<-princomp(X, cor=TRUE)
summary(pr, loadings=T)
#对数据做主成分分析

X<-as.matrix(X)
XX<-t(X)%*%X
eigen(XX)#求矩阵的特征值
X%*%eigen(XX)$vectors

#主成分分析
pr<-princomp(X, cor=TRUE)
summary(pr, loadings=T)
pre<-predict(pr);pre
screeplot(pr,npcs=10, type=c("lines"))#主成分碎石图
mtext("主成分", side=1,line=2)
mtext("特征值", side=2,line=2)
par(mfrow=c(1,2))
biplot(pr, xlab="第1主成分",ylab="第2主成分")#缺省的是第一、二主成分散点图
biplot(pr,choices=3:4,scale=1,xlab="第3主成分",ylab="第4主成分")
#biplot(pr,choices=5:6,scale=1)

#加权主成分聚类
data<-pre
PCAdata<-as.matrix(data[,1:6])
#各主成分方差贡献率作为权重
tezhz<-eigen(XX)$values
w<-tezhz[1:6]/sum(tezhz[1:6])
f<-PCAdata%*%w
F<-as.data.frame(cbind(PCAdata,f))
F<-cbind(DA[1:53,1],F)
datax<-F

#聚类
X<-data
rownames(X)<-datax[,1]
X<-X[,1:6]
data<-scale(X)

distance <- dist(data,method ="canberra")  #lance距离
data.hc <- hclust(distance) #最长距离法
plot(data.hc, hang = -1,ylab="高度") #绘画系谱图 
P <- rect.hclust(data.hc, k =3) #分为3类

data.hc <- hclust(distance,method="single") #最短距离法
plot(data.hc, hang = -1) #绘画系谱图 
P <- rect.hclust(data.hc, k = 3) #分为3类

data.hc <- hclust(distance,method="average") #类平均法
plot(data.hc, hang = -1) #绘画系谱图 
P <- rect.hclust(data.hc, k = 3) #分为3类

data.hc <- hclust(distance,method="centroid") #重心法
plot(data.hc, hang = -1) #绘画系谱图 
P <- rect.hclust(data.hc, k = 3) #分为3类

data.hc <- hclust(distance,method="median") #中间距离法
plot(data.hc, hang = -1) #绘画系谱图 
P <- rect.hclust(data.hc, k = 3) #分为3类

data.hc <- hclust(distance,method="ward.D") #离差平方和法
plot(data.hc, hang = -1) #绘画系谱图 
P <- rect.hclust(data.hc, k = 3) #分为3类

data.hc <- hclust(distance,method="mcquitty") #相似法
plot(data.hc, hang = -1) #绘画系谱图
P <- rect.hclust(data.hc, k = 3) #分为3类

######################################################各种图的画法
#聚类热力图
x_matrix <- data.matrix(data)#将数据框转化为矩阵
rc <- rainbow(nrow(X), start = 0, end = .3)#设置颜色
heatmap(
  x_matrix, col = rainbow(56), #选择颜色
  scale = "column",#标准化数据
  margins = c(5,10),RowSideColors = rc,#给每支股票赋予对应的颜色
  main = "",cexCol=1,Colv=NA#只在列上聚类
)


#plot.phylo函数的4种不同类型的聚类树形图
hc = hclust(dist(data))
#install.packages("ape")
library(ape)

#1.碎屑图
plot(as.phylo(hc), type = "cladogram", cex = 1.0, label.offset = 0.1)

#2.无根图
plot(as.phylo(hc), type = "unrooted",cex=1.0, label.offset = 0.1)

#3. 辐射状图
plot(as.phylo(hc), type = "radial")

#4. 扇形图
d <- dist(X, method = "euclidean")
# Hierarchical clustering using Ward's method
hc <- hclust(d, method = "complete")
# Plot the obtained dendrogram
mypal = c("#556270", "#4ECDC4", "#1B676B", "#FF6B6B", "#C44D58")
clus4 = cutree(hc,5)
# Size reflects miles per gallon
plot(as.phylo(hc), type = "fan", tip.color = mypal[clus4], label.offset =0.1, 
    cex =0.7,col = "red")