Example

dist函数默认使用欧式距离进行聚类,类似其他的距离定义可以用method来定义:

• method
  


• the distance measure to be used. This must be one of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski". Any unambiguous substring can be given.

hclust函数对数据进行聚类,连接准则使用“ave”

• 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).

dist(c(1:5))

##   1 2 3 4
## 2 1
## 3 2 1
## 4 3 2 1
## 5 4 3 2 1

## default:method = "euclidean"
## Show the distance of the vetor
hc <- hclust(dist(USArrests), "ave")
## hc is the answer of hierarchical clustering

plot(hc)

啰嗦几句,可以看到聚类的结果是一个树状图,也就是说,具体聚成几类不是层次聚类决定的,基本上是需要通过一个剪枝的过程来完成,比如,想要两类就划一条水平的线,刚好和树型图相交两次,化为左右两棵树,便是需要的两类.

K-means

James MacQueen,1967[2]最早提出kmeans的算法. kmeans的思路是假设每个类都有一个中心点,而计算的目标是最小化各个数据到它所在类中心点的距离:

[W(C)=\sum^k_{i=1}\sum_{x \in S_i}||x-\mu_i||^2]

基本算法:

1. 初始化k个中心点({ m_1,…,m_K})如果已经分为几类,则计算这几类的中心点


2. 对于一个给定的中心点,计算将所有数据分类到离它最近的那个中心点: [C(i)=argmin_k||x_i-x_k||^2]

重复1和2 直到类不发生任何改变

kmeans的算法实际上会耗费比较多的时间,具体的复杂度为(O(n^{dk+1}\log{n}))(复杂度的部分我都未验证,来源这里[3]),而现在常用的kmeans算法使用的是另一套较为高效的算法:Lloyds algorithm(O(n k d i)):

在给定一系列聚类初值的情况下,进行如下的步骤:

1. Assignment step: [S_i^{(t)} = \big { x_p : \big | x_p - m^{(t)}_i \big |^2 \le \big | x_p - m^{(t)}_j \big |^2 \ \forall j, 1 \le j \le k \big},] 因此(x_p)应当是归于(S_i^{(t)})类.




2. Update step:

   
   

   [m^{(t+1)}_i = \frac{1}{|S^{(t)}i|} \sum{x_j \in S^{(t)}_i} x_j]

Example

require(graphics)
set.seed(123)
# a 2-dimensional example
x <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
           matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
colnames(x) <- c("x", "y")
(cl <- kmeans(x, 2))

## K-means clustering with 2 clusters of sizes 50, 50
##
## Cluster means:
##         x       y
## 1 0.92383 1.01164
## 2 0.01032 0.04392
##
## Clustering vector:
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [36] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 8.139 7.396
##  (between_SS / total_SS =  74.0 %)
##
## Available components:
##
## [1] "cluster"      "centers"      "totss"        "withinss"
## [5] "tot.withinss" "betweenss"    "size"         "iter"
## [9] "ifault"

library(ggplot2)
cluster=as.factor(cl$cluster)
qplot(x[,1],x[,2], col = cluster)+
geom_point(aes(cl$centers[,1],cl$centers[,2]),size=4, col = 1:2, pch = 8, cex = 2)+
  labs(x = "X", y = "Y",title="kmeans")

Example

恩恩,mclust包里面的clustCombi函数提供了合适的渠道来快速GMM~

哈哈~省心..

library(mclust)

## Warning: package ‘mclust’ was built under R version 3.1.2

## Package ‘mclust’ version 4.4
## Type ‘citation("mclust")’ for citing this R package in publications.

output <- clustCombi(x)
output

##
##  EM/BIC Solution
##  —————
##
## Number of components: 2
## Model name: EII
##
## Component num.1:
##      proportion: 0.50
##      mean: 0.01  0.05
## Component num.2:
##      proportion: 0.50
##      mean: 0.92  1.01
##
##  Combining steps
##  —————
##
##   Step | Classes combined at this step | Classes labels after this step
## ——-|——————————-|——————————-
##    0   |              —              | 1   2
##    1   |            1  &  2            | 1
##
##  Classification for K classes: output$classification[[K]]
##  Combining matrix (K classes -> (K-1) classes): output$combiM[[K]]

plot(output, x,what="classification",reg=2) 

DBSCAN

DBSCAN的定义是基于density reachability(密度可达性)实现的,p和q两个点称为直接密度可达是指p和q之间的距离小于一个给定的距离(\epsilon),称p,q是密度可达是指:存在一个序列 (p_1,\ldots,p_n)使得(p_1 = p)并且(p_n = q) 其中每一个 (p_{i+1})都和(p_i)相连通.思路很像基于连通的聚类是不是?

对于一个聚类,具有以下两个性质(DBSCAN的游戏规则):

1. 一个聚类内的所有点都是密度可达的.


2. 如果一个点对于聚类中每个点都是密度可达的,那这个点就是这个类中的一个.

B和C是密度可达的:

DBSCAN需要两个参数来开始执行运算,一个是eps代表(\epsilon) 最短距离,MinPts代表一个聚类中最小点数.伪代码算法如下:

## pseudocode
DBSCAN(D, eps, MinPts)
   C = 0
   for each unvisited point P in dataset D
      mark P as visited
      NeighborPts = regionQuery(P, eps)
      if sizeof(NeighborPts) < MinPts
         mark P as NOISE
      else
         C = next cluster
         expandCluster(P, NeighborPts, C, eps, MinPts)

expandCluster(P, NeighborPts, C, eps, MinPts)
   add P to cluster C
   for each point P’ in NeighborPts
      if P’ is not visited
         mark P’ as visited
         NeighborPts’ = regionQuery(P’, eps)
         if sizeof(NeighborPts’) >= MinPts
            NeighborPts = NeighborPts joined with NeighborPts’
      if P’ is not yet member of any cluster
         add P’ to cluster C

regionQuery(P, eps)
   return all points within P’s eps-neighborhood (including P)

OPTICS

OPTICS和DBSCAN聚类相同,都是基于密度的聚类,但是,OPTICS的好处在于可以处理不同密度的类,结果有点像基于连通性的聚类,不过还是有些区别的. 上段伪代码:

 OPTICS(DB, eps, MinPts)
    for each point p of DB
       p.reachability-distance = UNDEFINED
    for each unprocessed point p of DB
       N = getNeighbors(p, eps)
       mark p as processed
       output p to the ordered list
       if (core-distance(p, eps, Minpts) != UNDEFINED)
          Seeds = empty priority queue
          update(N, p, Seeds, eps, Minpts)
          for each next q in Seeds
             N’ = getNeighbors(q, eps)
             mark q as processed
             output q to the ordered list
             if (core-distance(q, eps, Minpts) != UNDEFINED)
                update(N’, q, Seeds, eps, Minpts)

In update(), the priority queue Seeds is updated with the \varepsilon-neighborhood of p and q, respectively:

 update(N, p, Seeds, eps, Minpts)
    coredist = core-distance(p, eps, MinPts)
    for each o in N
       if (o is not processed)
          new-reach-dist = max(coredist, dist(p,o))
          if (o.reachability-distance == UNDEFINED) // o is not in Seeds
              o.reachability-distance = new-reach-dist
              Seeds.insert(o, new-reach-dist)
          else               // o in Seeds, check for improvement
              if (new-reach-dist < o.reachability-distance)
                 o.reachability-distance = new-reach-dist
                 Seeds.move-up(o, new-reach-dist)

Example

看起来EM算法相对来说是一套复杂的方法,我自己做了另外一套的估计方法,效果还不错,因为涉及到自己的算法,就不再展示代码了,直观的做一个图来说明我所做的结果:

两个线性模型生成的数据,然后用一个算法估计两条直线的参数…

Spectral clustering

谱聚类

谱聚类最早来自于图像的识别领域,最早提出谱聚类的人我到现在也没有考证到,不过目测最基本的算法出来是要被吐槽的(因为实际上没太多理论价值),而在后来,谱聚类演化到图像识别领域,并且和图论的算法相结合,得到了不错的进展,反过来,在统计和机器学习等领域大放异彩,值得参考的是数据之城的博文.和A Tutorial on Spectral Clustering[10].

首先来概括下最直观,最简单的谱聚类:主成分+Kmeans.说完了…所以这就是为什么我觉得最早谱聚类的雏形并没有很多引用,甚至没有被提到,当然,没有被提到分两种情况,第一种是因为论文实在是太厉害,所有人都知道里面的知识,所以不用提到,比如AIC准则的论文,多少论文里面都会用AIC,看看有几个会提到1957年的那篇论文,第二种呢,是因为论文所涉及的理论和方法并不是很fancy,再后面推广什么不提原文也不会有太多心理压力,我猜谱聚类是因为后者,来看下谱聚类的算法:

1. 给出(N)个数据点两两之间的相似性。也就是一个(NN)的相似性矩阵(A)，(A(i,j))代表(i)和(j)两个数据点的相似度，数值越大则表示越相似。注意(A(i,j)=A(j,i)),(A(i,i)=0).




2. 计算矩阵(D)，使它的对角元是(A)矩阵的对应的那一列(或行)的值之和，其余地方为0。也就是使得 [D(i,i)=\sum_jA(i,j)]




3. 令(B=D-A)




4. 求(B)矩阵的前(k)个特征值和特征向量，将B矩阵投影到一个(k)维空间。第(i)个特征值的第j个值，就表示第j个数据点在k维空间中第i维的投影。就是说如果把k个特征向量并成一个(Nk)的矩阵，则每一行代表一个数据点在(k)维空间的坐标。




5. 根据每个数据点的(k)维空间坐标，使用K-means或者其它聚类算法在k维空间对数据进行聚类。

A <- as.matrix(dist(x,diag = T))
D <- diag(colSums(A))
B <- D-A
eig <- eigen(B)
cum <- cumsum(eig$values)/sum(eig$values)
k = which(cum > 0.8)[1]
final_data <- B%%eig$vectors[,1:k]
cl <- kmeans(final_data, 3)
cluster=as.factor(cl$cluster)
qplot(x[,1],x[,2], col = cluster)+
  labs(x = "X", y = "Y",title="Spectral Clustering")

Chameleon

变色龙聚类

别看名字如此帅气,事实上,Chameleon的名字是这样的:Chameleon:Hierarchical Clustering Using Dynamic Modeling[11].明尼苏达大学的三位统计学家(George Karypis,Eui-Hong (Sam) Han,Vipin Kumar)在捣鼓层次聚类的结果,不过不得不说,相对于谱聚类,展现出的效果靠谱的多.好吧,虽然说出身是层次聚类(基于连通性的聚类),但是做出的效果绝对是前沿的高端结果.

这就是论文在1999年发表出来的时候上的效果图(MD我何时才能有出头之日啊!)

说一下我对CC(Chameleon Clustering)的理解吧,CC和层次聚类类似,是一个基于连通性的聚类方法,主要的流程是:

• 先将数据聚类为一个稀疏的图(类似Knn的方法).


• 将数据分割为一个个小块.


• 再把小块根据连通性和分割平均强度来进行划分.

Canopy

Canopy是一个大数据下的聚类框架,最早接触canopy是在mahout的机器学习库中看到的.canopy是一个计算复杂度在(O(nk))的聚类系统,与之相对的是Canopy较低的计算精度,二者是相互取舍得到这样的结果.

Canopy来自于Andrew McCallum, Kamal Nigam and Lyle Ungar in 2000.[9],Canopy的意思是华盖,苍穹(苍穹聚类是不是一个很有厚重感的名字?),基本思路是:

“对于每一个传统聚类,都会存在一个Canopy使得在这个聚类内所有观测值都它里面”

“For every traditional cluster, there exists a canopysuch that all elements of the cluster are in thecanopy.”

从这里就可以看出Canopy的创造者是彻头彻尾的计算机背景(三位作者是卡内基梅隆大学Carnegie Mellon University的CS系),数学背景会选择说存在一个neighbourhood(邻域),而统计学家可能会选择说存在一个region(区间)来包含观测值,是没有Canopy这样的名字来的文艺范十足.

另外要注意一点的是,Canopy并不一定是一个,也可以是一族,上图的虚线和实线都是代表着一个Canopy,而在他们之间,同时存在多个Cannopies,如果细分,虚线和实线的半径就是之后要定义的两个阈值,松阈值(loose distance)和紧阈值(tight distance).

Examples

这里写一个小程序来计算各个给定数据,t1和t2的时候计算canopy的各个中心

## Calculate the center of each canopy
canopy_cluster <- function(data=x,t1,t2){
  n = dim(data)[1]
  data = data[sample(n),]
  ## Sample the data
  distance = as.matrix(dist(data,diag = T))
  ## Calculate the distance matrix
  judge=rep(T,n)
  index=numeric(0)
  for ( i in 1:n){
    if(judge[i]==F) next
    index=c(index,i)
    pd = which(distance[i,]<t2)
    judge[pd] = rep(F,length(pd))
}

  data[index,]

}

## Draw the figure of the canopy
canopy_plot=function(data,t1,t2,original_data=x){
  t = seq(0,2pi,length.out = 100)
  n = dim(data)[1]
  t1_data = data.frame(x = 0,y = 0,num = 0)
  t2_data = data.frame(x = 0,y = 0,num = 0)
  for ( i in 1:n){
    t1_data[((i-1)100+1):(i100),1] = data[i,1]+t1cos(t)
    t1_data[((i-1)100+1):(i100),2] = data[i,2]+t1sin(t)
    t1_data[((i-1)100+1):(i100),3] = rep(i,100)
     t2_data[((i-1)100+1):(i100),1] = data[i,1]+t2cos(t)
     t2_data[((i-1)100+1):(i100),2] = data[i,2]+t2sin(t)
     t2_data[((i-1)100+1):(i*100),3] = rep(i,100)
  }

  p=qplot(x[,1],x[,2])+
geom_point(aes(data[,1],data[,2]),size=4, col = 1:n, pch = 8, cex = 2)+
  labs(x = "X", y = "Y",title="Canopy")+
  geom_path(data=t1_data,aes(x=x,y=y,group=num),linetype=2)+
  geom_path(data=t2_data,aes(x=x,y=y,group=num),linetype=1)

return(p)
}

set.seed(123)
data=canopy_cluster(x,t2=0.5)
canopy_plot(data,t1=1,t2=0.5)

data=canopy_cluster(x,t2=1)
canopy_plot(data,t1=1.5,t2=1)

data=canopy_cluster(x,t2=1.5)
canopy_plot(data,t1=1.7,t2=1.5)