《資料分析實戰-托馬茲.卓巴斯》讀書筆記第4章-聚類技巧（K均值、BIRCH、DBSCAN）

邀月發表於2020-03-16
筆記聚類

python資料分析個人學習讀書筆記-目錄索引

 

第4章　解釋了多種聚類模型；從最常見的k均值演算法開始，一直到高階的BIRCH演算法和DBSCAN演算法。

本章會介紹一些技術，幫助對一個銀行營銷電話的資料進行聚類。將學習以下主題：
·評估聚類方法的表現
·用k均值演算法聚類資料
·為k均值演算法找到最優的聚類數
·使用mean shift聚類模型發現聚類
·使用c均值構建模糊聚類模型
·使用層次模型聚類資料
·使用DBSCAN和BIRCH演算法發現潛在的訂閱者

4.1導論

分類問題中，我們知道每個觀測值的歸類（經常稱作監督學習或有老師的學習），聚類問題則不然，不需要標籤就能在資料中發現模式（稱作無監督學習）。
聚類方法根據觀測值之間的相似度將其歸到不同的桶中。這種分析在探索階段很有用，這個階段分析師想知道資料中是否天然存在某種模式。

 4.2評估聚類方法的表現

沒有真實的標籤，我們無法使用前一章介紹的指標。在本技巧中，我們會介紹三種幫我們評估聚類方法有效性的指標：Davis-Bouldin、偽F（有時也稱作Calinsk-Harabasz）和Silhouette值是內部評價指標。如果知道真實的標籤，我們就可以使用很多指標，比如調整後的Rand指數、同質性或是完整性。
參考scikit文件中關於聚類方法的部分，對聚類方法的多種外部評價指標有一個更深的瞭解：

http://scikit-learn.org/stable/modules/clustering.html#clustering-performance-evaluation；

內部聚類驗證方法的列表，參考http://datamining.rutgers.edu/publication/internalmeasures.pdf。
本示例需要裝好pandas、NumPy和Scikit。

前面提到的三個內部評價指標，Scikit只實現了Silhouette值；為了寫這本書，原作者實現了另兩個。要評估你的聚類模型的表現，你可以使用helper.py檔案中的printCluster Summary（...）方法：

def printClustersSummary(data, labels, centroids):
    '''
        Helper method to automate models assessment
    '''
    print('Pseudo_F: ', pseudo_F(data, labels, centroids))
    print('Davis-Bouldin: ',
        davis_bouldin(data, labels, centroids))
    print('Silhouette score: ',
        mt.silhouette_score(data, np.array(labels),
            metric='euclidean'))

原理:我們介紹第一個指標，偽F分值。Calinski-Harabasz指標是在模型上啟發式地得到聚類的數目，其最大值代表最優的聚類方式。
偽F指數計算方式，是用每個聚類的中心到整個資料集的幾何中心的平方距離，與聚類數減1的比值，再乘以每個聚類中觀測值的數目。這個數字再除以一個比值，這個比值是用每個點和聚類中心的平方距離，除以所有觀測值數與聚類數之差得到的。

 1 def pseudo_F(X, labels, centroids):
 2     '''
 3         The pseudo F statistic :
 4         pseudo F = [( [(T - PG)/(G - 1)])/( [(PG)/(n - G)])]
 5         The pseudo F statistic was suggested by
 6         Calinski and Harabasz (1974) in
 7         Calinski, T. and J. Harabasz. 1974.
 8             A dendrite method for cluster analysis.
 9             Commun. Stat. 3: 1-27.
10             http://dx.doi.org/10.1080/03610927408827101
11 
12         We borrowed this code from
13         https://github.com/scampion/scikit-learn/blob/master/
14         scikits/learn/cluster/__init__.py
15 
16         However, it had an error so we altered how B is
17         calculated.
18     '''
19     center = np.mean(X,axis=0)
20     u, count = np.unique(labels, return_counts=True)
21 
22     B = np.sum([count[i] * ((cluster - center)**2)
23         for i, cluster in enumerate(centroids)])
24 
25     X = X.as_matrix()
26     W = np.sum([(x - centroids[labels[i]])**2
27                  for i, x in enumerate(X)])
28 
29     k = len(centroids)
30     n = len(X)
31 
32     return (B / (k-1)) / (W / (n-k))

首先，我們計算整個資料集中心的幾何座標。從代數上講，這不過是計算每一列的平均值。使用NumPy的.unique（...），然後統計每個聚類中觀測值的數目。
往.unique（...）方法傳return_counts=True時，不僅會返回標籤向量中的去重值列表u，也會返回去重後每個值的個數c。
然後，我們計算每個聚類的中心到資料集中心的平方距離。使用.enumerate（...）方法遍歷中心列表中的每個元素，計算每個和資料集中心的平方差，再乘以聚類中觀測值的個數：c[i]。如名字所示，NumPy的.sum（...）方法，計算列表中所有元素的和。這個方法返回的是Calinski-Harabasz指數。
與偽F指數相反，Davis-Bouldin指標衡量的是最壞情況下聚類之間的異質性與聚類內部的同質性。也就是說，目標要找出使得這個指標最小的聚類數目。在後面的4.4節，我們會開發一種方法，透過最小化Davis-Bouldin指標來找到最優的聚類數目。
參考MathWorks（Matlab開發商）關於Davis-Bouldin指標的講解：http://au.mathworks.com/help/stats/clustering.evaluation.daviesbouldinevaluation-class.html。

 1  def davis_bouldin(X, labels, centroids):
 2     '''
 3         The Davis-Bouldin statistic is an internal evaluation
 4         scheme for evaluating clustering algorithms. It
 5         encompasses the inter-cluster heterogeneity and
 6         intra-cluster homogeneity in one metric.
 7 
 8         The measure was introduced by
 9         Davis, D.L. and Bouldin, D.W. in 1979.
10             A Cluster Separation Measure
11             IEEE Transactions on Pattern Analysis and
12             Machine Intelligence, PAMI-1: 2, 224--227
13 
14             http://dx.doi.org/10.1109/TPAMI.1979.4766909
15     '''
16     distance = np.array([
17         np.sqrt(np.sum((x - centroids[labels[i]])**2))
18         for i, x in enumerate(X.as_matrix())])
19 
20     u, count = np.unique(labels, return_counts=True)
21 
22     Si = []
23 
24     for i, group in enumerate(u):
25         Si.append(distance[labels == group].sum() / count[i])
26 
27     Mij = []
28 
29     for centroid in centroids:
30         Mij.append([
31             np.sqrt(np.sum((centroid - x)**2))
32             for x in centroids])
33 
34     Rij = []
35     for i in range(len(centroids)):
36         Rij.append([
37             0 if i == j
38             else (Si[i] + Si[j]) / Mij[i][j]
39             for j in range(len(centroids))])
40 
41     Di = [np.max(elem) for elem in Rij]
42 
43     return np.array(Di).sum() / len(centroids)

首先，我們計算每個觀測值和所屬聚類的中心之間的幾何距離，並統計每個聚類中觀測值的數目。
Si衡量的是聚類內部的同質性，即是每個觀測值到聚類中心的平均距離。Mij計算各個聚類中心之間的幾何距離，量化了聚類之間的異質性。Rij衡量的是兩個聚類之間切分得好不好，Di選取了最差情況下的切分。Davis-Bouldin指標是Di的平均值。
Silhouette分值（指數）是又一個聚類內部評價指標。我們沒有實現計算程式（因為Scikit已經提供了實現），但我們簡要描述了它的計算過程與衡量目標。一個聚類，其Silhouette分值衡量的是聚類中每個觀測值與其他觀測值之間的平均距離，將這個值與每個聚類內點與點之間的平均距離相比。理論上，這個指標的取值範圍在-1到1之間；-1意味著所有的觀測值都在不合適的聚類中，實際上，更合適的聚類可能就是臨近的某個。儘管理論上可行，但實際使用中，你更可能完全不考慮負的Silhouette分值。光譜的另一端，1代表完美的切分，所有的觀測值都落在合適的聚類中。這個值如果是0，代表著聚類之間一個完美重疊，意味著每個聚類都有應該歸屬於其他聚類的觀測值，並且這些觀測值的數目是相等的。
參考:對這些演算法（以及其他演算法）更深層次的描述，參考https://cran.r-project.org/web/packages/clusterCrit/vignettes/clusterCrit.pdf。

4.3用k均值演算法聚類資料
k均值聚類演算法可能是在聚類向量化資料方面最為知名的資料探勘技術了。它基於觀測察值之間的相似度，將它們分到各個聚類；決定因素就是觀測值和最近一個聚類中心的歐幾里得距離。
裝好pandas和Scikit。
Scikit在其聚類子模組中提供了多種聚類模型。這裡，我們使用.KMeans（...）來估算聚類模型（clustering_kmeans.py檔案）：

def findClusters_kmeans(data):
    '''
        Cluster data using k-means
    '''
    # create the classifier object
    kmeans = cl.KMeans(
        n_clusters=4,
        n_jobs=-1,
        verbose=0,
        n_init=30
    )

    # fit the data
    return kmeans.fit(da */ta)   

原理:和前一章的技巧一樣，我們從讀取資料開始，選擇我們想用來區別觀測值的特徵：

# the file name of the dataset
r_filename = '../../Data/Chapter04/bank_contacts.csv'

# read the data
csv_read = pd.read_csv(r_filename)

# select variables
selected = csv_read[['n_duration','n_nr_employed',
        'prev_ctc_outcome_success','n_euribor3m',
        'n_cons_conf_idx','n_age','month_oct',
        'n_cons_price_idx','edu_university_degree','n_pdays',
        'dow_mon','job_student','job_technician',
        'job_housemaid','edu_basic_6y']]

Tips:

/* The method findClusters_kmeans took 2.32 sec to run.
D:\Java2018\practicalDataAnalysis\helper.py:142: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  X = X.as_matrix()
Pseudo_F:  11515.72135543927
D:\Java2018\practicalDataAnalysis\helper.py:168: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  for i, x in enumerate(X.as_matrix())])
Davis-Bouldin:  1.2627990546726073
Silhouette score:  0.3198322783627837
 */

解決方案：修改X.as_matrix()為X.values

/*
The method findClusters_kmeans took 2.32 sec to run.
Pseudo_F:  11349.586839983096
Davis-Bouldin:  1.2608899499519972
 */


我們使用了與前一章相同的資料集，並只使用在分類過程中最具描述力的特徵。另外，我們仍然使用@hlp.timeit裝飾器來衡量估算模型的速度。
Scikit的.KMeans（...）方法有多個引數。引數n_clusters定義了期待資料中要有多少聚類，以及要返回多少聚類。在本書4.4節中，我們將寫一個迴圈方法，為k均值演算法找出最優的聚類數。
n_jobs引數是你機器上並行執行的任務數；指定-1意味著將並行執行的數目指定為處理器核的數目。你可以顯式傳入一個整數，比如8，來指定任務數目。
verbose引數控制著估算階段你能看到多少資訊；設定為1會列印出所有關於估算的細節。
n_init引數控制著要估算多少模型。k均值演算法的每一次執行，會隨機選擇每個聚類的中心，迴圈調整這些中心，直到模型聚類之間的隔離程度和聚類內部的相似程度達到最好。.KMeans（...）方法以不同的初始條件（為各中心隨機化初始種子），構建n_init引數指定的數目的模型，然後選擇準則表現最好的一個。
準則衡量的是聚類內部的差異。它是聚類內部的平方和。我們將聚類內部的平方和與平方總和的比率稱為解釋準則。

Scikit當然不是估算k均值模型的唯一方法；我們也可以使用SciPy（clustering_kmeans_alternative.py檔案）：

 1 def findClusters_kmeans(data):
 2     '''
 3         Cluster data using k-means
 4     '''
 5     # whiten the observations
 6     data_w = vq.whiten(data)
 7 
 8     # create the classifier object
 9     kmeans, labels = vq.kmeans2(
10         data_w,
11         k=4,
12         iter=30
13     )
14 
15     # fit the data
16     return kmeans, labels
17 
18 # the file name of the dataset
19 r_filename = '../../Data/Chapter04/bank_contacts.csv'
20 
21 # read the data
22 csv_read = pd.read_csv(r_filename)
23 
24 # select variables
25 selected = csv_read[['n_duration','n_nr_employed',
26         'prev_ctc_outcome_success','n_euribor3m',
27         'n_cons_conf_idx','n_age','month_oct',
28         'n_cons_price_idx','edu_university_degree','n_pdays',
29         'dow_mon','job_student','job_technician',
30         'job_housemaid','edu_basic_6y']]
31 
32 # cluster the data
33 #centroids, labels = findClusters_kmeans(selected.as_matrix())
34 centroids, labels = findClusters_kmeans(selected.values)
35 
36 hlp.printClustersSummary(selected, labels, centroids)

與使用Scikit時不同，使用SciPy時，我們需要先洗白資料。洗白類似於將資料標準化（參考本書1.14節），只不過它不刪除均值；.whiten（...）將所有特徵的方差統一為1。
傳給.kmeans2（...）方法的第一個引數是已清洗的資料集。k引數指定了將資料應用到多少個聚類上。我們的例子中，.kmeans2（...）最多執行30個迴圈（iter引數）；如果那時方法還沒有收斂，其也會停止，並返回最後一次的估算結果。
參考,MLPy也提供了一個估算k均值模型的方法：http://mlpy.sourceforge.net/docs/3.5/cluster.html#k-means。

4.4為k均值演算法找到最優的聚類數

很多時候，你不知道資料集中該有多少聚類。對於二維或三維資料，你可以畫一畫資料集，肉眼識別出聚類。然而，對於高維資料，在一張圖表中繪製都不太可能，所以這種方法就不太可行了。
在本技巧中，我們會展示如何為k均值模型找到最優的聚類數。我們將使用Davis-Bouldin指標來評估不同聚類數下k均值模型的表現。目標是當這個指標達到一個最小值時停止。
需裝好pandas、NumPy和Scikit。

為了找到最優的聚類數，我們編寫了findOptimalClusterNumber（...）方法。總體而言，估算k均值模型的演算法並沒有改變——只不過用findOptimalClusterNumber（...）方法替換了前面使用的findClusters_kmeans（...）方法（clustering_kmeans_search.py檔案）：

# find the optimal number of clusters; that is, the number of
# clusters that minimizes the Davis-Bouldin criterion
optimal_n_clusters = findOptimalClusterNumber(selected)

findOptimalClusterNumber（...）方法定義如下：

 1 def findOptimalClusterNumber(
 2         data, 
 3         keep_going = 1, 
 4         max_iter = 30
 5     ):
 6     '''
 7         A method that iteratively searches for the 
 8         number of clusters that minimizes the Davis-Bouldin
 9         criterion
10     '''
11     # the object to hold measures 
12     measures = [666]
13 
14     # starting point
15     n_clusters = 2 
16     
17     # counter for the number of iterations past the local 
18     # minimum
19     #超過區域性最小值的迴圈次數的計數器
20     keep_going_cnt = 0
21     stop = False   # flag for the algorithm stop
22     
23     def checkMinimum(keep_going):
24         '''
25             A method to check if minimum found
26         '''
27         global keep_going_cnt # access global counter
28 
29         # if the new measure is greater than for one of the 
30         # previous runs
31         if measures[-1] > np.min(measures[:-1]):
32             # increase the counter
33             keep_going_cnt += 1
34 
35             # if the counter is bigger than allowed 
36             if keep_going_cnt > keep_going:
37                 # the minimum is found
38                 return True
39         # else, reset the counter and return False
40         else:
41             keep_going_cnt = 0
42 
43         return False
44 
45     # main loop 
46     # loop until minimum found or maximum iterations（迴圈） reached
47     while not stop and n_clusters < (max_iter + 2):
48         # cluster the data
49         cluster = findClusters_kmeans(data, n_clusters)
50 
51         # assess the clusters effectiveness
52         labels = cluster.labels_
53         centroids = cluster.cluster_centers_
54 
55         # store the measures（指標）   
56         measures.append(
57             hlp.davis_bouldin(data,labels, centroids)
58         )
59 
60         # check if minimum found
61         stop = checkMinimum(keep_going)
62 
63         # increase the iteration
64         n_clusters += 1
65 
66     # once found -- return the index of the minimum
67     return measures.index(np.min(measures)) + 1

原理:使用findOptimalClusterNumber（...）方法時，傳入的第一個引數是資料集：顯然，沒有資料時，什麼模型也估算不出。keep_going引數指的是，如果當前模型的Davis-Bouldin指標比已有的最小值要大，那還要尋找多少個迴圈。透過這種方式，我們可以（一定程度上）避免找到區域性最小值而非全域性最小值的情況：
邀月工作室

max_iter引數指定了最多構建多少個模型；我們的方法以n_clusters=2構建k均值模型開始，然後迴圈執行，直到找到一個全域性最小值或是達到了最大的迴圈次數。
我們的findOptimalClusterNumber（...）方法以定義執行的引數開始。
measures列表用於儲存估算模型的Davis-Bouldin指標；由於我們的目標是找到最小的Davis-Bouldin指標，我們指定了一個大數，以免這個值成了我們的最小值。
n_clusters定義了起始點：第一個k均值模型將資料分成兩個聚類。
keep_going_cnt用來指定Davis-Bouldin指標比之前的最小值要大時，還要走多少個迴圈。以前面展示的圖為例，我們的方法不會終止在某個區域性最小值；儘管迴圈10和迴圈13的Davis-Bouldin指標（分別）比迴圈9和迴圈12小，但是聚類數為11和14的模型有更低的值。在這個例子中，我們在聚類數為15和16時終止，此時Davis-Bouldin指標更大。
標識變數stop用於控制主迴圈的執行。我們的while迴圈直到找到最小值或者達到了最大迴圈次數。
一個人讀程式碼時，看到當不為stop並且n_clusters<（max_iter+2）時while迴圈不會終止，他能輕鬆地將程式碼翻譯成自然語言，也能體會到Python語法之美。我覺得，這提高了程式碼的可讀性與可維護性，同時也使程式碼更準確。
迴圈以估算預定義聚類數的模型開始：

# cluster the data
  cluster = findClusters_kmeans(data, n_clusters)

估算出模型後，我們得到估算的標籤和中心，使用.davis_bouldin（...）方法計算Davis-Bouldin指標，如本章第一個技巧中那樣。這個指標透過append（...）方法附到指標列表。
現在，我們要檢查新估算的模型，其Davis-Bouldin指標是否比之前估算的模型都小。我們使用checkMinimum（...）方法。這個方法只能由findOptimalClusterNumber（...）方法訪問：

 1 def checkMinimum(keep_going):
 2     '''
 3         A method to check if minimum found
 4     '''
 5     global keep_going_cnt # access global counter
 6 
 7     # if the new measure is greater than for one of the 
 8     # previous runs
 9     if measures[-1] > np.min(measures[:-1]):
10         # increase the counter
11         keep_going_cnt += 1
12 
13         # if the counter is bigger than allowed 
14         if keep_going_cnt > keep_going:
15             # the minimum is found
16             return True
17     # else, reset the counter and return False
18     else:
19         keep_going_cnt = 0
20 
21     return False

首先，我們將keep_going_cnt定義為一個全域性變數。這樣我們不需要將變數傳入check Minimum（...）方法，並且可以全域性追蹤計數器。在我看來，這樣有益於程式碼的可讀性。
然後，我們將新近估算的模型與已有的最小值比較；我們使用NumPy的.min（...）方法得到最小值。如果當前的Davis-Bouldin指標比已有的最小值大，我們就給keep_going_cnt計數器增1，並檢查是否超過了keep_going這個迴圈次數限制——超過了就返回True。返回True意味著我們找到了全域性最小值（在keep_going限制的前提下）。然而，如果新估算的Davis-Bouldin指標更小，我們先將keep_going_cnt計數器重置為0，返回False。如果沒有超過keep_going的限制，我們也返回False。
知道了我們找沒找到全域性最小值之後，我們將n_clusters增1。
當checkMinimum（...）方法返回True或者n_clusters超過max_iter+2時迴圈中止。我們的演算法一開始用兩個聚類來估算k均值模型，所以max_iter要加上2。
一旦中止while迴圈，我們返回找到的最小值的索引加1這個值。要加上一個1，因為n_
一旦中止while迴圈，我們返回找到的最小值的索引加1這個值。要加上一個1，因為n_cluster=2時Davis-Bouldin指標的索引是1。
列表的索引值從0開始。
現在我們可以用最優的聚類數來估算模型，並列印指標：

# cluster the data
cluster = findClusters_kmeans(selected, optimal_n_clusters)

# assess the clusters effectiveness評估聚類的有效性
labels = cluster.labels_
centroids = cluster.cluster_centers_

hlp.printClustersSummary(selected, labels, centroids)
/* 

The method findClusters_kmeans took 1.97 sec to run.
The method findClusters_kmeans took 0.56 sec to run.
The method findClusters_kmeans took 0.76 sec to run.
The method findClusters_kmeans took 0.90 sec to run.
The method findClusters_kmeans took 1.01 sec to run.
The method findClusters_kmeans took 1.16 sec to run.
The method findClusters_kmeans took 1.31 sec to run.
The method findClusters_kmeans took 1.34 sec to run.
The method findClusters_kmeans took 1.55 sec to run.
The method findClusters_kmeans took 1.69 sec to run.
The method findClusters_kmeans took 1.82 sec to run.
The method findClusters_kmeans took 2.08 sec to run.
The method findClusters_kmeans took 2.11 sec to run.
The method findClusters_kmeans took 2.40 sec to run.
The method findClusters_kmeans took 2.49 sec to run.
The method findClusters_kmeans took 3.03 sec to run.
The method findClusters_kmeans took 2.98 sec to run.
The method findClusters_kmeans took 3.05 sec to run.
Optimal number of clusters: 17
The method findClusters_kmeans took 2.96 sec to run.
Pseudo_F:  11493.774263351304
Davis-Bouldin:  0.9516597767023216
*/


之前提到，要繪製資料的聚類並不容易，不過，探索資料視覺化時，畫成對的圖有時會很有幫助（clustering_kmeans_search_alternative.py檔案）：

 1 def plotInteractions(data, n_clusters):
 2     '''
 3         Plot the interactions between variables
 4     '''
 5     # cluster the data
 6     cluster = findClusters_kmeans(data, n_clusters)
 7 
 8     # append the labels to the dataset for ease of plotting
 9     data['clus'] = cluster.labels_
10 
11     # prepare the plot
12     ax = sns.pairplot(selected, hue='clus')
13 
14     # and save the figure
15     ax.savefig(
16         '../../Data/Chapter04/k_means_{0}_clusters.png' \
17         .format(n_clusters)
18     )

我們使用Seaborn和Matplotlib繪製互動。首先，我們用預先定義的n_clusters估算模型。然後，我們將聚類加到資料集中。使用Seaborn的.pairplot（...）方法，我們建立了一張圖表，遍歷每個特徵，畫出響應的互動。最後，我們儲存這張圖表。
14個聚類的結果（為了可讀性，限定3個特徵），看起來類似這樣：
邀月工作室
Tips:

/*
The method findClusters_kmeans took 3.43 sec to run.
D:\Java2018\practicalDataAnalysis\Codes\Chapter04\clustering_kmeans_search_alternative.py:37: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['clus'] = cluster.labels_
D:\tools\Python37\lib\site-packages\statsmodels\nonparametric\kde.py:487: RuntimeWarning: invalid value encountered in true_divide
  binned = fast_linbin(X, a, b, gridsize) / (delta * nobs)
D:\tools\Python37\lib\site-packages\statsmodels\nonparametric\kdetools.py:34: RuntimeWarning: invalid value encountered in double_scalars
  FAC1 = 2*(np.pi*bw/RANGE)**2
   */

關於DataFrame的loc方法：

/*    
   Compare these two access methods:
In [340]: dfmi['one']['second']
Out[340]:
0    b
1    f
2    j
3    n
Name: second, dtype: object
In [341]: dfmi.loc[:, ('one', 'second')]
Out[341]:
0    b
1    f
2    j
3    n
Name: (one, second), dtype: object
 */

即使維度很低，透過肉眼要看出資料集中應該有多少聚類也不容易。對於其他資料，建立這樣一張圖表也許會有好處。這個方法的問題在處理離散或類別的變數時會特別明顯。

使用mean shift聚類模型發現聚類

mean shift模型是一個類似於尋找中心（或最大密度）的方法。與k均值不同，這個方法不要求指定聚類數——這個模型根據資料中找到的密度中心的數目返回聚類的數目。
需裝好pandas和Scikit。
我們估算的開啟方式類似之前的模型——從讀入資料集和限制特徵數開始。然後我們用findClusters_meanShift（...）方法估算模型（clustering_meanShift.py檔案）：

 1 def findClusters_meanShift(data):
 2     '''
 3         Cluster data using Mean Shift method
 4         使用Mean Shift方法聚類資料
 5     '''
 6     bandwidth = cl.estimate_bandwidth(data,
 7         quantile=0.25, n_samples=500)
 8 
 9     # create the classifier object
10     meanShift = cl.MeanShift(
11         bandwidth=bandwidth,
12         bin_seeding=True
13     )
14 
15     # fit the data
16     return meanShift.fit(data)
/*
The method findClusters_meanShift took 20.53 sec to run.
Pseudo_F:  1512.7080324697126
Davis-Bouldin:  1.736359716302394
Silhouette score:  0.2172934930829158
*/

原理:我們首先估算頻寬（estimate_bandwidth（...））。頻寬會用在這個方法的RBF核中。
我們使用SVM分類技巧中的徑向基函式：參考本書3.5節。
estimate_bandwidth（...）方法至少要傳入資料作為引數。
由於方法中使用的演算法，隨著觀測值數目的增長，呈二次方複雜度，所以對於有很多觀測值的樣本，最好限制記錄數目；使用n_samples吧。
quantile引數決定了從什麼地方切分傳入.MeanShift（...）方法的核的樣本。
quantile取0.5就意味著中位數。
現在可以估算模型了。我們傳入bandwidth和（可選的）bin_seeding引數，構建模型物件。
如果bin_seeding置為True，核初始位置就設到了所有資料點的（使用bandwidth）離散化的群。這個引數置為True會加快演算法，因為要初始化的核種子變少了；隨著觀測值成百上千地增長，這個會帶來顯著的加速。
談到了估算的速度，這個方法比任何k均值都會慢得多（即便有數十個聚類）。在我們的機器上，估算會比k均值慢5到20秒：平均值在13秒左右。
要是對Mean Shift演算法的原理感興趣，作者推薦這篇文件：http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/TUZEL1/MeanShift.pdf。
參考本書3.5節。

4.6
使用c均值構建模糊聚類模型

k均值和mean shift聚類演算法將觀測值歸到截然不同的聚類：一個觀測值可以且只能被歸到一個相似樣本的聚類中。對於離散可分的資料集，這個也許是對的，但資料如果有重疊，將它們放進同一個桶裡就很困難了。畢竟，我們的世界並不是非黑即白，我們能看到無數色彩。
c均值聚類模型允許每個觀測值屬於不止一個聚類，各從屬關係有一個權重：對每個觀測值來說，它所屬的所有聚類對應的權重之和必須為1。
需要pandas和Scikit-Fuzzy模組。Scikit-Fuzzy模組一般不在Anaconda發行版中，所以你需要自行安裝。
>pip install Scikit-Fuzzy

/* Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
.......
Successfully built Scikit-Fuzzy
Installing collected packages: decorator, networkx, Scikit-Fuzzy
Successfully installed Scikit-Fuzzy-0.4.2 decorator-4.4.1 networkx-2.4
FINISHED
 */

Scikit已經提供估算c均值的演算法。呼叫方法和之前的略有不同（clustering_cmeans.py檔案）：

 1 import skfuzzy.cluster as cl
 2 import numpy as np
 3 
 4 @hlp.timeit
 5 def findClusters_cmeans(data):
 6     '''
 7         Cluster data using fuzzy c-means clustering
 8         algorithm
 9     '''
10     # create the classifier object
11     return cl.cmeans(
12         data,
13         c = 5,          # number of clusters 聚類數
14         m = 2,          # exponentiation factor 指數因子
15         
16         # stopping criteria
17         error = 0.01,
18         maxiter = 300
19     )

原理:如同之前的技巧，我們首先載入資料，選擇我們想用來估算模型的列。
然後呼叫findClusters_cmeans（...）方法估算模型。之前的方法都是建立一個未訓練的模型，然後用.fit（...）方法應用資料，與此不同，.cmeans（...）方法是建立時就用傳入的資料（及其他引數）估算模型了。
c引數指定了要應用多少個聚類，m引數指定了在每次迴圈時應用到成員關係函式的指數因子。
注意你要指定的聚類數，要是指定得太多，計算指標時會因為有些聚類沒有成員導致錯誤！如果你遇到了IndexError，估算中止，那麼在不丟失精確度的前提下減少聚類數吧。
我們明確了終止條件。如果當前迴圈與前一次迴圈（在成員關係函式上的變化）差異小於0.01，我們的模型就會終止迴圈。
這個模型返回一系列物件。centroids是五個聚類的座標。u存著每個觀測值的成員值；結構如下：

/*
The method findClusters_cmeans took 2.25 sec to run.
[[0.20158559 0.09751069 0.07842359 ... 0.16641284 0.16997486 0.18780077]
 [0.14041632 0.05272408 0.04176255 ... 0.25766411 0.13334598 0.13312636]
 [0.15019893 0.05824498 0.04623189 ... 0.14150986 0.2692789  0.14128872]
 [0.13702899 0.05074051 0.04020164 ... 0.28432326 0.27960719 0.38820815]
 [0.37077018 0.74077974 0.79338032 ... 0.15008993 0.14779308 0.14957599]]
Pseudo_F:  8340.883112129994
Davis-Bouldin:  1.3062853046748828
Silhouette score:  0.35108693925427953
 */

如果將每一行的所有列值加起來，你會發現（果然）一直是1。
剩下的返回物件我們不用太操心：u0是每個觀測值的初始成員種子，d是最終的歐幾里得距離矩陣，jm是目標函式的變更歷史，p是估算模型所用的迴圈數，fpc是模糊劃分係數（fuzzy partition coefficient）。
為了計算我們的指標，我們需要將聚類賦給每個觀測值。我們用NumPy的.argmax（...）方法來做這事：

 # assess the clusters effectiveness
labels = [
    np.argmax(elem) for elem in u.transpose()
]

這個方法返回列表中最大元素的索引。我們的u矩陣的初始分佈是u_cluster x n_sample，我們需要先用.transpose（）方法轉置u矩陣，以便按行迴圈，每一行都是我們樣本中的觀測值。

4.7
使用層次模型聚類資料

層次聚類模型的目標是構建聚類的分層。從概念上講，你可以理解成是一棵聚類的決策樹：基於聚類之間的相似度（或相異度），它們聚合（或切分）成更抽象（或更具體）的聚類。聚合的做法常被稱作自底向上，而切分的做法被稱作自頂向下。
需裝好pandas、SciPy和PyLab。
>pip install pylab

/* Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
......
Successfully built pyusb PyVISA
Installing collected packages: pyusb, PyVISA, PyVISA-py, pylab
Successfully installed PyVISA-1.10.1 PyVISA-py-0.3.1 pylab-0.0.2 pyusb-1.0.2
FINISHED */

聚合演算法的複雜度是O（n3），對於大規模的資料集來說，層次聚類可能會極慢。估算我們的模型時，我們使用的是單連結演算法，其演算法複雜度更好，是O（n2），不過資料集很大時依然會很慢（clustering_hierarchical.py檔案）：

1 import pylab as pl
2     def findClusters_link(data):
3     '''
4         Cluster data using single linkage hierarchical
5         clustering
6         使用單連結分層聚類資料
7     '''
8     # return the linkage object
9     return cl.linkage(data, method='single')

 

原理:程式碼簡單之極：只需呼叫SciPy的.linkage（...）方法，傳入資料，指定方法。
詳細說來，.linkage（...）方法是基於選擇的方法，根據特定的距離維度聚合（或切分）聚類；method='single'根據兩個聚類之間每個點的最小距離（所以是O（n2））聚合聚類。
關於所有指標的列表，參考SciPy文件：http://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html#scipy.cluster.hierarchy.linkage。
我們可以將聚合（或拆分）以樹的形式具象化：

 1 ####import pylab as pl
 2 import matplotlib.pylab as pl
 3 
 4 # cluster the data
 5 cluster = findClusters_ward(selected)
 6 
 7 # plot the clusters
 8 fig  = pl.figure(figsize=(16,9))
 9 ax   = fig.add_axes([0.1, 0.1, 0.8, 0.8])
10 dend = cl.dendrogram(cluster, truncate_mode='level', p=20)
11 ax.set_xticks([])
12 ax.set_yticks([])
13 
14 # save the figure
15 fig.savefig(
16     '../../Data/Chapter04/hierarchical_dendrogram.png',
17     dpi=300
18 )

首先，我們使用PyLab和.add_axes（...）建立一幅圖。傳入函式的列表是圖的座標：[x，y，width，height]，數值已調整到0和1之間，1意味著100%。.dendrogram（...）方法輸入所有的連結，聚類並建立圖形。我們不希望輸出所有的連結，所以我們在p＝20層時切掉了樹。我們也不需要座標標籤，所以我們使用.set_{x，y}ticks（...）方法關掉了它們。
對於我們的資料集，你的樹看起來類似這樣：
邀月工作室
你可以看出來，越接近頂端，你得到的聚類越不清晰（發生了更多的聚合）。
我們也可以使用MLPy估算分層聚類模型（clustering_hierarchical_alternative.py檔案）：

 1 def findClusters_ward(data):
 2     '''
 3         Cluster data using Ward's hierarchical clustering
 4     '''
 5     # create the classifier object
 6     ward = ml.MFastHCluster(
 7         method='ward'
 8     )
 9 
10     # fit the data
11     ward.linkage(data)
12 
13     return ward
14 
15 # the file name of the dataset
16 r_filename = '../../Data/Chapter04/bank_contacts.csv'
17 
18 # read the data
19 csv_read = pd.read_csv(r_filename)
20 
21 # select variables
22 selected = csv_read[['n_duration','n_nr_employed',
23         'prev_ctc_outcome_success','n_euribor3m',
24         'n_cons_conf_idx','n_age','month_oct',
25         'n_cons_price_idx','edu_university_degree','n_pdays',
26         'dow_mon','job_student','job_technician',
27         'job_housemaid','edu_basic_6y']]
28 
29 # cluster the data
30 cluster = findClusters_ward(selected)
31 
32 # assess the clusters effectiveness
33 labels = cluster.cut(20)
34 centroids = hlp.getCentroids(selected, labels)

我們使用.MFastHCluster（...）方法，傳入ward作為引數。分層聚類的沃德方法使用的是沃德方差最小化演算法。.linkage（...）方法可看作k均值.fit（...）方法的等價物；它將模型用到資料上，找出所有的連結。
估算完成後，我們使用.cut（...）方法在第20層剪掉這棵樹，以得到可管理的聚類數。返回的物件包括資料集中每個觀測值的聚類成員索引。
沃德聚類方法不返回聚類的中心（因為隨著我們切樹的層次不同，中心會變化），為了能使用評估方法，我們寫了一個方法，可返回指定標籤的中心：

 1 def getCentroids(data, labels):
 2     '''
 3         Method to get the centroids of clusters in clustering
 4         models that do not return the centroids explicitly
 5     '''
 6     # create a copy of the data
 7     data = data.copy()
 8 
 9     # apply labels
10     data['predicted'] = labels
11 
12     # and return the centroids
13     return np.array(data.groupby('predicted').agg('mean'))

首先，我們為資料建立一份本地的副本（因為我們不希望修改傳入方法的這個物件），並建立一個新列，叫作predicted。然後為predicted列去重後的每一層計算平均值，並投射為np.array（...）。
現在有了中心，我們可以評估模型的表現了。
如果你對分層聚類感興趣，我建議閱讀這篇文章：http://www.econ.upf.edu/~michael/stanford/maeb7.pdf。

4.8使用DBSCAN和BIRCH演算法發現潛在的訂閱者

DBSCAN（Density-based Spatial Clustering of Applications with Noise，基於密度的空間聚類）和BIRCH（Balanced Iterative Reducing and Clustering using Hierarchies，利用層次方法的平衡迭代規約和聚類）是首批能有效處理帶噪音的資料的演算法。這裡的噪音可以理解為，資料集中與其他部分相比，所放地方不恰當的點；DBSCAN將這種觀測值放到一個未歸類的桶，而BIRCH認為它們是離群值，直接挪走。
需裝好pandas和Scikit。
兩種演算法都能在Scikit中找到。要使用DBSCAN，使用clustering_dbscan.py檔案中的程式碼：

 1 import sklearn.cluster as cl
 2 import numpy as np
 3 
 4 @hlp.timeit
 5 def findClusters_DBSCAN(data):
 6     '''
 7         Cluster data using DBSCAN algorithm 使用DBScan演算法 聚類資料
 8     '''
 9     # create the classifier object
10     dbscan = cl.DBSCAN(eps=1.2, min_samples=200)
11 
12     # fit the data
13     return dbscan.fit(data)
14 
15 # cluster the data
16 cluster = findClusters_DBSCAN(selected)
17 
18 # assess the clusters effectiveness
19 labels = cluster.labels_ + 1
20 centroids = hlp.getCentroids(selected, labels)

邀月工作室

至於BIRCH模型，檢視clustering_birch.py檔案中的程式：

 1 import sklearn.cluster as cl
 2 import numpy as np
 3 
 4 @hlp.timeit
 5 def findClusters_Birch(data):
 6     '''
 7         Cluster data using BIRCH algorithm
 8     '''
 9     # create the classifier object
10     birch = cl.Birch(
11         branching_factor=100,
12         n_clusters=4,
13         compute_labels=True,
14         copy=True
15     )
16 
17     # fit the data
18     return birch.fit(data)
19 
20 
21 # cluster the data
22 cluster = findClusters_Birch(selected)
23 
24 # assess the clusters effectiveness
25 labels = cluster.labels_
26 centroids = hlp.getCentroids(selected, labels)
/* 
The method findClusters_Birch took 0.99 sec to run.
Pseudo_F:  1733.8811847551842
Davis-Bouldin:  1.2605820279777722
*/

原理:要估算DBSCAN，我們使用Scikit的.DBSCAN（...）方法。這個方法可以接受很多引數。eps控制的是認為兩個樣本鄰接的最大距離。min_samples引數控制的是鄰接關係；一旦超過這個閾值，至少有這麼多鄰居的點就被當成一個核心點。
如前所述，DBSCAN中的殘餘值（離群值）會歸到一個未歸類的桶，返回時標籤為-1。在可以評估聚類方法的表現前，我們需要移動標籤，使它們從0開始（0是未聚類的資料點）。這個方法不返回中心的座標，所以我們使用前一技巧中的.getCentroids（...）方法。
估算BIRCH也不費什麼勁（就我們這部分來說）。我們使用Scikit的.Birch（...）方法。branching_factor引數控制的是樹的父節點中最多有多少子聚類（或觀測值）；超過這個值時，聚類（以及子聚類）會遞迴拆分。n_clusters引數告訴方法我們想將資料聚成多少個類。將compute_label設定成True是讓方法在結束時為每個觀測值準備並儲存標籤。
BIRCH演算法直接丟棄離群值。如果你不希望方法修改原始資料，可以將copy引數設定為True：這會複製原始資料。
參考:
DBSCAN的論文：https://www.aaai.org/Papers/KDD/1996/KDD96-037.pdf。
BIRCH的論文：http://www.cs.sfu.ca/CourseCentral/459/han/papers/zhang96.pdf。
建議也閱讀.DBSCAN（...）和.Birch（...）方法的文件：

http://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html#sklearn.cluster.DBSACN；

http://scikit-learn.org/stable/modules/generated/sklearn.cluster.Birch.html#sklearn.cluster.Birch。

 

  第4章完。

python資料分析個人學習讀書筆記-目錄索引

 

隨書原始碼官方下載：
http://www.hzcourse.com/web/refbook/detail/7821/92

 

相關文章