（數據科學學習手札08）系統聚類法的Python源碼實現（與Python，R自帶方法進行比較）

聚類分析是數據挖掘方法中應用非常廣泛的一項，而聚類分析根據其大體方法的不同又分為系統聚類和快速聚類，其中系統聚類的優點是可以很直觀的得到聚類數不同時具體類中包括了哪些樣本，而Python和R中都有直接用來聚類分析的函數，但是要想掌握一種方法就得深刻地理解它的思想，因此自己從最底層開始編寫代碼來實現這個過程是最好的學習方法，所以本篇前半段是筆者自己寫的代碼，如有不細致的地方，望指出。

一、僅使用numpy包進行系統聚類的實現：

'''以重心法為距離選擇方法搭建的系統聚類算法原型'''
# @Feffery 
# @說明：目前僅支持維度為2，重心法的情況

import numpy as np
import time

price = [1.1,1.2,1.3,1.4,10,11,20,21,33,34]
increase = [1 for i in range(10)]
data = np.array([price,increase],dtype='float32')


class Myhcluster():

    def __init__(self):
        print('開始進行系統聚類')
    '''系統聚類法的啟動函數，有輸入變量和距離計算方法兩個輸入參數'''
    def prepare(self,data,method='zx'):
        if method == 'zx':
            self.zx(data)
    '''重心法進行系統聚類'''
    def zx(self,data):
        token = len(data[0,:])
        flu_data = data.copy()
        classfier =[[] for i in range(len(data[1,]))]
        LSdist = np.array([0 for i in range(token ** 2)], dtype='float32').reshape([len(data[0, :]), token])
        index = 0
        while token > 1:
            '''計算距離矩陣'''
            for i in range(len(data[0,:])):
                for j in range(len(data[0,:])):
                        LSdist[i,j] = round(((flu_data[0,i]-flu_data[0,j])**2+(flu_data[1,i]-flu_data[1,j])**2)**0.5,4)

            '''將距離矩陣中的0元素替換為NAN'''
            for i in range(len(data[0,:])):
                for j in range(len(data[0,:])):
                    if LSdist[i,j] == 0:
                        LSdist[i,j] = np.nan

            '''保存該次系統聚類中最短距離對應的兩個樣本的標號'''
            T = set([np.argwhere(LSdist == np.nanmin(LSdist))[0,0],np.argwhere(LSdist == np.nanmin(LSdist))[0,1]])
            TT = [i for i in T]

            '''針對該次聚類情況進行產生新子類亦或是歸入舊子類的選擇'''
            RQ = TT
            for x in range(len(classfier)):
                if classfier[0] == []:#判斷是否為n個樣本中第一次迭代產生新類
                    classfier[0] = TT
                    index = 0
                    break
                elif classfier[-2] != []:#判斷是否已在理論最大歸類次數前完成所有樣品的聚類
                    print('最后一次分類，獲得由樣本{}組成的新類'.format([__ for __ in range(len(data[1,]))]))
                    return 0
                elif TT[0] in classfier[x] or TT[1] in classfier[x]:
                    if classfier[x+1]==[]:
                        classfier[x+1] = list(set(classfier[x]).union(set(RQ)))
                        index = x+1
                        break
                    else:
                        RQ = list(set(classfier[x]).union(set(RQ)))
                        classfier[len(data[1,])-token] = RQ
                        continue
                elif x == len(data[1,])-1:
                    classfier[len(data[0,:])-token] = TT
                    index = len(data[0,:])-token
            print('第{}次分類，獲得由樣本{}組成的新類'.format(str(len(data[0,:])-token+1),set(classfier[index])))
            #求得重心並對原數據進行覆蓋
            for k in set(classfier[index]):
                flu_data[0,k] = np.mean([data[0,_] for _ in set(classfier[index])])
                flu_data[1,k] = np.mean([data[1, _] for _ in set(classfier[index])])
            token -= 1


a = time.clock()
dd = Myhcluster()#進行算法封裝的類的傳遞
dd.prepare(data)#調用類中的系統聚類法（默認重心法）
print('自己編寫的系統聚類算法使用了'+str(round(time.clock()-a,3))+'秒')

 

與Scipy中系統聚類方法進行比較：

'''與Scipy中自帶的層次聚類方法進行比較'''
import scipy.cluster.hierarchy as sch
import numpy as np

a = time.clock()
disMat = sch.distance.pdist(data.T,'euclidean')

Z=sch.linkage(disMat,method='average')

sch.dendrogram(Z)
print('Scipy中的系統聚類算法用了'+str(round(time.clock()-a,3))+'秒')

 

與R自帶系統聚類算法進行比較：

> #系統聚類法的R實現
> rm(list=ls())
> a <- Sys.time()
> price <- c(1.1,1.2,1.3,1.4,10,11,20,21,33,34)
> increase <- rep(1,10)
> data <- data.frame(price,increase)#生成樣本數據框
> d <- dist(data)#創建樣本距離陣
> hc <- hclust(d,'centroid')#用重心法進行系統聚類
> cbind(hc$merge, hc$height)#展示分類過程
[,1] [,2] [,3]
[1,] -1 -2 0.10000
[2,] -3 -4 0.10000
[3,] 1 2 0.15000
[4,] -5 -6 1.00000
[5,] -7 -8 1.00000
[6,] -9 -10 1.00000
[7,] 3 4 8.93750
[8,] 5 6 12.50000
[9,] 7 8 17.18056
> Sys.time() - a
Time difference of 0.007000923 secs
> plot(hc)#繪制層次聚類圖