Kmeans算法

1、概述

該方法屬於無監督學習算法（無y值）。根據已有的數據，利用距離遠近的思想將目標數據集聚為指定的k個簇。簇內樣本越相似，聚類的效果越好。需要注意的是如若數據存在量綱上的差異，必須先進行標簽化處理。或者數據集中含有離散型字符變量，需先設置成啞變量或進行數值化。對於未知簇個數的數據集，需要先確定簇數，再進行聚類操作。

應用方面常用在聚類模型中。比如電商平台歷史交易數據，划分不同的價值等級（如VIP，高價值，潛在價值，低價值）等等，還可用於異常值的監控（遠離任何簇的異常點），如（信用卡交易異常，社交平台點擊異常（如遇爬蟲腳本），電商網站交易異常（一張銀行卡被用於多個Id支付，送貨地址大致集中在相近區域））

2、過程

• 任選k個對象作為初始簇的中心
• 計算各對象到初始簇的中心的距離，重新划分簇
• 計算各個簇的均值，為新的簇中心
• 再次計算各對象到新的簇中心的距離，重新划分簇，然后再重新計算均值，划分簇，repeat...until 新簇（簇中心）不在發生變化。

3、缺點

• 結果容易受異常值（離群點）影響。（中心點通過樣本均值確定）
• 該方法不適合樣本中非球形的簇。

4、實現Kmeans聚類方式

• 原始代碼的方式
• sklean

模擬kmeans實現案例：

數據集：

 

 　　

 

 

手動模擬 實現過程：

 

 

代碼實現：

import numpy as np 

def kmeans(X, k, maxIt):   # X數據集矩陣,maxIt為最多迭代的次數
	numPoints, numDim = X.shape

	dataSet = np.zeros((numPoints, numDim + 1))
	dataSet[:,:-1] = X

	# centroids = dataSet[np.random.randint(numPoints, size = k), :]   # 隨機選2行作為初始簇中心點
	centroids = dataSet[0:2,:]  # 選擇1、2行作為初始簇中心
	centroids[:,-1] = range(1,k+1)  #划分簇類別

	iterations = 0
	OldCentriods = None

	while not shouldStop(OldCentriods, centroids, iterations, maxIt):
		print("iteration:\n", iterations)
		print("dataSet:\n", dataSet)
		print("centroids:\n", centroids)

		OldCentriods = np.copy(centroids)
		iterations +=1

		updateLabels(dataSet, centroids)  # 對中心點和數據集重新歸類
		centroids = getCentroids(dataSet, k)  # 根據均值生成新中心點

	return dataSet


def shouldStop(OldCentriods, centroids , iterations , maxIt):
	if iterations > maxIt:
		return True 
	return np.array_equal(OldCentriods,centroids)  

def updateLabels(dataSet, centroids):
	numPoints , numDim = dataSet.shape
	for i in range(0, numPoints):
		dataSet[i, -1] = getLabelFromClosestCentroid(dataSet[i, :-1], centroids)

# 遍歷數據中的每個對象，每個對象跟2個中心點距離比較以后，獲取新的label將數據集中的對象label更新
def getLabelFromClosestCentroid(dataSetRow, centroids):
	label = centroids[0,-1]
	minDist = np.linalg.norm(dataSetRow - centroids[0, :-1])

	for i in range(1, centroids.shape[0]):
		dist = np.linalg.norm(dataSetRow - centroids[i, :-1])
		if dist < minDist:
			minDist = dist
			label = centroids[i,-1]

	print('minDist', minDist)
	return label 

# 求均值中心點
def getCentroids(dataSet, k):
	result = np.zeros((k, dataSet.shape[1]))

	for i in range(1,1+k):
		oneCluster = dataSet[dataSet[:,-1] ==i , :-1] 
		result[i-1,:-1]  = np.mean(oneCluster,axis = 0)
		result[i-1,-1]  = i

	return result


x1 = np.array([1,1])
x2 = np.array([2,1])
x3 = np.array([4,3])
x4 = np.array([5,4])

testX = np.vstack((x1,x2,x3,x4))
result = kmeans(testX,2, 10)

print('final result:' ,result)

　　結果迭代2次，打印結果與以上一致。

 

 若改變數據集為10個點

A1 = np.array([0.5,2])
A2 = np.array([0.8,3])
A3 = np.array([1.2,0.6])
A4 = np.array([1.6,2.2])
A5 = np.array([2.4,3.6])
A6 = np.array([2.5,2.8])
B1 = np.array([2.2,1.8])
B2 = np.array([3,2.5])
B3 = np.array([2.8,1.6])
B4 = np.array([4,1])

　　以上代碼會生成不同的聚類效果：

 

 取k=3時結果：

5、k值的確定方法

常采用探索法確定k值：給定不同的k值，對比某些評估指標的變動情況。進而選擇合理的k值。常用方法如下

• 簇內離差平方和拐點法
• 輪廓系數法
• 間隔統計量法

1）拐點法

 使用sclearn子模塊cluster中KMeans類，思路是使簇內離差平方和總和最小。（隨着簇的增加，離差平方和之和趨於穩定（波動小於閥值），可以通過繪制各簇離差平方和之和與簇的個數之間的折線圖，觀察k值拐點變化的位置（斜率變化的點）。

離差反映的是變量與數學期望偏離的程度，集x-E(x)，其平方和為：

 

實現函數：

def k_SSE(X, clusters):
    K = range(1, clusters+1) # 設定簇k的取值范圍
    TSSE = []   # 存儲總的簇內離差平方和
    for k in K :
        SSE = []  # 存儲各個簇內離差平方和
        kmeans = KMeans(n_clusters=k)
        kmeans.fit(X)
        labels = kmeans.labels_  # 返回簇標簽
        centers = kmeans.cluster_centers_   # 返回簇中心
        for label in set(labels):
            SSE.append(np.sum((X.loc[labels == label,]-centers[label,:])**2))  # 計算個簇內離差平方和並存於列表SSE
        TSSE.append(np.sum(SSE)) # 計算k值下總離差平方和
        
    plt.style.use('ggplot')
    plt.plot(K,TSSE,'b*-')
    plt.xlabel('簇的個數')
    plt.ylabel('簇內離差平方和之和')
    plt.show()

X = pd.DataFrame(np.concatenate([ np.array([x1,y1]),np.array([x2,y2]),np.array([x3,y3])],axis = 1).T)
k_SSE(X,15)

 

用此方法驗證一下上面10個數據集的最佳k值如下：

 

從圖中可以發現k為2時圖線斜率明顯變化，且變化最大，所以取最佳k值為2。

2）輪廓系數法

 使用的是sklearn子模塊metrics中的silhouette_score函數，該方法要求簇的個數大於等於2（k>=2），輪廓系數計算公式如下：

 

 s(i)為輪廓系數，a(i)為樣本i與其簇內其他樣本距離的平均值。反應了簇內的密集性。b(i)為樣本i與其他每個簇簇內距離的平均值的最小值。

 

如圖表示的i與C1中其他樣本點距離的平均值為a(i)，i分別與C2、C3、C4各自簇內樣本點距離的平均值m1,m2,m3，比較3個值取最小值mx為b(i)。

s(i)的大小：s(i)值越接近-1，說明i分配不合理，聚類效果不佳，輪廓系數接近0說明i落在了模糊地帶，可能在簇的邊界處。輪廓系數接近1說明i的分配是合理的，聚類效果不錯。

實現函數：

from sklearn import metrics

def k_silhouette(X , clusters):
    K = range(2, clusters+1)
    S = []  # 空列表用於存放不同簇下的輪廓系數
    for k in K:
        kmeans = KMeans(n_clusters = k)
        kmeans.fit(X)
        labels = kmeans.labels_
        S.append(metrics.silhouette_score(X,labels, metric = 'euclidean')) #調用子模塊metrics中函數計算輪廓系數
    
    plt.style.use('ggplot')
    plt.plot(K, S, 'b*-')
    plt.xlabel('簇的個數')
    plt.ylabel('輪廓系數')
    plt.show()

k_silhouette(X,9)

　圖形：

 

 還是以上10個數據集，可見k=2時，輪廓系數大於0，且最接近1，所以k=2時分類效果最好。

3）間隔統計量法 （Gap Statistic）

公式：

 

 

 

各參數Dk：簇內樣本點之間的歐式距離

nk為第k個簇內的樣本量。uk為第k個簇內的樣本均值，Wk為Dk的標准化結果。為參照組數據集的Wk向量。

  的值為，通過參照數據集的數學期望與實際數據集的對數的差比較，找到實際數據集下降最快的k值。

 

 首次出現Gap(k)≥Gap(k+1)-sk+1時的k值即為最佳k值。

 

 

 

 

 案例：

import pandas as pd
import numpy as np  
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn import metrics
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus']=False  

# 隨機生成三組二元正態分布隨機數 
np.random.seed(1234)
mean1 = [0.5, 0.5]
cov1 = [[0.3, 0], [0, 0.3]]
x1, y1 = np.random.multivariate_normal(mean1, cov1, 1000).T

mean2 = [0, 8]
cov2 = [[1.5, 0], [0, 1]]
x2, y2 = np.random.multivariate_normal(mean2, cov2, 1000).T

mean3 = [8, 4]
cov3 = [[1.5, 0], [0, 1]]
x3, y3 = np.random.multivariate_normal(mean3, cov3, 1000).T

plt.scatter(x1,y1)
plt.scatter(x2,y2)
plt.scatter(x3,y3)
plt.show()


# 構建數據集
X = pd.DataFrame(np.concatenate([np.array([x1,y1]),np.array([x2,y2]),np.array([x3,y3])], axis = 1).T)


# 自定義函數，計算簇內任意兩樣本之間的歐氏距離
def short_pair_wise_D(each_cluster):
    mu = each_cluster.mean(axis = 0)
    Dk = sum(sum((each_cluster - mu)**2)) * 2.0 * each_cluster.shape[0]
    return Dk

# 計算簇內的Wk值
def compute_Wk(data, classfication_result):
    Wk = 0
    label_set = set(classfication_result)   # 過濾所有數據集標簽值中的相同項，比如2個簇過濾后為{0,1}，10個簇為{0,1,...9}
    for label in label_set:
        each_cluster = data[classfication_result == label, :]   # 數據集中等於標簽值如1的那些值
        Wk = Wk + short_pair_wise_D(each_cluster)/(2.0*each_cluster.shape[0])   # Dk標准化
    return Wk   

# 計算GAP統計量 
def gap_statistic(X, B=10, K=range(1,11), N_init = 10):

    X = np.array(X)      # 將輸入數據集轉換為數組

    shape = X.shape
    tops = X.max(axis=0)  # 列方向求最大值
    bots = X.min(axis=0)
    dists = np.matrix(np.diag(tops-bots))  # 以每列最大值減最小值的值作為對角線元素其他元素為0組成矩陣
    rands = np.random.random_sample(size=(B,shape[0],shape[1]))  # 生成3維隨機數數組
    
    # 建立參照矩陣
    for i in range(B):
        rands[i,:,:] = rands[i,:,:]*dists+bots   # 對隨機矩陣機右乘對角矩陣（對應列乘以對角矩陣對角線上對應值），再加最小值
    
    gaps = np.zeros(len(K))
    Wks = np.zeros(len(K))
    Wkbs = np.zeros((len(K),B))

    for idxk, k in enumerate(K):        # 循環不同的k值，1~k個簇求Kmeans
        k_means =  KMeans(n_clusters=k)
        k_means.fit(X)
        classfication_result = k_means.labels_    # 對應於不同簇值（k）的每個對象的標簽值集合

        Wks[idxk] = compute_Wk(X, classfication_result)    # 存儲wk值


        for i in range(B):           # 計算每一個參照數據集下的各簇Wk值
            Xb = rands[i,:,:]
            k_means.fit(Xb)
            classfication_result_b = k_means.labels_
            Wkbs[idxk,i] = compute_Wk(Xb,classfication_result_b)  

    gaps = (np.log(Wkbs)).mean(axis = 1) - np.log(Wks)      
    sd_ks = np.std(np.log(Wkbs), axis=1)  
    sk = sd_ks*np.sqrt(1+1.0/B)    # sk為參照組數據集的無偏標准差

    gapDiff = gaps[:-1] - gaps[1:] + sk[1:]    # 用於判別最佳k的標准，當gapDiff首次為正時，對應的k即為目標值
    
    plt.style.use('ggplot')
    plt.bar(np.arange(len(gapDiff))+1, gapDiff, color = 'steelblue')
    plt.xlabel('簇的個數')
    plt.ylabel('k的選擇標准')
    plt.show()
    
gap_statistic(X)    

　　圖形效果如下：

1）散點圖

 

 

2）k值-gapdiff趨勢圖

 

 

結果可以看出k=3時，條件Gap(k)-Gap(k+1)-sk+1的值首次大於0，固最佳k值取3。注意樣本數據量太小時，此方法結果可能存在一定的問題。可以其他2種方法比較確定最佳k值。