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

時間  2020-09-26
標籤 西瓜 習題 简体版
原文   原文鏈接

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

目錄

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}| \]

證實：python

  • 非負性、同一性、對稱性
    對於全部的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）不等式的證實過程以下：app

    首先須要證實一個引理和Holder不等式，而後再證實閔可夫斯基不等式。dom

    引理 \(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()
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程