kmeans聚類

本文轉載自查看原文 2018-06-30 17:59 1803 Kmeans

聚類算法介紹

k-means算法介紹

k-means聚類是最初來自於信號處理的一種矢量量化方法，現被廣泛應用於數據挖掘。k-means聚類的目的是將n個觀測值划分為k個類，使每個類中的觀測值距離該類的中心（類均值）比距離其他類中心都近。

k-means聚類的一個最大的問題是計算困難，然而，常用的啟發式算法能夠很快收斂到局部最優解。這通常與高斯分布的期望最大化算法相似，這兩種算法都采用迭代求精的方法。此外，它們都使用聚類中心來對數據進行建模

k-means算法的提出與發展

詹姆斯·麥奎恩（James MacQueen）1967年第一次使用這個術語“k-means”，雖然這個想法可以追溯到1957年的雨果·斯坦豪斯（Hugo Steinhaus）。
標准算法首先由Stuart Lloyd在1957年提出，作為脈沖編碼調制技術，盡管直到1982年才發布在貝爾實驗室以外。
1965年，E. W. Forgy發表了基本相同的方法，這就是為什么它有時被稱為Lloyd-Forgy。

k-means算法的優勢適應問題

k-means算法優點

是解決聚類問題的一種經典算法，簡單、快速。
對處理大數據集，該算法是相對可伸縮和高效率的。因為它的復雜度是$O(nkt)$, 其中, n 是所有對象的數目, k 是簇的數目, t 是迭代的次數。通常$k\ll n$且$t\ll n$ 。
當結果簇是密集的，而簇與簇之間區別明顯時, 它的效果較好。

k-means算法缺點

在簇的平均值被定義的情況下才能使用，這對於處理符號屬性的數據不適用。
必須事先給出k（要生成的簇的數目），而且對初值敏感，對於不同的初始值，可能會導致不同結果。
它對於“躁聲”和孤立點數據是敏感的，少量的該類數據能夠對平均值產生極大的影響。

k-means算法的思想介紹

(1) 選定某種距離作為數據樣本件的相似性度量

由於k-means聚類算法不適合處理離散型數據，因此在計算個樣本距離時，可以根據實際需要選擇歐氏距離、曼哈頓距離或者明可夫斯基距離中的一種來作為算法的相似性度量。

假設給定的數據集$X=\{x_m|m=1,2,...,total\}$,$X$中的樣本用d個屬性$A_1,A_2,...,A_d$來表示，並且d個描述屬性都是連續型數據。數據樣本$x_i=(x_{i1},x_{i2},x_{id}),x_j=(x_{j1},x_{j2}...,x_{jd})$，其中，$x_{i1},x_{i2},x_{id}$和$x_{j1},x_{j2}...,x_{jd}$分別是樣本$x_i$和$x_j$對應的d個描述屬性$A_1,A_2,...,A_d$的具體取值。樣本$x_i$和$x_j$之間的相似度通常用他們之間的距離d($x_i,x_j$)來表示，距離越小，樣本$x_i$和$x_j$越相似，差異度越小。距離越大，樣本$x_i$和$x_j$越不相似，差異越大。
歐氏距離公式如下：

\[d(x_i,x_j)=\sqrt{\sum_{k=1}^d(x_{ik}-x_{jk})^2} \]

曼哈頓距離如下：

\[d(x_i,x_j)=\sum_{k=1}^d|x_{ik}-x_{jk}| \]

明可夫斯基距離如下：

\[d(x_i,x_j)=\sqrt[p]{\sum_{k=1}^d|x_{ik}-x_{jk}|^p} \]

當$p=1$時，明氏距離即為曼哈頓距離，當$p=2$時，明氏距離即為歐氏距離，當$p=\infty$時，明氏距離即為切比雪夫距離。

(2) 選擇評價聚類性能的准則函數

k-means聚類算法使用誤差平方和准則函數來評價聚類性能。給定數據集X，其中只包含描述屬性，不包含類別屬性。假設X包含k個聚類子集$X_1,X2,...,X_k$,各個聚類子集中的樣本量分別為$n_1,n_2,...,n_k$，各個聚類子集均值代表點分別為$m_1,m_2,...,m_k$，則誤差平方和准則函數公式為：

\[E=\sum_{i=1}^k \sum_{p\in X_i}(p-m_i)^2 \]

(3) 相似度的計算根據一個簇中對象的平均值來進行

將所有對象隨機分配到k個非空的簇中。
計算每個簇的平均值，並用該平均值代表相應的簇。
根據每個對象與各個簇中心的距離，分配給最近的簇。
然后轉2，重新計算每個簇的均值。這個過程不斷重復知道滿足某個准則函數為止。

k-means實現流程：分步驟寫

k-means算法2個核心問題，一是度量記錄之間的相關性的計算公式，此處采用歐氏距離。一是更新簇內質心的方法，此處用平均值法，即means。
此時的輸入數據為簇的數目k和包含n個對象的數據庫——通常在軟件中用數據框或者矩陣表示。輸出k個簇，使平方誤差准則最小。
下面為實現k-means聚類的Python代碼

# (1)選擇初始簇中心。
# (2)對剩余的每個對象，根據其與各個簇中心的距離，將它賦給最近的簇。
# (3)計算新的簇中心。
# (4)重復(2)和(3)，直至准則函數不再明顯變小為止。
from numpy import *

#定義加載數據的函數。如果數據以文本形式存儲在磁盤內，可以用此函數讀取
def loadDataSet(fileName)
dataMat = []
fr = open(fileName)
for line in fr.readlines()
curLine = line.strip().split('\t')
fltLine = list(map(float,curLine))
dataMat.append(list(fltLine))
return dataMat

#該函數計算兩個向量的距離，即歐氏距離
def distEclud(vecA,vecB)
return sqrt(sum(power(vecA - vecB,2)))

#曼哈頓距離
def Manhattan(vecA,vecB)
return sum(abs(vecA-vecB))

#明考夫斯基距離
def MinKowski(vecA,vecB,p)
return power(sum(power(abs(vecA-vecB),p)),1/p)

#此處為構造質心，而不是從數據集中隨機選擇k個樣本點作為質心，
#也是比較合理的方法
def randCent(dataSet,k)
    n = shape(dataSet)[1]  #n為dataSet的列數
    centroids = mat(zeros((k,n)))  #構造k行n列的矩陣，就是k個質心的坐標
    for j in range(n)  #
        minJ = min(dataSet[,j])  #該列數據的最小值
        maxJ = max(dataSet[,j])  #該列數據的最大值
        rangeJ = float(maxJ-minJ)  #全距
        #隨機生成k個數值，介於minJ和maxJ之間，填充J列
        centroids[,j] = minJ + rangeJ * random.rand(k,1)  
      
    return centroids
   
#定義kMeans函數 
def kMeans(dataSet,k,distMeas=distEclud,createCent=randCent)
    m = shape(dataSet)[0]  #m為原始數據的行數
    #clusterAssment包含兩列，一列記錄簇索引值，一列存儲誤差
    clusterAssment = mat(zeros((m,2)))  
    centroids = createCent(dataSet,k)  #構造初始的質心
    clusterChanged = True  #控制變量
    while clusterChanged  #當控制變量為真時，執行下述循環
        clusterChanged = False  
        for i in range(m)   
            minDist = inf  #首先令最小值為無窮大 
            minIndex = -1  #令最小索引為-1
            for j in range(k)  #下面是一個嵌套循環
             #計算點數據點I和簇中心點J的距離
                distJI=distMeas(centroids[j,],dataSet[i,])
                if distJI < minDist  
                    minDist = distJI
                    minIndex = j
            if clusterAssment[i,0] != minIndex
                clusterChanged = True
            clusterAssment[i,] = minIndex,minDist**2
        print(centroids)
        #更新質心的位置
        for cent in range(k)
        #mat.A意味着將矩陣轉換為數組，即matrix-->array
            ptsInClust = dataSet[nonzero(clusterAssment[,0].A==cent)[0]]
            centroids[cent,] = mean(ptsInClust,axis=0)
    return centroids,clusterAssment

k-means與EM算法

EM是機器學習十大算法之一，是一種好理解，但是推導過程比較復雜的方法。
下面將原英文版的EM算法介紹翻譯一遍，在翻譯的過程中也加深一點自己的理解。

Jensen不等式

假設$f$是一個定義域為實數的函數。回憶一下，如果對於所有的$x \in R,f''(x) \geq 0$，則$f$是一個凸函數。如果$f$的輸入值是一個向量，則當$x$的海塞矩陣(hessian矩陣$H$)是半正定矩陣時，$f$是凸函數。如果對於所有的$x$，$f''(x)>0$恆成立，那么我們說$f$是嚴格凸函數（如果$x$是一個向量，則當$H$是嚴格半正定矩陣時，寫做$H>0$，則$f$是嚴格凸函數）。
設$f$是一個凸函數，並且$X$是一個隨機變量，則：

\[E[f(X)] \geq f(EX). \]

當且僅當$x=$常數時，$E[f(X)]=f(EX)$$.$
進一步說，如果$f$是一個嚴格凸函數，則當P{X=E(X)}=1時，滿足$E[f(X)] = f(EX).$

由於我們在寫某一隨機變量的期望時，習慣上是不寫方括號的，因此在上述式子中，$f(EX)=f(E[X])$
為了解釋上述理論，可以用圖1幫助理解。

如圖1所示，實黑線表示凸函數$f$，$X$是一個隨機變量，取值為$a$和$b$的概率均為0.5。因此，$a$和$b$的均值$E(X)$在二者之間。

與此同時，我們在y軸上可以看見$f(x)$,$f(b)$和$f(EX)$的值，並且，$E[f(X)]$$f(a)$和$f(b)$之間。從這個例子可以看出，因為$f$是一個凸函數，所以肯定滿足$E[f(X)]\ge f(EX)$。
一般來說，很多人會忘記上式的不等號方向，那么，記住這個圖，就很易能夠想起來上面的公式了。

注意，如果$f$是一個(嚴格的)凹函數($f''(x)\le 0$或者$H\le0$)，Jensen不等式仍然成立，只是方向反過來而已，即$E[f(X)]\le f(EX).$}

EM算法

假設我們有包含m個獨立樣本的訓練集$\{x^{(1)},x^{(2)},...,x^{(m)}\}$,基於該樣本和模型$p(x,z)$估計待估參數$\theta$,似然函數如下：

\[\ell (\theta)=\sum _{i=1}^mlogp(x;\theta)=\sum _{i=1}^mlog\sum _zp(x,z;\theta). \]

但是，明確找出$\theta$的最大似然估計值恐怕非常困難。因為在這里，$z^{(i)}$是隱含隨機變量，通常情況下如果已知$z^{(i)}$，那么求上述最大似然估計值會比較容易。

在這種情況下，EM算法給出了一種有效的求最大似然估計值的方法。直接求$\ell (\theta)$的最大似然估計值也許會很困難，我們的策略是在$\ell (\theta)$下面構造一個下界(E-步驟)，然后優化這個下界，使其不斷逼近求$\ell (\theta)$的最大似然估計值(M-步驟)。

對於每一個$i$，令$Q_i$是$z^{(i)}$上關於某一分布的概率，($\sum _zQ_i(z)=1,Q_i(z)\ge 0$).}得到下面的式子.如果z是連續型隨機變量，那么$Q_i$就是密度函數，$Q_i$在z上的總和就是$Q_i$在z上的積分：

\[\sum _i logp(x^{(i)};\theta)=\sum _ilog\sum _{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \]

\[=\sum _i log \sum _{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

\[\ge \sum _i \sum _{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\]

最后一步應用了Jensen不等式。特別地，在這里$f(x)=logx$是一個凹函數，因為在實數范圍內

\[f''(x)=-1/x^2<0 \]

式子

\[\sum _{z^{(i)}}Q_i(z^{(i)})[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}] \]

是式子的$[p(x^{(i)},z^{(i)};\theta)/Q_i(z^{(i)})]$關於$z^{(i)}$在$Q_i$給出的分布上的期望。根據Jensen不等式，有：

\[f(E_{z^{(i)}\sim Q_i}[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}])\ge E_{z^{(i)}\sim Q_i}[f(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})})] \]

下標$z^{(i)}\sim Q_i$表示期望是對來自於分布$Q_i$的隨機變量$z^{(i)}$所求的期望。這個條件使式子(2)到式子(3)得以成立。

現在，對於來自於$Q_i$的任意一個集合，式子(3)給出了一個$\ell (\theta)$的一個下界。$Q_i$有很多種可能的選擇，我們應該選哪個呢?如果我們當前對參數$\theta$(基於已有的條件)已經有了一些猜測，那么使這個下界與$\theta$的值接近看起來就比較自然了。換句話說，我們可以使上述的不等式在$\theta$的某個特定值的條件下變成等式。(稍后我們就能看到上述是何如使$\ell (\theta)$在EM的迭代下單調遞增的。)

為了使下界盡可能地靠近$\theta$的某個特定值，我們需要在步驟中加入上面推導的Jensen不等式以得到等式。為了實現上述步驟，我們知道只需要使期望變成常量——固定的隨機變量就行了。例如，我們需要滿足：

\[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c \]

而這個常數$c$不依賴於$z^{(i)}$。這個很簡單就可以實現，只需要使

\[Q_i(z^{(i)})\propto p(x^{(i)},z^{(i)};\theta) \]

實際上，因為我們已知$\sum _zQ_i(z^{(i)})=1$,所以有：

\[Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum _zp(x^{(i)},z^{(i)};\theta)} \]

\[=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \]

\[=p(z^{(i)}|x^{(i)};\theta) \]

因此，我們簡單假設$Q_i$是在給定$x^{(i)}$和參數$\theta$的情況下$z^{(i)}$的后驗概率.

現在，對於上述的$Q_i$，公式(3)給出了我們想取極大值的似然函數$\ell$的下界。這是E步驟。在算法里的M步驟，我們又對公式(3)取極大值以獲得新的$\theta$.重復執行上述步驟，即：

start循環直到收斂

(E-step)對於每一個i,令：

\[Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta) \]

M-step 令：

\[\theta=arg\ max_\theta \sum _i \sum _{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

我們怎么知道這個算法最終收斂呢？假設$\theta^{(t)}$和$\theta^{(t+1)}$是EM算法連續迭代兩次以后得到的參數。下面證明$\theta^{(t)} \le \theta^{(t+1)}$，即證明EM算法會使log似然概率單調遞增。證明上述的關鍵是對於$Q_i$的選擇。在以$\theta^{(t)}$為起始標志的EM算法迭代過程中，我們會選擇$Q_i^{(t)}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta^{(t)})$.我們在前面已經看到了，上述選擇會滿足Jensen不等式，像公式(3)的應用一樣，能夠得到收斂的等式，即：

\[\ell (\theta^{(t)})=\sum _i \sum _{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \]

最大化等式右邊得到參數$\theta^{(t+1)}$，因此：

\[\ell (\theta^{(t+1)})\ge \sum _i \sum _{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \]

\[ \ge \sum _i \sum _{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} =\ell (\theta^{(t)}) \]

第一個不等式來自於：

\[ \ell(\theta) \ge \sum _i \sum _{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

上式對於$Q_i$和$\theta$的任意值，尤其是$Q_i=Q_i^{(t)}$,$\theta=\theta^{(t+1)}$時成立。為了得到不等式(5)，我們依據的是$\theta^{(t+1)}$應該等於：

\[ arg\ max_\theta \ \sum _i \sum _{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

因此根據上式計算的$\theta^{(t+1)}$肯定是大於等於$\theta^{(t)}$的。最后，自然而然就推到了公式(6).

因此，EM使得似然函數單調收斂。我們們關於EM算法的描述中，我們說最終似然函數會收斂。盡管我們的結果已經展示了上述結論，但是一個令人信服的收斂性檢驗將會檢查在給定的參數容忍度的情況下，

$\ell(\theta)$會隨着EM算法的迭代而不斷增加，如果EM使得$\ell(\theta)$的上升速度慢到一定程度，就意味着收斂。

如果我們定義

\[ J(Q,\theta)=\sum _i \sum _{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} ,\]

在我們前面的推導中已經知道了$\ell (\theta)\ge J(Q,\theta)$。EM算法還可以看作使畸變函數J不斷上升的過程。E-step固定Q，最大化J；M-step固定$\theta$最大化J。而Q和$\theta$都是在迭代過程中不斷更新的。

k-means與EM

在1.8.4節對k-means算法的迭代思想進行了說明，可以看出與EM算法具有很大的相似性，在原理上幾乎相同。

k-means舉例及R實現

本文以R自帶數據集iris為例，實現聚類的k-means算法。首先，看一下iris數據集的基本結構，如表1所示。

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
1	5.10	3.50	1.40	0.20
2	4.90	3.00	1.40	0.20
3	4.70	3.20	1.30	0.20
4	4.60	3.10	1.50	0.20
5	5.00	3.60	1.40	0.20
6	5.40	3.90	1.70	0.40

該數據集共有5個變量，表1顯示的是前6行的數據，完整數據一共有150行。五個變量的中文名稱分別是：萼片長度、萼片寬度、花瓣長度、花瓣寬度和品種。以此樣本作為聚類的原始數據，屬於有標簽聚類，可以將聚類結果與真實情況進行對比，以得出聚類效果。
下面是實現k-means聚類的R語言代碼：

> iris_test<-iris[,14]
> iris_test_cl<-kmeans(x = iris_test,centers = 3,iter.max = 20)
> iris_test_cl
K-means clustering with 3 clusters of sizes 50, 62, 38

Cluster means
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1     5.006000    3.428000     1.462000    0.246000
2     5.901613    2.748387     4.393548    1.433871
3     6.850000    3.073684     5.742105    2.071053

Clustering vector
  [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 1 1 1 1 1 1 1 1 1
 [39] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [77] 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3 3 2 3 3 3 3 3 3 2
[115] 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 3 3 3 3 2 3 3 3 2 3 3 3 2 3 3 2

Within cluster sum of squares by cluster
[1] 15.15100 39.82097 23.87947
 (between_SS / total_SS =  88.4 %)

Available components

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      
> table(iris_test_cl$cluster,iris$Species)
   
    setosa versicolor virginica
  1     50          0         0
  2      0         48        14
  3      0          2        36

根據輸出結果可以看出，setosa類分類正確率為100%，versicolor類有48個分類正確，有兩個錯誤分到了virginica中，viginica有36個分類分類正確，有14個錯誤分到了versicolor中。總體來說，聚類的正確率為89.3%

k-means算法的改進

k-mode算法

k-mode算法實現了對離散型數據的快速聚類，保留了k-means算法的效率的同時，將k-means的應用范圍擴大到了離散型數據。

k-mode算法是按照k-means算法的核心內容進行修改，針對分類屬性的度量和更新質心的問題而改進。具體如下：

度量記錄之間的相關性D的計算公式是比較兩記錄之間，屬性相同為0，不同為1，並把所有的值相加。因此D越大，就說明兩個記錄之間的不相關性越強，也可以理解為距離越大，與歐氏距離代表的意義是一樣的。
更新modes，使用每個簇的每個屬性出現頻率最大的那個屬性值作為代表簇的屬性值。

k-prototype算法

k-prototype算法可以對離散型和數值型兩種混合的數據進行聚類，在k-prototype中定義了一個對數值型和離散型屬性都計算的相異性度量標准。

k-prototype是結合k-means和k-mode算法，針對混合屬性的，解決兩個核心問題如下：

度量具有混合屬性的方法是，數值型屬性采用k-means方法得到P1，分類屬性采用k-mode方法得到P2，那么D=P1+a*P2，a是權重，如果覺得分類屬性重要，則增加權重的值。當a=0時，只有數值型屬性。
更新一個簇的中心的方法是結合k-means和k-mode的方法。

k-中心點算法

k-中心點算法是針對k-means算法對於孤立點敏感所提出的改進方法。為了解決上述問題，不采用簇中的平均值作為參照點，可以選用簇中位置最中間的對象，即中心點作為參照點。這樣划分方法仍然基於最小化所有對象與其參照點的相異度之和的原則來執行的。
對於算法的實現來說，就是在更新質心的時候，不是計算所有點的平均值，而是中位數來代表其中心。然后進行迭代，直到質心不再變化為止。

Enhanced k-means

下面討論一下k-means聚類的收斂性。

任意生成k個初始類中心$\mu_1,\mu_2,...,\mu_k \in R$.

重復如下步驟直到收斂：

對於每一個i，令

\[c^{(i)}=arg \ {min_j ||x^{(i)}-\mu_j||^2} \]

對於每一個j，令

\[\mu_j=\frac{\sum _{i=1}^m1\{c^{(i)}=j\}x^{(i)}} {\sum _{i=1}^m1\{c^{(i)}=j\}} \]

那么，k-means算法是否能夠保證一定收斂呢？答案是肯定的。
定義畸變函數

\[J(c,\mu)=\sum _{i=1}^n||x^{(i)}-\mu _c^{(i)}||^2 \]

k-means算法的目的就是使J降至最小。在內循環中，固定$\mu$,可以通過調整$c$來使J減小；同樣的，固定$c$，調整每個類的質心$\mu$也可以使J減小。當J減到最小時，$\mu$和$c$也同時收斂。

但是，需要注意的是，J函數是一個非凸函數，所以J不能保證收斂到全局最優，而可以保證收斂到局部最優。通常，k-means達到局部最優的結果差不多是全局最優。但是如果擔心陷入了很嚴重的局部最優，可以多運行幾次k-means算法（(給出不同的初始類中心)，選出使J達到最小值的那個模型的$\mu$和$c$。

為了克服k-means上述缺點，可以在k-mean的結果上對其進行改進(enhanced k-means)。該算法步驟如下：

根據k-means算法得出k個類和相應的類中心。

計算每一個樣本點與所有類中心的歐氏距離。

假設$x_i$在第$r$個類中，$n_r$表示第$r$個類包含的樣本點數目，$d_{ir}^2$表示$x_i$與第$r$個類的中心點的歐氏距離。如果存在類$s$，使得

\[\frac{n_r}{n_r-1}d_{ir}^2>\frac{n_s}{n_s+1}d_{is}^2, \]

則將$x_i$移到類$s$中。
如果有多個類滿足上述不等式，則將$x_i$移動到使得
$\frac{n_s}{n_s+1}d_{ir}^2$最小的那個類中。
重復步驟2$\sim$4，直到沒有變化為止。

二分K-均值算法

為了克服k-means算法收斂於局部最小值的問題，有人提出了使用另一個稱為二分-K均值(bisecting K-means)的算法。該算法首先將所有的點作為一個簇，然后將該簇一分為二。之后選擇其中一個簇繼續划分，選擇哪一個簇進行划分取決於對其划分是否可以最大程度降低SSE的值。上述基於SSE的划分過程不斷重復，直到得到用戶指定的簇數目為止。

二分K-均值算法的步驟如下：

將所有的點看成一個簇
當簇數目小於k時，對於每一個簇

計算總誤差
在給定的簇上面進行K-means聚類(k=2)
計算將該簇一分為二后的總誤差

選擇使得誤差最小的那個簇進行划分操作。

另一種方法是選擇SSE最大的那個簇進行划分。該算法的R代碼如下。

#計算誤差平方和SSE，group接受一個數值型矩陣，centroid為其均值向量


SSE<-function(group,centroid)
{
  rows= nrow(group)
  group = as.matrix(group)
  centroid = matrix(rep(centroid,rows),nrow = rows,byrow = TRUE)
  group = group - centroid
  group = group * group
  sums = apply(X = group,MARGIN = 1,FUN = sum)
  SSE = sum(sums)
  return(SSE)
}
#定義二分K-means聚類函數，dataSet為需要聚類的數據集，k為指定聚類個數
biKmeans<-function(dataSet,k)  
{
  if(!is.matrix(dataSet))
    dataSet = as.matrix(dataSet)
  if(k==1)
    return(dataSet)
  else
  {
    #定義一個空矩陣來存放不需要再次處理的已經聚好的類
    clusteredDataSet = matrix(nrow = 0,ncol = ncol(dataSet))
    #定義一個空向量用來存放上述類的標簽
    clustered = vector(length = 0)
    #計數器clusters
    clusters = 1
    #i用來改變每次聚類的標簽
    i = 0
    while(clusters < k)
    {
      #使用k-means用來二分聚類
      cl = kmeans(dataSet,2)
      #類標簽
      cluster = cl$cluster+i
      #計算聚好的兩個類的與類中心的離差平方和
      SSE1 = SSE(group = dataSet[cluster==(1+i),],centroid = cl$centers[1,])
      SSE2 = SSE(group = dataSet[cluster==(2+i),],centroid = cl$centers[2,])
      #將誤差平方和大的拿出來，下次繼續二分
      if(SSE1 > SSE2)
        index = (cluster == (1+i))
      else
        index = (cluster ==(2+i))
      #將不需要再次聚類的類放在clusterDataSet例
      clusteredDataSet = rbind(clusteredDataSet,dataSet[!index,])
      #將上述類的標簽存放在clustered里
      clustered = c(clustered,cluster[!index])
      #更新dataSet，下次繼續二分
      dataSet = dataSet[index,]
      #計數器+1
      clusters = clusters +1
      #i+2，因為每次二分的結果是1和2，在結果的基礎上加上2，則類標簽可以明確區分
      i = i + 2
    }#end while
    #得出最終聚類結果
    clusteredDataSet = rbind(clusteredDataSet,dataSet)
    #類標簽
    clustered = c(clustered,cluster[index])
    #轉換為數據框
    clusteredDataSet = as.data.frame(clusteredDataSet)
    #添加類標簽列
    clusteredDataSet$cluster = clustered
    return(clusteredDataSet)
  }
}
> iris_test_bik<-biKmeans(iris_test,3)
> dim(iris_test_bik)
[1] 150   5
> table(iris_test_bik$cluster,iris$Species)
   
    setosa versicolor virginica
  1     50          3         0
  3      0          9        50
  4      0         38         0

由上述最終的聚類結果顯示，聚類正確率約為92%，比單純使用k-means算法正確率高。

Enhanced k-means算法更加穩健的原因

k-means算法具有明顯的兩個缺點，即容易陷入局部最優和出現“超級類”。局部最優比較好理解，“超級類”是指聚類結果中有一個包含許多樣本點的大類，其他的類包含的樣本點都較少。enhanced k-means同時考慮了樣本點與聚類中心的距離和類中樣本點的個數，保留了k-means聚類算法的優點，同時又克服了出現"超級類"的問題。

Enhanced k-means舉例及R實現

k-中心點算法可以使用R中cluster包中的pam函數(pam的全稱為Partitioning Around Medoids，即圍繞中心點分割)，同樣以iris數據集為例，說明pam函數的使用。

> library(cluster)
> pam.cl<-pam(iris_test,k = 3)
> table(pam.cl$clustering,iris$Species)
   
    setosa versicolor virginica
  1     50          0         0
  2      0         48        14
  3      0          2        36
> par(mfrow=c(1,2))
> plot(pam.cl)

圖2將原始數據進行了降維，提取了兩個主成分，並且這兩個主成分能夠解釋原始數據方差的95.81%，說明降維效果很好。將這兩個成分繪制成散點圖，就得到了圖2.結合table函數的輸出結果和圖2可知，setosa類與其他兩個類的界限比較明顯，而其他兩個類之間具有交叉。

kmeans_4.PNG
kmeans_4.PNG

聚類指標評價

多種聚類結果的比較

聚類性能的度量大致有兩類，一類是將聚類結果與某個“參考模型”進行比較，稱為“外部指標”；另一類是直接考察聚類結果而不利用任何參考模型，稱為“內部指標”。

外部指標

對數據集$D=\{x_1,x_2,...,x_m\}$,假定通過聚類給出的簇划分為

\[C=\{C_1,C_2,...,C_k\}, \]

參考模型給出的簇划分為

\[C^+=\{C_1^+,C_2^+,...,C_s^+\} \]

相應地，令 $\lambda$ 與 $\lambda^+$ 分別表示與 $C$ 和 $C^+$ 對應的簇標記向量，我們將樣本兩兩配對考慮，定義

\[a=|SS|,SS=\{(x_i,x_j)|\lambda _i=\lambda _j,\lambda _i^+=\lambda _j^+,i<j \} \]

\[b=|SD|,SD=\{(x_i,x_j)|\lambda _i =\lambda _j,\lambda _i^+ \ne \lambda _j^+,i<j \} \]

\[c=|DS|,DS=\{(x_i,x_j)|\lambda _i \ne \lambda _j,\lambda _i^+=\lambda _j^+,i<j \} \]

\[d=|DD|,DD=\{(x_i,x_j)|\lambda _i \ne \lambda _j,\lambda _i^+ \ne \lambda _j^+,i<j \} \]

其中集合$SS$包含了在$C$中隸屬於相同簇且在$C^+$中也隸屬於相同簇的樣本對，集合$SD$包含了在$C$中隸屬於相同簇但在$C^+$中隸屬於不同簇的樣本對，……由於每個樣本對 $(x_i,x_j)(i<j)$ 僅能出現在一個集合中，因此有 $a+b+c+d=m(m-1)/2$成立。

基於以上給出的信息，可以導出下面常用的聚類性能度量外部指標。

Jaccard系數(Jaccard Coefficient,簡稱JC)

\[JC=\frac{a}{a+b+c}. \]

FM指數(Fowlkws and Mallows Index,簡稱FMI)

\[FMI=\sqrt{\frac{a}{a+b}\times \frac{a}{a+c}}. \]

Rand指數(Rand Index,簡稱RI)

\[RI=\frac{a+d}{m(m-1)/2}=\frac{a+d}{C_m^2}. \]

顯然，上述性能度量的結果均在[0,1]之間，值越大越好。}

內部指標

考慮聚類結果的划分$C=\{C_1,C_2,...,C_k\}$,定義：

\[avg(C)=\frac{2}{|C|(|C|-1)} \sum _{1\le i<j \le |C|}dist(x_i,x_j), \]

\[diam(C)=max_{1\le i<j \le |C|}dist(x_i,x_j), \]

\[d_{min}(C_i,C_j)=min_{x_i\in C_i,x_i\in C_j}dist(x_i,x_j), \]

\[d_{cen}(C_i,C_j)=dist(\mu _i,\mu _j), \]

其中，$dist(·,·)$ 用於計算兩個樣本之間的距離，$\mu$代表$C$的中心點$\mu =\frac{1}{|C|}\sum _{1\le i \le |C|}x_{i}.$顯然，$avg(C)$對應於$C$內樣本間的平均距離，$diam(C)$對應於簇$C$內樣本間的最遠距離，$d_{min}(C_i,C_j)$對應於簇$C_i$與簇$C_j$最近樣本間的距離，$d_{cen}(C_i,C_j)$對應於簇$C_i$與簇$C_j$中心點的距離。

基於上式可以導出下面這些常用的聚類性能度量內部指標。

DB指數(Davies-Bouldin Index,簡稱DBI)

\[DBI=\frac{1}{k}\sum _{i=1}^k max_{j\ne i}(\frac{avg(C_i)+avg(C_j)}{d_{cen}(\mu _i,\mu _j)}). \]

Dunn指數(Dunn Index，簡稱DI)

\[DI=min_{1\le i \le k}\{min_{j\ne i}(\frac{d_{min}(C_i,C_j)}{max_{1\le l \le k}\ diam(C_l)})\} \]

顯然，DBI的值越小越好，而DI則相反，值越大越好。

共表型相關系數

共表型相關系數(Cophenetic Correlation Coefficient，簡稱CPCC)

內容參考：http//people.revoledu.com/kardi/tutorial/Clustering/index.html.

另外，共表相關系數也屬於內部指標，這里單獨拿出來是因為其計算比較復雜，且僅僅用在層次聚類方法上。共表型相關系數是用於衡量層次聚類聚類性能的指標，為了清楚說明共表相關系數的計算，下面先從層次聚類說起。
層次聚類的標准輸出結果是一個樹狀圖，樹狀圖的橫坐標顯示的是樣本點，縱坐標是樹的高度。樹狀圖可以看作是層次聚類的可視化結果。

使用樹狀圖，我們可以很容易地找出決定聚類個數的分割點。比如說，在圖4中，如果將分割高度設置為2，則將6個樣本聚成2類(一類包含了A、B樣本，一類包含了C、E、D、F)；如果將分割高度設置為1.2，則將6個樣本聚成3類(一類包含A、B，一類包含C，一類包含E、D、F)。

層次聚類算法一共有兩種實現方法。第一種方法是$agglomerative \ approach$，先從底部將所有的樣本點單獨看作一類，然后將兩個最近的樣本點聚在一起，形成一個新的類。接着計算所有類之間的聚類，將兩個最近的類聚在一起。上述步驟重復下去，直到把所有的點聚為一類。該方法也稱為自底向上的層次聚類方法。第二種方法是$divisive\ approach$，即先從頂部開始，將所有的樣本點看作一個大類。然后將該大類分成兩個小類，直到所有的樣本點自成一類結束。上述層次聚類的一個可行的方法是先將上述所有樣本生成一個最小生成樹(例如使用Kruskal算法)，然后通過最大距離將該樹分割。這種方法也被稱為自頂向下的層次聚類方法。

這里使用第一種方法演示層次聚類方法，步驟如下：

計算各個變量之間的聚類，生成距離矩陣

將每個樣本點視為一個類

迭代，直到類的個數為1

將最近的兩個類合並

更新距離矩陣

為了以圖示方法演示層次聚類算法，本文引入了一個小例子。假設我們有6個對象，每個對象有兩個特征，具體數據如表2所示。我們可以繪制X1與X2的散點圖觀察對象之間的關系，如圖5所示。可以看到，A與B的距離較近，D、F、E的距離較近，C有點離群。可以通過計算樣本間的距離矩陣來反應樣本之間的關系。一般來說，計算距離會采用上面所說的歐氏距離。注意，距離矩陣的輸出結果肯定是一個對稱矩陣。這里的輸出結果如表3所示。

-	X1	X2
A	1.00	1.00
B	1.50	1.50
C	5.00	5.00
D	3.00	4.00
E	4.00	4.00
F	3.00	3.50

-	A	B	C	D	E	F
A	0.00	0.71	5.66	3.61	4.24	3.20
B	0.71	0.00	4.95	2.92	3.54	2.50
C	5.66	4.95	0.00	2.24	1.41	2.50
D	3.61	2.92	2.24	0.00	1.00	0.50
E	4.24	3.54	1.41	1.00	0.00	1.12
F	3.20	2.50	2.50	0.50	1.12	0.00

因為距離矩陣是一個對稱矩陣，所以一般只采用上三角或者下三角數據，本文采用的是下三角矩陣。對角線上的元素全為0，因為計算的是樣本點與自己的距離。總的來說，如果有m個樣本點，那么下三角矩陣包含有$\frac{1}{2}m(m-1)$的數據。在我們的例子里，一共有$1/2\times6\times(6-1)=15$個元素。可以看出，在距離矩陣中最小的元素是0.5，即D與F之間的距離。如果樣本點包含的還有分類數據，則可以采用其他方法計算樣本點之間的距離。

如何將對象聚在一個類中是層次聚類的關鍵。給出一個距離矩陣，對象之間的距離可以通過計算類之間的距離的一些標准來確定。最基礎和最常用的幾個方法如下。

Single Linkageminimun distance criteron

\[d_{A\rightarrow B}=min_{\forall i \in A,\forall j \in B}(d_{ij}) \]

Complete Linkagemaximun distance criteron

\[d_{A\rightarrow B}=max_{\forall i \in A,\forall j \in B}(d_{ij}) \]

Average Groupaverage distance criteron

\[d_{A\rightarrow B}=average_{\forall i \in A,\forall j \in B}(d_{ij}) \]

Centroid distance criteron

\[d_{A\rightarrow B}=||c_A-c_B||=\frac{1}{n_in_j}\sum_ {\forall i \in A,\forall j \in B}(d_{ij}) \]

已存在的包含有$n_k$個樣本的類k與新聚成的類$(r,s)$之間的距離公式如下：

\[d_{k \rightarrow (r,s)} =\alpha_r d_{k\rightarrow r}+\alpha_s d_{k\rightarrow s}+\beta d_{r\rightarrow s}+ \gamma |d_{k\rightarrow r}-d_{k\rightarrow s}| \]

具體的參數如表4所示。

clustering method	$\alpha_r$	$\alpha_s$	$\beta$	$\gamma$
Single Link	1/2	1/2	0	-1/2
Complete Link	1/2	1/2	0	1/2
Unweighted pair group method average(UPGMA)	$\frac{n_r}{n_r+n_s}$	$\frac{n_s}{n_r+n_s}$	0	0
weighted pair group method average(WPGMA)	1/2	1/2	0	0
unweighted pair group method centroid(UPGMC)	$\frac{n_r}{n_r+n_s}$	$\frac{n_s}{n_r+n_s}$	$\frac{-n_sn_s}{(n_r+n_s)^2}$	0
weighted pair group method centroid(WPGMC)	1/2	1/2	-1/4	0
Ward's method	$\frac{n_r+n_k}{n_r+n_s+n_k}$	$\frac{n_s+n_k}{n_r+n_s+n_k}$	$\frac{-n_k}{n_r+n_s+n_k}$	0

最小距離層次聚類法也被稱為單鏈聚類或者最近鄰聚類。兩個類之間的距離定義為兩個類中對象之間距離的最小值。比如對於該例來說，最小的是D和E的距離0.50，選擇先把D和E聚為一類(D,F)。然后其他單個樣本之間的距離不變，計算(D,F)和其他樣本之間的距離，采用最小距離法，得到的結果為：$d((D,F),A)=d(F,A)=3.20,d((D,F),B)=d(B,D)=2.92,d((D,F),C)$
$=d(D,C)=2.24,d((D,F),E)=d(D,E)=1.00$,然后可以與原來的樣本間距離放在一起，發現0.71——A與B之間的距離最小，因此將A與B聚為一類，得到(A,B).依次類推，最終將所有的樣本點聚成一個大類。總結一下計算過程如下：

在一開始，我們有6個類，即A，B，C，D，E，F，

將距離最小為0.5的D和F聚為一類(D,F)

將距離最小為0.71的A和B聚為一類(A,B)

將距離最小為1.00的(D,F)和E聚為一類(D,F,E)

將距離最小為1.41的(D,F,E)和C聚為一類(D,F,E,C)

將距離最小為2.50的(D,F,E,C)和(A,B)聚為一類(D,F,E,C,A,B)

最后一個類包含了所有的樣本點，聚類過程結束

根據上述過程，我們可以輕松的把層次聚類的樹狀圖畫出來。

那么，這個聚類的結果性能如何呢？有一個被稱為共表相關系數的指標可以用來展示聚類結果的擬合優度。

為了計算層次聚類的共表相關系數，我們需要以下兩個信息：

距離矩陣
共表矩陣

我們在上文已經計算出來了樣本之間的距離矩陣，現在只需要計算共表矩陣就行了。我們需要使用在合並類的過程中使用的最小距離來填充距離矩陣的下三角。上面我們已經將合並時所使用的距離進行了總結，利用上面的總結信息，可以得到如表所示的共表矩陣。

-	A	B	C	D	E
A
B	0.71
C	2.50	2.50
D	2.50	2.50	1.41
E	2.50	2.50	1.41	1.00
F	2.50	2.50	1.41	0.50	1.00

將距離矩陣和共表矩陣一一對應生成兩個向量，計算這兩個向量之間的相關系數，就會得到共表相關系數。在這里，最終CPCC=86.339%,可以看到，這個聚類的效果比較好。

這里以R中自帶數據集iris作為樣本，使用cluster包中的hclust函數進行層次聚類，分析聚類結果。

> library(cluster)
> iris_test <- iris[,14]
> d <-dist(iris_test)
> hc1 <- hclust(d,'complete')  #最大距離法
> d1 <- cophenetic(hc1)
> cor(d,d1)
[1] 0.7269857
> 
> hc2 <- hclust(d,'single')  #最小距離法
> d2 <- cophenetic(hc2)
> cor(d,d2)
[1] 0.8638787

從上述R中的運行結果來看，single方法(最小距離法)的結果要優於complete法(最大距離法)。

聚類個數的選擇

輪廓系數(Silhouette Coefficient,簡稱SC)

對於數據集$x=\{x_1,x_2,...,x_n\}$,使用某種聚類方法得到類的集合$C=\{C_1,C_2,...,C_k\}$,每個類的樣本量為$num=\{n_1,n_2,...,n_k\}$}

\[a_i=\frac{1}{n_k-1}\sum _{x_i,x_j\in C_k}dist_{i\ne j}(x_i,x_j) \]

\[b_i=min(\frac{1}{n_l}\sum _{x_i\in C_k,x_j\in C_l,l\ne k}dist(x_i,x_j)) \]

\[S_i=\frac{b_i-a_i}{max(b_i,a_i)} \]

\[SC=\frac{1}{n}\sum _{i=1}^n S_i \]

SC即為輪廓系數。該值處於-1-1之間，值越大，表示聚類效果越好。

Gap統計量(Gap Statistic,簡稱GS)}

下面簡單說一下Gap統計量原理及計算方法。}
Tibshirani等人在2001年提出了使用Gap統計量來確定聚類過程中最佳聚類個數k，這個統計量除了可以使用在最常用的k-means聚類方法上，對於任何聚類方法都適用。

假設有數據集$\{x_{ij}\},i=1,2,3...,n,j=1,2,...,p$，該數據有n個觀測點和p個特征。令$d_{ii'}$表示觀測點$i$和觀測點$i'$之間的距離，最常用的是歐氏距離。

假設我們已經將數據聚成$k$個類，分別是$C_1,C_2,...,C_k$，其中$C_r$表示第$r$個類的下標，令$n_r=|C_r|.$即$n_r$表示類$C_r$中的觀測點的個數。令

\[D_r=\sum _{i,i'\in C_r}d_{ii'}=\sum _{i,i'\in C_r}|x_i-x_i'|^2=2n_r\sum _{i\in C_r}|xi-\bar x|^2 \]

\[W_k=\sum _{r=1}^k\frac{1}{2n_r}D_r \]

也可以使用$W_k$來確定聚類個數$k$，尋找使得$W_k$突然變小的$k$.但是使用$W_k$來確定聚類個數有以下兩個問題：}

沒有相關的類進行比較。
$W_k-W_{k-1}$的值沒有進行標准化，以至於不同的模型無法進行比較。

Gap統計量將$W_k$進行了標准化，公式為

\[Gap_n(k)=E_n^+\{log(W_k)\}-log(W_k) , \]

尋找使得$Gap(k)$達到最大的那個$k$，即為最佳的聚類個數。
其中，$E_n^+$表示根據某一相關分布構造的樣本量為n的數據，其
$log(W_k)$的期望值。
通常，構造上述的樣本具有以下兩種方法，即：

計算出每個特征（或變量）的取值范圍，然后在該范圍內隨機生成n個服從均勻分布的數字，作為構造的樣本觀測值。
如果$X$是一個$n\times p$的數值矩陣，假設每一列的均值都是0，然后對$X$進行奇異值分解，

\[X=UDV^T. \]

將$X$進行轉換，得到$X'$

\[X'=XV. \]

然后利用方法1根據$X'$得到每一列的數值$Z'$。
最后，將$Z'$進行轉換，

\[Z=Z'V^T \]

最終得到我們需要的數據$Z$。

計算Gap統計量的步驟如下。

將數據聚類，類的個數分別取$k=1,2,...,K$，得出$W_k,k=1,2,...,K.$
根據均勻分布，使用上面的方法生成B個數據集。計算出每個數據集的
$W_{kb}^+,b=1,2,...,B,k=1,2,...,K$，計算出估計的Gap統計量的值：

\[Gap(k)=(1/B)\sum _b log(W_{kb}^+)-log(W_k) \]

令

\[\bar l = (1/B)\sum _b log(W_{kb}^+)$$，計算標准誤差 $$sd_k=[(1/B)\sum _b \{log(W_{kb}^+)-\bar l\}^2]^{1/2}\]

定義$s_k=sd_k\sqrt {1+1/B}$，最后，計算出最佳的聚類個數$\hat k$，

\[\hat k=smallest\ k\ such\ that\ Gap(k) \le Gap(k+1)-s_{k+1} \]

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 kmeans 聚類 Kmeans聚類與層次聚類 Python之聚類（KMeans,KMeans++） SparkMLlib聚類學習之KMeans聚類 tensorflow kmeans 聚類 kmeans聚類理論篇聚類學習筆記——kmeans 淺談Kmeans聚類 Kmeans聚類_R代碼 Kmeans 聚類及其python實現