西瓜書習題試答-第09章-聚類

試答系列：「西瓜書」-周志華《機器學習》習題試答
系列目錄
[第01章：緒論]
[第02章：模型評估與選擇]
[第03章：線性模型]
[第04章：決策樹]
[第05章：神經網路]
[第06章：支援向量機]
[第07章：貝葉斯分類器]
[第08章：集成學習]
第09章：聚類
第10章：降維與度量學習
第11章：特徵選擇與稀疏學習
第12章：計算學習理論(暫缺)
第13章：半監督學習
第14章：概率圖模型
(後續章節更新中…)

目錄

9.1 試證明：p≥1時，閔可夫斯基距離滿足距離度量的四條基本性質；0≤p<1時，閔可夫斯基距離不滿足直遞性，但滿足非負性、同一性、對稱性；p趨向無窮大時，閔可夫斯基距離等於對應分量的最大絕對距離，即

\[\lim_{p\rightarrow \infty}(\sum_{u=1}^n|x_{iu}-x_{ju}|^p)^{\frac{1}{p}}=\max_u |x_{iu}-x_{ju}|
\]

證明：

  • 非負性、同一性、對稱性
    對於所有的p≥0,由於距離公式中的絕對值形式，決定了最終閔可夫斯基距離的\(dist_p(x_i,x_j)\)的非負性、同一性、對稱性：\(|x_{iu}-x_{ju}|≥0\)，當且僅當\(x_{iu}=x_{ju}\)時等於0，而且\(|x_{iu}-x_{ju}|=|x_{ju}-x_{iu}|\)。

  • 直遞性
    基於閔可夫斯基（Minkowski）不等式，可以證得p≥1時滿足直遞性，0≤p≤1時不滿足。

    定理：閔可夫斯基（Minkowski）不等式
    對任意p≥1以及\(x,y\in R^n\),有

    \[|x+y|_p≤|x|_p+|y|_p
    \]

    其中,\(|x|_p=(\sum|x_i|_p)^{1/p}=distm_k(x,0)\)。
    如果，\(x_1,…,x_n,y_1,…,y_n>0，p<1\)，則”≤”可以變為”≥”。

    令\(x=x_i-x_k，y=x_k-x_j\)，即得：p≥1時的直遞性滿足，p<1時，直遞性不滿足。
    p≥1時的閔可夫斯基（Minkowski）不等式的證明過程如下：

    首先需要證明一個引理和Holder不等式，然後再證明閔可夫斯基不等式。

    引理 \(a^λb^{1-λ}≤λa+(1-λ)b\)，其中a,b≥0, 0≤λ≤1。
    證明思路：兩邊取對數，由於對數ln函數是凸函數，利用Jensen不等式可證得。

    定理(Holder不等式)
    對任意的\(1\leq p,q\leq \infty,\frac{1}{p}+\frac{1}{q}=1\)，以及\(x,y\in R^n\)，有

    \[\sum_{i=1}^n|x_i y_i|\leq|x|_p |x|_q
    \]

    當p=q=2時，Holder不等式退化為柯西-施瓦茨不等式(文字可表述為：兩個向量的內積不大於兩個模長之積)。
    證明：由上面的引理有

    \[\frac{|x_i||y_i|}{|x|_p|y|_p}=\left(\frac{|x_i|^p}{|x|_p^p}\right)^{1/p}\cdot\left(\frac{|y_i|^q}{|y|_q^q}\right)^{1/q}\leq\frac{1}{p}\frac{|x_i|^p}{|x|_p^p}+\frac{1}{q}\frac{|y_i|^q}{|y|_q^q}
    \]

    不等式兩邊對i求和便有：

    \[\frac{1}{|x|_p|y|_p}\sum_{i=1}^n |x_i y_i|\leq \frac{1}{p}+\frac{1}{q}=1
    \]

    定理得證。

    定理(Minkowski不等式)
    內容見前文
    證明：只需考慮1<p<∞的情況，p=1或者∞的情形易證。當1<p<∞時有

    \[\begin{aligned}
    |x+y|_p^p&=\sum_i|x_i+y_i|^p\\
    &=\sum_i |x_i+y_i|^{p-1}|x_i+y_i|\\
    &\leq\sum_i |x_i+y_i|^{p-1}(|x_i|+|y_i|)
    \end{aligned}\]

    由Holder不等式

    \[\begin{aligned}
    \sum_i |x_i+y_i|^{p-1}|x_i|&\leq|(x+y)^{p-1}|_q|x|_p\\
    &=\left(\sum_i |x_i+y_i|^{(p-1)q}\right)^{1/q}|x|_p\\
    &=\left(\sum_i |x_i+y_i|^p\right)^{1/q}|x|_p\\
    &=|x+y|_p^{p/q}|x|_p
    \end{aligned}\]

    第三行的利用了關係(p-1)q=p，同理有：

    \[\sum_i |x_i+y_i|^{p-1}|y_i|\leq|x+y|_p^{p/q}|y|_p
    \]

    結合以上三個不等式有：

    \[|x+y|_p^p\leq|x+y|_p^{p/q}(|x|_p+|y|_p)
    \]

    不等式兩邊同除\(|x+y|_p^{p/q}\)即得閔可夫斯基不等式。

    以上證明過程參考了這篇博文，對於p<1時的情況，暫未找到證明方法。

  • \(p\rightarrow\infty\)時的閔可夫斯基距離
    設第u個特徵差(坐標差)：\(|x_{iu}-x_{ju}|=Δu\)，最大特徵差為\(\max_uΔ_u=ΔM\)，則有

    \[\begin{aligned}
    \lim_{p\rightarrow\infty}dist_{mk}(x_i,x_j)&=\lim_{p\rightarrow\infty}[\sum_u (\Delta_u)^p]^{1/p}\\
    &=\Delta M\cdot \lim_{p\rightarrow\infty}[\sum_u (\frac{\Delta_u}{\Delta M})^p]^{1/p}\\
    &=\Delta M\cdot [\sum_u Ⅱ(\Delta_u = \Delta M)]^0\\
    &=\Delta M\\
    &=\max_u |x_{iu}-x_{ju}|
    \end{aligned}\]

    其中\(∑_uⅡ(\Delta_u=\Delta_M)\)代表特徵差等於最大特徵差的個數，其結果至少有1個。

9.2 同一樣本空間中的集合X與Z之間的距離可通過「豪斯多夫距離」(Hausdorff distance)計算：

\[dist_H (X,Z)=max(dist_h (X,Z),dist_h(Z,X))\tag{9.44}
\]

其中

\[dist_h(X,Z)=\max_{x\in X}\min_{z\in Z}||x-z||_2\tag{9.45}
\]

試證明：豪斯多夫距離滿足距離度量的四條基本性質。
證明：首先，我們來直觀的理解一下豪斯多夫距離。
在這裡插入圖片描述
下面表示兩個點的歐式距離時，為了簡便，右下角的角標2略去不寫，比如將\(|x-y|_2\)表示為\(|x-y|\)。
如上圖所示，關於\(dist_h(X,Z)\)，我們把X想像成一個發光體(就像太陽)，把Z想像成一個受光體(就像地球)。太陽上的某個發光點x發射的光線只要到達地球上任意一個點，我們即認為完成了光線傳播，在Z上最先接收到光線的點\(z^*\)稱之為受光點。於是，\(z^*=argmin_z|x-z|\)，將距離\(|x-z^*|\)稱作為x點到Z的傳播距離，那麼\(dist_h(X,Z)=\max_x\min_z |x-z|\)表示X到Z的最遠光線傳播距離。
而\(dist_h(Z,X)\)考察的是Z到X的光線傳播距離，此時，Z是發光體，X是受光體。\(dist_h(X,Z)\)和\(dist_h(Z,X)\)未必相等，兩者中較大者即為豪斯多夫距離\(dist_H(X,Z)\)，代表著X和Z集合之間的最遠光線傳播距離。

豪斯多夫距離的非負性和對稱性顯而易見，下面只證明同一性和直遞性。
同一性：如果兩個集合的豪斯多夫距離為零，則兩個集合相等。
\(dist_H(X,Z)=0\)→ \(dist_h(X,Z)=0\) →任意\(x\in X，min_z|x-z|=0\) →任意\(x\in X\)，都存在\(z=x\) →\(X\subseteq Z\)，同理有\(Z\subseteq X\)，因此X=Z。
直遞性：(證明過程參考了這篇博文，在理解原文基礎上進行了修改)
在這裡插入圖片描述
觀察上圖，將X，Y，Z故意畫成不同的形狀，為了方便繪圖和討論，這裡假設集合內元素連續分布，其結果同樣適用於離散元素分布的情況。\(dist_h(X,Z)\)對應的距離為\(|x_1-z_1|\)（\(x_1\)和\(z_1\)分別為發光點和受光點，後面類似），\(dist_h(X,Y)\)對應的距離為\(|x_2-y_1|\)，\(dist_h(Y,Z)\)對應的距離為\(|y_2-z_2|\)，可以證明\(dist_h(X,Z)≤dist_h(X,Y)+ dist_h(Y,Z)\)：
設\(x_1\)點到Y受光點為\(y_3\)，亦即\(x_1\)到Y的最短距離為\(|x_1-y_3|\)；類似的，\(y_3\)到Z的受光點為\(z_3\)，那麼有：

\[\begin{aligned}
dist_h(X,Z)&=|x_1-z_1|\\
&\leq|x_1-z_3|\\
&\leq|x_1-y_3|+|y_3-z_3|\\
&\leq|x_2-y_1|+|y_2-z_2|\\
&=dist_h(X,Y)+dist_h(Y,Z)
\end{aligned}\]

上式中，第2行的不等式是由於\(|x_1-z_1|\)是\(x_1\)到Z的最近距離，第3行是應用了距離的三角不等式性質，第4行是根據定義，\(dist_h(A,B)\)是A到B的最遠光線傳播距離。
同理，可以證明 \(dist_h(Z,X)≤dist_h(Z,Y)+ dist_h(Y,X)\)。
不失一般性，假設\(dist_h(Z,X)≥dist_h(X,Z)\)，則有：
\(dist_H(Z,X)=dist_h(Z,X)≤dist_h(Z,Y)+ dist_h(Y,X)≤dist_H(Z,Y)+ dist_H(Y,X)\)
豪斯多夫距離的直遞性得證。

9.3 試析k均值演算法能否找到最小化式(9.24)的最優解。

答：k均值演算法的運行結果依賴於初始選擇的聚類中心，找到的結果是局部最優解，未必是全局最優解。

9.4 試編程實現k均值演算法，設置三組不同的k值、三組不同初始中心點，在西瓜數據集4.0上進行實驗比較，並討論什麼樣的初始中心有利於取得好結果。

答：編程程式碼附後。
對於k值，考慮三種情況：2，3，4。
對於初始中心點，從直覺上判斷，初始中心點應該盡量分散一些較好。結合西瓜數據集4.0的具體數據分布，通過手動選取方法來確定三組不同的初始點：
相互靠攏的幾個點：5,7,17,18
相互分散，靠近數據區域外周的幾個點：29,15,10,25
相互分散，位於數據區域中部的幾個點：27,17,12,21。
在這裡插入圖片描述
運行結果如下(小黑點是數據點，藍色圓圈是初始中心，紅色十字點是最終得到的聚類中心，紅色虛線是聚類邊界)：
在這裡插入圖片描述

上圖中第一列為初始聚類中心比較靠近，第二列的初始中心比較分散而靠近外周，第三列的初始中心分散而居中。
討論：在k=2和3時，不同的聚類中心所得結果相差不明顯，所得結果相近；k=4時，初始中心選擇分散時，最終的平方誤差(按9.24式計算)更小一些。因此，就當前所嘗試的有限情況而言，貌似可以得到結論：當k較小時，不同的初始中心選取差別不大，當k較大時，選取，選擇較分散的初始中心更有利於取得較好結果.

9.5 基於DBSCAN的概念定義，若\(x\)為核心對象，由\(x\)密度可達的所有樣本構成的集合為\(X\)。試證明：\(X\)滿足連接性(3.39)與最大性(9.40)。

證明：這個結論顯然的，設\(x_i, x_j \in X\)，由於\(X\)是由所有\(x\)密度可達的樣本構成，因此\(x_i, x_j\)都可由\(x\)密度可達，於是\(x_i, x_j\)密度可連，連接性成立。
若\(x_j\)由\(x_i\)密度可達，而\(x_j\)又由\(x\)密度可達，於是\(x_j\)由\(x\)密度可達，因此，\(x_j\in X\)，最大性成立。

9.6 試析AGNES演算法使用最小距離和最大距離的區別。

答：
參考1 (from 百度文庫)：

  • 最小距離又叫單鏈接，容易造成一種叫做Chaining的效果，兩個cluster明明從「大局」上離得比較遠，但是由於其中個別的點距離比較近就被合併了，並且這樣合併之後的- Chaining效應會進一步擴大，最後會得到比較鬆散的cluster。
  • 最大距離又叫全鏈接，則是單鏈接的反面極端，其效果剛好相反，限制非常大，兩個cluster即使已經很接近了，但是只要有不配合的點存在，就頑固到底，老死不相合併，也不是太好的辦法。
  • 這兩種方法的共同問題就是只考慮了某個有特點的數據，而沒有考慮類內數據的整體特點。
  • 平均距離或者平均鏈接，相對能得到合適一點的結果。平均鏈接的一個變種就是取兩兩距離的中值(中位數)，與取均值相比更加能夠解除個別偏離樣本對結果的干擾。

參考2 (from 百度文庫)：

  • 單鏈技術可以處理非橢圓形狀的簇，但對雜訊和離群點很敏感。
  • 全鏈接對雜訊和離群點不敏感，但是可能使大的簇破裂，偏好球型簇。
  • 平均距離則是一種單鏈與全鏈的折中演算法。優勢是對雜訊和極端值影響小，局限是偏好球型簇。

個人理解：
前面兩處參考中主要考慮了不同距離規則下受雜訊數據的影響大小。下面不考慮雜訊的影響，考慮不同距離規則下發生合併的趨勢。
考慮下圖所示的情況，這裡有三個簇，有一個較大的簇和兩個較小的簇，假如把它們看成是中國戰國時期的秦、韓、魏三個諸侯國。
在這裡插入圖片描述
按最小距離演算法，首先合併秦韓；按最大距離演算法，首先合併韓魏。
因此，在最小距離演算法下，趨向於「遠交近攻」的蠶食原則，相互靠近的諸侯國首先進行合併；在最大距離演算法下，大國反應較為遲鈍，小國家之間優先進行「合縱抗衡」。

9.7 聚類結果中若每個簇都有一個凸包(包含簇樣本的凸多面體)，且這些凸包不相交，則稱為凸聚類。試析本章介紹的哪些聚類演算法只能產生凸聚類，哪些能產生非凸聚類。

答：

  1. 原型聚類(k-means，LVQ，高斯混合)所得結果均是凸聚類。比如，在k-means的聚類分介麵線段必然為某兩個聚類中心的中垂線，考慮如下圖所示情況：
    在這裡插入圖片描述
    聚類中心1的聚邊界為非凸包，右上角存在凹口。設其中一個分界線段為中心點1和中心點2的中垂線，那麼上圖中藍色區域距離中心點2更近，不可能劃分到聚類1。因此，k-means只能產生凸聚類。
  2. 基於密度聚類的DBSCAN能夠產生非凸聚類，因為該演算法基於樣本的分布密度對密度相連的樣本進行聚類，最終的結果通過聚類簇內的樣本來表徵。如果聚類簇內的樣本分布呈非凸型分布，那麼聚類結果也將是非凸包。比如，教材圖9.10中的情形。
  3. 屬於層次聚類的AGNES也能產生非凸聚類，比如如下情形中，採用dmin進行聚類，結果會出現包圍形式的非凸聚類結果。而當採用dmax時，趨向於產生凸型聚類結果。
    在這裡插入圖片描述

9.8 試設計一個聚類性能度量指標，並與9.2節中的指標比較。（暫缺）

答：

9.9 試設計一個能用於混合屬性的非度量距離。（暫缺）

答：

9.10 試設計一個能自動確定聚類數的改進k均值演算法，編程實現並在西瓜數據集4.0上運行。

答：

  • 方案1：最小化平方誤差(9.24)式來找到最佳k值。
    然而，該方案不可行，極限思考一下，令\(k=m，u_i=x_i, i=1~\sim m\)，此時的平方誤差為0，顯然，這樣就沒有什麼意義了。
  • 方案2：在(9.24)式的基礎上增加一個懲罰項作為最小化的目標函數，比如令\(J=erro+0.1*k=\sum_c\sum_{x\in C} |x-u_c|_2^2+0.1*k\)。
    試驗一下，結果如下：
    在這裡插入圖片描述
    觀察可見，平方誤差erro隨k值增加而遞減，以此來確定最佳k值將得到\(k_{best}=m\)的結論，與預期相符。
    而目標函數\(J\)曲線則呈上開口拋物線型，當k=4時，目標函數\(J\)最小。上圖為某一次運行結果，多次運行結果較穩定，該方案可行，可以用來確定最佳k值。
  • 方案3：以衡量聚類性能的內部指標DBI和DI指數來找到最佳k值，見(9.12)和(9.13)式，其中DBI指數越小越好，DI指數越大越好。在西瓜數據集4.0上試驗一下。結果如下：
    在這裡插入圖片描述
    觀察上圖，DBI曲線大體上呈拋物線型，k=4時最佳，多次運行較穩定，可用於確定最佳k值。
    DI指數大體上呈下開口拋物線，上圖所示最佳k值為6，但是多次運行時，DI曲線表現不穩定，不太適合用於確定最佳k值。
    從DBI和DI的表達式可見，DI的計算式中\(d_{min}(C_i,C_j)\)表示兩個簇之間的最短樣本距離，受個別雜訊樣本的影響較大，而DBI計算式中\(avg(C_i)\)和\(avg(C_j)\)表示簇內平均距離，是個統計平均量，受雜訊影響相對較小一些。

附：編程程式碼

習題9.4 (Python)

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 13 13:00:25 2020

@author: Administrator
"""

import numpy as np
import matplotlib.pyplot as plt

#========首先編寫兩個函數==============
# k_means：實現k-means聚類演算法
# points：用於繪圖的輔助函數，根據樣本點的坐標計算外圍凸多邊形的頂點坐標

def k_means(X,k,u0=None,MaxIte=30):
    # k-means聚類演算法
    # 輸入：
    #    X:樣本數據，m×n數組
    #    k:聚類簇數目
    #    u0:初始聚類中心
    #    MaxIte:最大迭代次數
    # 輸出:
    #    u:最終的聚類中心
    #    C:各個樣本所屬簇號
    #    erro:按(9.24)式計算的平方方差結果
    #    step:迭代次數
    
    m,n=X.shape   #樣本數和特徵數
    if u0==None:  #隨機選取k個樣本作為初始中心點
        u0=X[np.random.permutation(m)[:k],:] 
    u=u0.copy()
    step=0
    while True:
        step+=1
        u_old=u.copy()            #上一次迭代的中心點
        dist=np.zeros([m,k])      #存儲各個樣本到中心點的距離
        for kk in range(k):       #計算距離
            dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
        C=np.argmin(dist,axis=1)  #以距離最小的中心點索引號最為簇類號
        for kk in range(k):       #更新聚類中心
            u[kk]=X[C==kk,:].mean(axis=0)  
        if (u==u_old).all() or step>MaxIte:
            break                 #如果聚類中心無變化，或者超過最大
                                  #迭代次數，則退出迭代
    #=====計算平方誤差(9.24)
    erro=0
    for kk in range(k):
         erro+=((X[C==kk]-u[kk])**2).sum()
    return u,C,erro,step

def points(X,zoom=1.2):
    # 為了繪製出教材上那種凸聚類簇效果
    # 本函數用於計算凸多邊形的各個頂點坐標
    # 輸入:
    #     X:簇類樣本點的坐標
    #     zoom:縮放因子(最外圍樣本點向外擴展係數)
    # 輸出:
    #     ps:凸多邊形的頂點坐標
    
    X=X[:,0]+X[:,1]*1j  #將坐標複數化
    u=np.mean(X)        #聚類中心
    X=X-u               #原點移至聚類中心
    # 尋找凸多邊形的各個頂點坐標
    indexs=[]                         #存儲頂點坐標的索引號
    indexs.append(np.argmax(abs(X)))  #首先將距離距離中心最遠的點作為起始頂點  
    while True:
        p=X[indexs[-1]]     #當前最新確定的頂點
        X1=1E-5+(X-p)/(-p)  #以p點為原點，並且以u-p為x軸(角度為0)
        #上式中加1E-5的小正數是因為p點自己減去自己的坐標有時候會出現
        #(-0+0j)的情況，angle(-0+0j)=-180°,但是希望結果為0
        indexs.append(np.argmax(np.angle(X1)))  #新找到頂點
        if indexs[-1]==indexs[0]:               #如果這些頂點首尾相連了，則停止
            break
    # 將複數坐標還原成直角坐標
    ps=np.c_[np.real(X)[indexs],np.imag(X)[indexs]]
    ps=ps*zoom+[np.real(u),np.imag(u)]
    return ps

#================================================
#                  主程式
#================================================

#==============西瓜數據集4.0======================
D=np.array(
        [[0.697,0.460],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
         [0.403,0.237],[0.481,0.149],[0.437,0.211],[0.666,0.091],[0.243,0.267],
         [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.360,0.370],
         [0.593,0.042],[0.719,0.103],[0.359,0.188],[0.339,0.241],[0.282,0.257],
         [0.748,0.232],[0.714,0.346],[0.483,0.312],[0.478,0.437],[0.525,0.369],
         [0.751,0.489],[0.532,0.472],[0.473,0.376],[0.725,0.445],[0.446,0.459]])
m=D.shape[0]
#=============繪製樣本數據點及其編號===============
plt.figure()
plt.scatter(D[:,0],D[:,1],marker='o',s=250,c='',edgecolor='k')
for i in range(m):
    plt.text(D[i,0],D[i,1],str(i),
             horizontalalignment='center',
             verticalalignment='center')
plt.show()

#選取三種不同的初始聚類中心點
centers=np.array([[5,7,17,18],     #相互靠近的點
                  [29,15,10,25],   #分散而外周的點
                  [27,17,12,21]])  #分散而中間的點

#======運行k-means聚類，設置三組不同的k值、三組不同初始中心點=======
for i,k in enumerate([2,3,4]):    #不同的k值
    for j in range(3):            #3次不同初始值
        #=====k-means演算法
        u0=D[centers[j][:k],:]
        u,C,erro,step=k_means(D,k,u0)
        #=====畫圖======
        plt.subplot(3,3,i*3+j+1)
        plt.axis('equal')
        plt.title('k=%d,step=%d,erro=%.2f'%(k,step,erro))
        # 畫樣本點
        plt.scatter(D[:,0],D[:,1],c='k',s=1)
        # 畫聚類中心
        plt.scatter(u0[:,0],u0[:,1],marker='o',c='',s=50,edgecolors='b')
        plt.scatter(u[:,0],u[:,1],marker='+',c='r',s=50)#,c='',s=80,edgecolors='g')
        for kk in range(k):
            plt.annotate('',xy=u[kk],xytext=u0[kk],arrowprops=dict(arrowstyle='->'),ha='center')
        # 畫聚類邊界
        for kk in range(k):
            ps=points(D[C==kk])
            plt.plot(ps[:,0],ps[:,1],'--r',linewidth=1)
plt.show()

習題9.10 (Python)

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 13 13:00:25 2020

@author: Administrator
"""

import numpy as np
import matplotlib.pyplot as plt

#========首先編寫兩個函數==============
# k_means：實現k-means聚類演算法
# points：用於繪圖的輔助函數，根據樣本點的坐標計算外圍凸多邊形的頂點坐標

def k_means(X,k,u0=None,MaxIte=30):
    # k-means聚類演算法
    # 輸入：
    #    X:樣本數據，m×n數組
    #    k:聚類簇數目
    #    u0:初始聚類中心
    #    MaxIte:最大迭代次數
    # 輸出:
    #    u:最終的聚類中心
    #    C:各個樣本所屬簇號
    #    erro:按(9.24)式計算的平方方差結果
    #    step:迭代次數
    
    m,n=X.shape   #樣本數和特徵數
    if u0==None:  #隨機選取k個樣本作為初始中心點
        u0=X[np.random.permutation(m)[:k],:] 
    u=u0.copy()
    step=0
    #print('D維度',X.shape,'u=',u)
    while True:
        
        step+=1
        u_old=u.copy()            #上一次迭代的中心點
        dist=np.zeros([m,k])      #存儲各個樣本到中心點的距離
        for kk in range(k):       #計算距離
            dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
        C=np.argmin(dist,axis=1)  #以距離最小的中心點索引號最為簇類號
        #print(u,C)
        for kk in range(k):       #更新聚類中心
            u[kk]=X[C==kk,:].mean(axis=0)  
        if (u==u_old).all():      #如果聚類中心無變化，則退出迭代
            break
        if step>MaxIte:           #如果超過最大迭代次數，則退出迭代
            print('超過最大迭代次數')
            break
    #=====計算平方誤差(9.24)
    erro=0
    for kk in range(k):
         erro+=((X[C==kk]-u[kk])**2).sum()
    return u,C,erro,step


def Inner_Index(X,u):
    #計算衡量聚類性能的內部指標DBI和DI指數
    m,n=X.shape
    k=u.shape[0]
    dist=np.zeros([m,k])      #存儲各個樣本到中心點的距離
    for kk in range(k):       #計算距離
        dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
    C=np.argmin(dist,axis=1)  #以距離最小的中心點索引號最為簇類號
    #===計算avg和diam
    avg=[]                    #存儲簇內平均距離
    diam=[]                   #存儲簇內最遠距離
    for kk in range(k):
        Xkk=X[C==kk]          #第kk個聚類簇內的樣本點
        mk=len(Xkk)
        if mk<=1:
            avg.append(0)
            diam.append(0)
            continue
        dist=[]               #存儲簇內兩個不同樣本之間的距離
        for i in range(mk):
            for j in range(i+1,mk):
                dist.append(np.sqrt(((Xkk[i]-Xkk[j])**2).sum()))
        avg.append(np.mean(dist))
        #print('dist=',dist,'mk=',mk)
        diam.append(np.max(dist))
    #===計算DB指數
    DBIij=np.zeros([k,k])    #存儲 [avg(Ci)+avg(Cj)]/dcen(ui,uj)
    for i in range(k):
        for j in range(i+1,k):
            DBIij[i,j]=(avg[i]+avg[j])/np.sqrt(((u[i]-u[j])**2).sum())
            DBIij[j,i]=DBIij[i,j]
    DBI=np.mean(np.max(DBIij,axis=1))
    #===計算Dunn指數
    dmin=[]    #存儲Ci和Cj之間的最短距離
    for i in range(k):
        for j in range(i+1,k):
            Xi=X[C==i]
            Xj=X[C==j]
            dmin_ij=[]
            for xi in Xi:
                for xj in Xj:
                    dmin_ij.append(np.sqrt(((xi-xj)**2).sum()))
            dmin.append(min(dmin_ij))
    DI=min(dmin)/max(diam)
    
    return DBI,DI


#================================================
#                  主程式
#================================================

#==============西瓜數據集4.0======================
D=np.array(
        [[0.697,0.460],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
         [0.403,0.237],[0.481,0.149],[0.437,0.211],[0.666,0.091],[0.243,0.267],
         [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.360,0.370],
         [0.593,0.042],[0.719,0.103],[0.359,0.188],[0.339,0.241],[0.282,0.257],
         [0.748,0.232],[0.714,0.346],[0.483,0.312],[0.478,0.437],[0.525,0.369],
         [0.751,0.489],[0.532,0.472],[0.473,0.376],[0.725,0.445],[0.446,0.459]])
m=D.shape[0]

#======考察平方誤差(9.24式)隨k值的變化規律=======
# 結果如預期，平方誤差隨k值增加而減小
# 如果加上一個懲罰項 0.1*k，出現拋物線
erros=[]
ks=range(2,8+1)         #考察k=2~8
for k in ks:    
    erro=[]
    for i in range(5):  #對於每個k值，初始5次聚類，取其中最佳者
        u,C,erro_i,step=k_means(D,k)
        erro.append(erro_i)
    erros.append(min(erro))
plt.figure()
plt.plot(ks,erros,'o-')
plt.plot(ks,0.1*np.array(ks)+np.array(erros),'o-')
plt.xlabel('k')
plt.legend(['erro','erro+0.1*k'])
plt.show()


#======考察DBI評分(9.12式)隨k值的變化規律=======
DBIs=[]
DIs=[]
ks=range(2,8+1)         #考察k=2~8
for k in ks:    
    DBI=[]
    DI=[]
    for i in range(5):  #對於每個k值，初始5次聚類，取其中最佳者
        u,C,erro,step=k_means(D,k)
        DBi,Di=Inner_Index(D,u)
        DBI.append(DBi)
        DI.append(Di)
    DBIs.append(min(DBI))
    DIs.append(max(DI))
# 繪製計算結果
ax1=plt.figure().add_subplot(111)
ax2=ax1.twinx()    #雙y坐標

ax1.plot(ks,DBIs,'o-k',label='DBI')
ax2.plot(ks,DIs,'o-r',label='DI')

ax1.set_xlabel('k')
ax1.set_ylabel('DBI')
ax2.set_ylabel('DI')

ax1.legend(loc=(0.5,.9))
ax2.legend(loc=(0.5,.8))

plt.show()
Tags:

VirMach 便宜 VPS

QNews

QNews