Cluster analysis or clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense or another) to each other than to those in other groups (clusters). It is a main task of exploratory data mining, and a common technique for statistical data analysis, used in many fields, including machine learning, pattern recognition, image analysis, information retrieval, bioinformatics, data compression, and computer graphics.

——from wikipedia

聚类分析起源于分类学，在考古的分类学中，人们主要依靠经验和专业知识来实现分类。随着生产技术和科学的发展，人类的认识不断加深，分类越来越细，要求也越来越高，有时光凭经验和专业知识是不能进行确切分类的，往往需要定性和定量分析结合起来去分类，于是数学工具逐渐被引进分类学中，形成了数值分类学。后来随着多元分析的引进，聚类分析又逐渐从数值分类学中分离出来而形成一个相对独立的分支。

值得提出的是将聚类分析和其它方法联合起来使用，如判别分析、主成分分析、回归分析等往往效果更好。

聚类方法的分类

系统聚类法(层次化聚类算法)
划分式聚类算法
- k-means聚类法
基于密度和网格的聚类算法
机器学习的聚类算法（SVM和神经网络）
其它聚类算法：
- 聚类预报法
- 模糊聚类法

聚类算法分类图

系统聚类法：

层次聚类又称为系统聚类(Hierarchical Clustering)是一种常见的聚类算法,层次聚类算法将数据分层建立簇,形成一棵树。如果按自底向上进行层次分解,则称为凝聚(agglomerative)的层次聚类,而按自顶向下的进行层次分解,则称为分裂法(divisive)层次聚类。这种算法的基本思路是首先将所有对象看成独立的个体类,通过计算类间的距离来选择最小距离的两个类合并成一个新类,再重新计算新类和其它类之间的距离,选择最小距离的两个类合并,依次迭代合并直到没有合并为止。

其中在计算个体间的距离时：定义“距离”的统计量包括了欧氏距离 (euclidean)、马氏距离(manhattan)、两项距离(binary)、明氏距离(minkowski；。还包括相关系数和夹角余弦。前者用于样品间的分类，后者用于指标间的分类！

将每个样本单独作为一类，然后将不同类之间距离最近的进行合并，合并后重新计算类间距离。这个过程一直持续到将所有样本归为一类为止。在计算类间距离时比较传统的有六种不同的方法，分别是最短距离法(single)、最长距离法(complete ,默认)、类平均法、重心法(centroid)、中间距离法(median)、离差平方和法(ward)。还有加权平均法(Weighted pair group method average,WPGMA),无加权平均法(Unweighted pair group method average,UPGMA),无加权中心法(Unweighted pair group method centroid,UPGMC)等。当然还有其他的新的方法。

层次聚合算法的计算复杂性为O(n^2),适合于小型数据集的分类. 层次聚类方法有两个主要的缺陷:(1)难于确定聚类结果的类个数;(2)聚类过程中合并或是分裂均不可回复,因此影响了聚类的结果,好处是可以减少尝试的次数,加快计算过程。确定聚类的类个数是一个算法难题,目前还没有很好的方式来确定。

下面利用R来进行实践

个体、样本之间的距离

type function *dist()*

dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2) **下面是各参数的意义！**


x   
a numeric matrix, data frame or "dist" object.

method  
the distance measure to be used. This must be one of "euclidean（欧式距离）", "maximum", "manhattan（马氏距离）", "canberra（兰式距离，适用于值大于0）", "binary（0,1）" or "minkowski（明氏距离，最一般化的形式）". Any unambiguous substring can be given.

diag
logical value indicating whether the diagonal of the distance matrix should be printed by print.dist.

upper
logical value indicating whether the upper triangle of the distance matrix should be printed by print.dist.

p
The power（幂） of the Minkowski distance.当值为1时，是绝对距离；2的话就是欧式距离

m
An object with distance information to be converted to a "dist" object. For the default method, a "dist" object, or a matrix (of distances) or an object which can be coerced to such a matrix using as.matrix(). (Only the lower triangle of the matrix is used, the rest is ignored).

digits, justify
passed to format inside of print().

right, ...
further arguments, passed to other methods.

# see the examples given by the docu. 
#这个是对行，也就是对样品进行聚类
x <- matrix(rnorm(100), nrow = 5)
dist(x)

##          1        2        3        4
## 2 6.368985                           
## 3 5.046238 6.168687                  
## 4 5.365449 6.696279 5.013462         
## 5 6.908890 4.346960 6.201976 6.099895

t<-dist(x)
dist(x, diag = TRUE)#看diag参数是什么意思

##          1        2        3        4        5
## 1 0.000000                                    
## 2 6.368985 0.000000                           
## 3 5.046238 6.168687 0.000000                  
## 4 5.365449 6.696279 5.013462 0.000000         
## 5 6.908890 4.346960 6.201976 6.099895 0.000000

dist(x, upper = TRUE)#看upper参数是什么意思

##          1        2        3        4        5
## 1          6.368985 5.046238 5.365449 6.908890
## 2 6.368985          6.168687 6.696279 4.346960
## 3 5.046238 6.168687          5.013462 6.201976
## 4 5.365449 6.696279 5.013462          6.099895
## 5 6.908890 4.346960 6.201976 6.099895

dist(x,method = "manhattan")#默认是用欧式距离来计算；用马氏距离看看，再看看结果！

##          1        2        3        4
## 2 24.41463                           
## 3 18.42940 22.92726                  
## 4 19.03564 20.50912 20.55415         
## 5 24.98495 16.43086 22.22202 19.38073

类间的距离

可以利用hclust()函数和heapmap图来进行选择与计算

type function *hclust()*

hclust(d, method = "complete", members = NULL)

Arguments

d       a dissimilarity structure as produced by dist.
method      the agglomeration method to be used. This should be (an unambiguous abbreviation of) one of "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
members     NULL or a vector with length size of d. See the ‘Details’ section.

m <- as.matrix(dist(x))
d <- as.dist(m)
plot(hclust(d)) # to see a dendrogram of clustered variables,默认的类间聚类方法是compelete

plot(hclust(d,method = "average"))

plot(hclust(d,method = "mcquitty"))

plot(hclust(d,method = "median"))

plot(hclust(d,method = "centroid"))

plot(hclust(d,method = "single"))

plot(hclust(d,method = "ward.D"))

plot(hclust(d,method = "ward.D2"))

heatmap(x)

##　heatmap函数

type function *heatmap()*

Usage:

heatmap(x, Rowv = NULL, Colv = if(symm)"Rowv" else NULL,
        distfun = dist, hclustfun = hclust,
        reorderfun = function(d, w) reorder(d, w),
        add.expr, symm = FALSE, revC = identical(Colv, "Rowv"),
        scale = c("row", "column", "none"), na.rm = TRUE,
        margins = c(5, 5), ColSideColors, RowSideColors,
        cexRow = 0.2 + 1/log10(nr), cexCol = 0.2 + 1/log10(nc),
        labRow = NULL, labCol = NULL, main = NULL,
        xlab = NULL, ylab = NULL,
        keep.dendro = FALSE, verbose = getOption("verbose"), ...)

Arguments

x   
numeric matrix of the values to be plotted. 你需要画热图的一个矩阵
Rowv    决定是否对行进行聚类
determines if and how the row dendrogram should be computed and reordered. Either a dendrogram or a vector of values used to reorder the row dendrogram or NA to suppress any row dendrogram (and reordering) or by default, NULL, see ‘Details’ below.
Colv    决定是否对列进行聚类
determines if and how the column dendrogram should be reordered. Has the same options as the Rowv argument above and additionally when x is a square matrix, Colv = "Rowv" means that columns should be treated identically to the rows (and so if there is to be no row dendrogram there will not be a column one either).
distfun  相当于dist()函数中对个体的聚类
function used to compute the distance (dissimilarity) between both rows and columns. Defaults to dist.
hclustfun    相当于hclust函数中对类之间关系的调整
function used to compute the hierarchical clustering when Rowv or Colv are not dendrograms. Defaults to hclust. Should take as argument a result of distfun and return an object to which as.dendrogram can be applied.
scale    是否进行中心化
character indicating if the values should be centered and scaled in either the row direction or the column direction, or none. The default is "row" if symm false, and "none" otherwise.

下面这个是对指标，也就是对列进行聚类！

## Use correlations between variables "as distance" 

dd <- as.dist((1 - cor(USJudgeRatings))/2)
round(1000 * dd) # (prints more nicely)

##      CONT INTG DMNR DILG CFMG DECI PREP FAMI ORAL WRIT PHYS
## INTG  567                                                  
## DMNR  577   18                                             
## DILG  494   64   82                                        
## CFMG  432   93   93   21                                   
## DECI  457   99   98   22    9                              
## PREP  494   61   72   11   21   21                         
## FAMI  513   66   79   21   32   29    5                    
## ORAL  506   44   47   23   25   26    8    9               
## WRIT  522   46   53   20   29   27    7    5    3          
## PHYS  473  129  106   94   60   64   76   78   54   72     
## RTEN  517   31   28   35   36   38   25   29    9   16   47

plot(hclust(dd)) # to see a dendrogram of clustered variables

## example of binary and canberra distances.
x <- c(0, 0, 1, 1, 1, 1)
y <- c(1, 0, 1, 1, 0, 1)
dist(rbind(x, y), method = "binary")

##     x
## y 0.4

## answer 0.4 = 2/5
dist(rbind(x, y), method = "canberra")

##     x
## y 2.4

## answer 2 * (6/5)

再来看一个例子

下面我们用iris数据集来进行聚类分析，在R语言中所用到的函数为hclust。首先提取iris数据中的4个数值变量，然后计算其欧氏距离矩阵。然后将矩阵绘制热图，从图中可以看到颜色越深表示样本间距离越近，大致上可以区分出三到四个区块，其样本之间比较接近。

data=iris[,-5]
dist.e=dist(data,method='euclidean')
heatmap(as.matrix(dist.e),labRow = F, labCol = F)

然后使用hclust函数建立聚类模型，结果存在model1变量中，其中ward参数是将类间距离计算方法设置为离差平方和法。使用plot(model1)可以绘制出聚类树图。如果我们希望将类别设为3类，可以使用cutree函数提取每个样本所属的类别。

model1=hclust(dist.e,method='ward.D')
plclust(model1)

## Warning: 'plclust' is deprecated.
## Use 'plot' instead.
## See help("Deprecated")

plot(model1)

#rect.hclust(model1)  #用矩形画出分为3类的区域
result=cutree(model1,k=3)
plot(result)

为了显示聚类的效果，我们可以结合多维标度和聚类的结果。先将数据用MDS进行降维，然后以不同的的形状表示原本的分类，用不同的颜色来表示聚类的结果。可以看到setose品种聚类很成功，但有一些virginica品种的花被错误和virginica品种聚类到一起。

k-means聚类法

对于一个给定的 n 个数据对象的数据集,采用目标函数最小化的策略,初始时选择一定量的聚类中心或数据点,通过某种原则把数据划分到各个组中,每个组为一个簇。最典型的划分式聚类算法就是 k-means 算法和 k-medoids 算法,这两种算法的改进算法非常多,应用也很广泛。

优缺点：

k-means 算法实现非常简单,运算效率也非常的高,适合对大型数据集进行分析处理。缺点是聚类结果不能重复,聚类结果跟初始点的选择有很大的关系,且不能作用于非凸集的数据。k-means 算法对类球形且大小差别不大的类簇有很好的表现,但不能发现形状任意和大小差别很大的类簇,且聚类结果易受噪声数据影响。

算法基本原理：

kmeans的计算方法如下：

*1 随机选取k个中心点

*2 遍历所有数据，将每个数据划分到最近的中心点中

*3 计算每个聚类的平均值，并作为新的中心点

*4 重复2-3，直到这k个中线点不再变化（收敛了），或执行了足够多的迭代

时间复杂度：O(Ink*m)

空间复杂度：O(n*m)

其中m为每个元素字段个数，n为数据量，I为跌打个数。一般I,k,m均可认为是常量，所以时间和空间复杂度可以简化为O(n)，即线性的。

聚类结果评价：

轮廓系数

轮廓系数（Silhouette Coefficient）结合了聚类的凝聚度（Cohesion）和分离度（Separation），用于评估聚类的效果。该值处于-1~1之间，值越大，表示聚类效果越好。具体计算方法如下：

对于第i个元素x_i，计算x_i与其同一个簇内的所有其他元素距离的平均值，记作a_i，用于量化簇内的凝聚度。选取x_i外的一个簇b，计算x_i与b中所有点的平均距离，遍历所有其他簇，找到最近的这个平均距离,记作b_i，用于量化簇之间分离度。对于元素x_i，轮廓系数s_i = (b_i – a_i)/max(a_i,b_i) 计算所有x的轮廓系数，求出平均值即为当前聚类的整体轮廓系数从上面的公式，不难发现若s_i小于0，说明x_i与其簇内元素的平均距离小于最近的其他簇，表示聚类效果不好。如果a_i趋于0，或者b_i足够大，那么s_i趋近与1，说明聚类效果比较好。

K值的选取：

在实际应用中，由于Kmean一般作为数据预处理，或者用于辅助分类贴标签。所以k一般不会设置很大。可以通过枚举，令k从2到一个固定值如10，在每个k值上重复运行数次kmeans(避免局部最优解)，并计算当前k的平均轮廓系数，最后选取轮廓系数最大的值对应的k作为最终的集群数目。

library(fpc) # install.packages("fpc")
data(iris)
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

# 0-1 正规化数据
min.max.norm <- function(x){
    (x-min(x))/(max(x)-min(x))
}
raw.data <- iris[,1:4]
norm.data <- data.frame(sl = min.max.norm(raw.data[,1]),
                                        sw = min.max.norm(raw.data[,2]),
                                        pl = min.max.norm(raw.data[,3]),
                                        pw = min.max.norm(raw.data[,4]))
head(norm.data)

##           sl        sw         pl         pw
## 1 0.22222222 0.6250000 0.06779661 0.04166667
## 2 0.16666667 0.4166667 0.06779661 0.04166667
## 3 0.11111111 0.5000000 0.05084746 0.04166667
## 4 0.08333333 0.4583333 0.08474576 0.04166667
## 5 0.19444444 0.6666667 0.06779661 0.04166667
## 6 0.30555556 0.7916667 0.11864407 0.12500000

对iris的4个feature做数据正规化，每个feature均是花的某个不为的尺寸。

# k取2到8，评估K
K <- 2:8
round <- 30 # 每次迭代30次，避免局部最优
rst <- sapply(K, function(i){
    print(paste("K=",i))
    mean(sapply(1:round,function(r){
        print(paste("Round",r))
        result <- kmeans(norm.data, i)
        stats <- cluster.stats(dist(norm.data), result$cluster)
        stats$avg.silwidth
    }))
})

## [1] "K= 2"
## [1] "Round 1"
## [1] "Round 30"

plot(K,rst,type='l',main='轮廓系数与K的关系', ylab='轮廓系数')

评估k，由于一般K不会太大，太大了也不易于理解，所以遍历K为2到8。由于kmeans具有一定随机性，并不是每次都收敛到全局最小，所以针对每一个k值，重复执行30次，取并计算轮廓系数，最终取平均作为最终评价标准，可以看到如上的示意图，

我们看到当k=2时，有最大的轮廓系数，虽然实际上有三个种类。

下面进行降纬度观察: 多维标度分析(MDS)是一种将多维空间的研究对象简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。

设想一下如果我们在欧氏空间中已知一些点的座标，由此可以求出欧氏距离。那么反过来，已知距离应该也能得到这些点之间的关系。这种距离可以是古典的欧氏距离，也可以是广义上的“距离”。MDS就是在尽量保持这种高维度“距离”的同时，将数据在低维度上展现出来。从这种意义上来讲，主成分分析也是多维标度分析的一个特例。

MDS与PCA的区别：DMS是对高纬度的距离进行展示,而PCA是对多个变量进行降维进行可视化，多用来进行数据重复性的查看！

#cmdscale是传统的
old.par <- par(mfrow = c(1,2))
k = 2 # 根据上面的评估 k=2最优
clu <- kmeans(norm.data,k)
mds = cmdscale(dist(norm.data,method="euclidean"),eig=T)
sum(abs(mds$eig[1:2]))/sum(abs(mds$eig))

## [1] 0.9588785

sum((mds$eig[1:2])^2)/sum((mds$eig)^2)

## [1] 0.9982746

plot(mds$points, col=clu$cluster, main='kmeans聚类 k=2', pch = 19)
plot(mds$points, col=iris$Species, main='原始聚类', pch = 19)

par(old.par)

参考内容：

1. 网页:http://www.cnblogs.com/bourneli/p/3645049.html
1. 《聚类算法研究》，孙吉贵等
1. 《聚类分析及其应用研究》唐东明博士论文

study the cluster analysis!

cluster

YuDalang

2016年6月13日