多元统计分析r代码-2

#图2.1马氏距离例
rm(list=ls())
x=seq(from=-5,to=10,by=0.01)
y1=dnorm(x,0,1) #计算密度
y2=dnorm(x,4,2) #计算密度
plot(x,y1,type="l",xlab="", ylab="",axes=F,xlim=c(-5,10), ylim=c(0,0.4))
lines(x,y2,type="l")
axis(2)
rect(-5.8, 0, 10, 0.4)
rect(0, 0, 0, 0.4)
rect(1.658, 0, 1.658, 0.1)
rect(4, 0, 4, 0.2)
text(1.658,-0.01,expression(A),cex=1)
text(0,-0.008,expression(0),cex=0.9)
text(4,-0.01,expression(4),cex=1)


#例2.1
x<-c(1,1,3,4,2,0,1,2,3,5); dim(x)<-c(5,2)
d<-dist(x)
hc1<-hclust(d, "single"); hc2<-hclust(d, "complete")
hc3<-hclust(d, "average"); hc4<-hclust(d, "ward")

opar <- par(mfrow = c(2, 2))
plot(hc1,hang=-1); plot(hc2,hang=-1)
plot(hc3,hang=-1); plot(hc4,hang=-1)
par(opar)

dend1<-as.dendrogram(hc1)
opar <- par(mfrow = c(2, 2),mar = c(4,3,1,2))
plot(dend1)
plot(dend1, nodePar=list(pch = c(1,NA), cex=0.8, lab.cex = 0.8),
type = "t", center=TRUE)
plot(dend1, edgePar=list(col = 1:2, lty = 2:3), dLeaf=1, edge.root = TRUE)
plot(dend1, nodePar=list(pch = 2:1,cex=.4*2:1, col = 2:3), horiz=TRUE)
par(opar)


#例2.2
## 输入相关矩阵.
x<- c(1.000, 0.846, 0.805, 0.859, 0.473, 0.398, 0.301, 0.382,
0.846, 1.000, 0.881, 0.826, 0.376, 0.326, 0.277, 0.277,
0.805, 0.881, 1.000, 0.801, 0.380, 0.319, 0.237, 0.345,
0.859, 0.826, 0.801, 1.000, 0.436, 0.329, 0.327, 0.365,
0.473, 0.376, 0.380, 0.436, 1.000, 0.762, 0.730, 0.629,
0.398, 0.326, 0.319, 0.329, 0.762, 1.000, 0.583, 0.577,
0.301, 0.277, 0.237, 0.327, 0.730, 0.583, 1.000, 0.539,
0.382, 0.415, 0.345, 0.365, 0.629, 0.577, 0.539, 1.000)
names<-c("身高 x1", "手臂长 x2", "上肢长 x3", "下肢长 x4", "体重 x5",
"颈围 x6", "胸围 x7", "胸宽 x8")
r<-matrix(x, nrow=8, dimnames=list(names, names))
## 作系统聚类分析,
## as.dist()的作用是将普通矩阵转化为聚类分析用的距离结构.
d<-as.dist(1-r); hc<-hclust(d); dend<-as.dendrogram(hc)
## 写一段小程序, 其目的是在绘图命令中调用它, 使谱系图更好看.
nP<-list(col=3:2, cex=c(2.0, 0.75), pch= 21:22,
bg= c("light blue", "pink"),
lab.cex = 1.0, lab.col = "tomato")
addE <- function(n){
if(!is.leaf(n)) {
attr(n, "edgePar") <- list(p.col="plum")
attr(n, "edgetext") <- paste(attr(n,"members"),"members")
}
n
}
## 画出谱系图.
op<-par(mfrow=c(1,1), mar=c(4,3,0.5,0))
de <- dendrapply(dend, addE); plot(de, nodePar= nP)
par(op)

#分成三类的程序和计算结果
plclust(hc, hang=-1); re<-rect.hclust(hc, k=3)


#例2.3续例1.2，海南板块的股票数据
read.csv("hn-2.csv",header=T)->DA
r<-DA[,2:12]
https://www.wendangku.net/doc/ba16977641.html,s<-DA[,1]
X<-data.frame(r, https://www.wendangku.net/doc/ba16977641.html,s)

HN<-dist(scale(r))
hc1<-hclust(HN, "complete")
hc2<-hclust(HN, "average")
hc3<-hclust(HN, "centroid")
hc4<-hclust(HN, "ward")

opar<-par(mfrow=c(2,2), mar=c(5.2,4,0,0),cex=0.8)
plclust(hc1,hang=-1)
re1<-rect.hclust(hc1,k=5,border="red")
plclust(hc2,h

ang=-1)
re2<-rect.hclust(hc2,k=5,border="red")
plclust(hc3,hang=-1)
re3<-rect.hclust(hc3,k=5,border="red")
plclust(hc4,hang=-1)
re4<-rect.hclust(hc4,k=5,border="red")
par(opar)

#例2.4动态聚类
km <- kmeans(HN, 5, nstart = 20); km
sort(km$cluster)