EM 算法-對鳶尾花數據進行聚類

公號：碼農充電站pro
主頁：https://codeshellme.github.io

以前介紹過K 均值算法，它是一種聚類算法。今天介紹EM 算法，它也是聚類算法，但比K 均值算法更加靈活強大。

EM 的全稱爲 Expectation Maximization，中文爲指望最大化算法，它是一個不斷觀察和調整的過程。

1，和麪過程

咱們先來看一下和麪的過程。

一般狀況下，若是你事先不知道面與水的比例，和麪過程多是下面這樣：

1. 先放入一些面和水。
2. 將麪糰揉拌均勻。
3. 觀察麪糰的稀稠程度：若是麪糰比較稀，則加入少量面；若是麪糰比較稠，則加入少量水。
4. 如此往復第2，3步驟，直到麪糰的稀稠程度達到預期。

這個和麪過程，就是一個EM 過程：

• 先加入一些面和水，將麪糰揉拌均勻，並觀察麪糰的稀稠程度。這是E 過程。
• 不斷的加入水和麪（調整水和麪的比例），直到達到預期麪糰程度。這是M 過程。

2，再看K 均值算法

在介紹K 均值 聚類算法時，展現過一個給二維座標點進行聚類的例子。

咱們再來看一下這個例子，以下圖：

上圖是一個聚類的過程，共有6 個步驟：

1. 初始時散點（綠色點）的分佈。
2. 隨機選出兩個中心點的位置，紅色x 和藍色x。
3. 計算全部散點分別到紅色x 和藍色x的距離，距離紅色x 近的點標紅色，距離藍色x近的點標藍色。
4. 從新計算兩個中心點的位置，兩個中心點分別移動到新的位置。
5. 從新計算全部散點分別到紅色x 和藍色x的距離，距離紅色x 近的點標紅色，距離藍色x近的點標藍色。
6. 再次計算兩個中心點的位置，兩個中心點分別移動到新的位置。中心點的位置幾乎再也不變化，聚類結束。

通過以上步驟就完成了一個聚類過程。

實際上，K 均值算法也是一個EM 過程：

• 肯定當前各種中心點的位置，並計算各個散點到現有的中心點的距離。這是E 過程。
• 將各個散點歸屬到各個類中，從新計算各個類的中心點，直到各種的中心點再也不改變。這是M 過程。

3，EM 算法

將二維數據點的聚類過程，擴展爲通常性的聚類問題，EM 算法是這樣一個模型：對於待分類的數據點，EM 算法讓計算機經過一個不斷迭代的過程，來構建一個分類模型。

EM 算法分爲兩個過程：

• E 過程：根據現有的模型，計算各個數據輸入到模型中的計算結果，這稱爲指望值計算過程，即 E 過程。
• M 過程：從新計算模型參數，以最大化指望值，這稱爲最大化過程，即M 過程。

以二維數據點的聚類過程爲例，咱們定義：

• 同一類中各個點到該類中心的平均距離爲 d；
• 不一樣類之間的平均距離爲 D。

那麼二維數據點聚類的M 過程，就是尋求最大化的D 和 -d。咱們但願的聚類結果是，同一類的點距離較近，不一樣類之間距離較遠。

EM 算法不是單個算法，而是一類算法。只要知足EM 這兩個過程的算法均可以被稱爲EM 算法。常見的EM 算法有GMM 高斯混合模型和HMM 隱馬爾科夫模型。

4，最大似然估計

高等數學中有一門課叫作《機率論與數理統計》，其中講到了參數估計。

統計推斷是數理統計的重要組成部分，它是指利用來自整體的樣本提供的信息，對整體的某些特徵進行估計或推斷，從而認識總體。

統計推斷分爲兩大類：參數估計和假設檢驗。

咱們假設，對於某個數據集，其分佈函數的基本形式已知，但其中含有一個或多個未知參數。

參數估計就是討論如何根據來自整體的樣本提供的信息對未知參數作出估計。參數估計包括點估計和區間估計。其中，點估計中有兩種方法：矩估計法和最大似然估計法。

最大似然估計是一種經過已知結果，估計未知參數的方法。

5，EM 算法原理

EM 算法使用的是最大似然估計的原理，它經過觀察樣本，來找出樣本的模型參數。

下面經過一個投硬幣的例子，來看下EM 算法的計算過程。

這個例子來自《Nature》（天然）期刊的論文《What is the expectation maximization algorithm？》（什麼是指望最大化算法？）。

假定有兩枚不一樣的硬幣 A 和 B，它們的重量分佈 θ_A​ 和 θ_B​ 是未知的，則能夠經過拋擲硬幣，計算正反面各自出現的次數來估計θ_A​ 和 θ_B​。

方法是在每一輪中隨機抽出一枚硬幣拋擲 10 次，一樣的過程執行 5 輪，根據這 50 次投幣的結果來計算 θ_A​ 和 θ_B​​ 的最大似然估計。

投擲硬幣的過程，記錄以下：

第1 到5 次分別投擲的硬幣是 B，A，A，B，A。H 表明正面，T 表明負面。將上圖轉化爲表格，以下：

次數	硬幣	正面數	負面數
1	B	5	5
2	A	9	1
3	A	8	2
4	B	4	6
5	A	7	3

經過這個表格，能夠直接計算 θ_A​ 和 θ_B​​，以下：

顯然，若是知道每次投擲的硬幣是A 仍是B，那麼計算θ_A​ 和 θ_B 是很是簡單的。

可是，若是不知道每次投擲的硬幣是A 仍是B，該如何計算θ_A​ 和 θ_B 呢？

此時咱們將上面表格中的硬幣一列隱藏起來，這時硬幣就是隱變量。因此咱們只知道以下數據：

次數	正面數	負面數
1	5	5
2	9	1
3	8	2
4	4	6
5	7	3

這時想要計算 θ_A​ 和 θ_B，就要用最大似然估計的原理。

計算過程以下圖：

第一步

先爲 θ_A​ 和 θ_B 設定一個初始值，好比 θ_A = 0.6，θ_B = 0.5。

第二步

咱們知道每一輪投幣的正 / 負面的次數：

• 第1輪：5 正 5 負，計算出現這種結果的機率：
  • 若是是A 硬幣，那麼P(H₅T₅|A) = 0.6^5 * 0.4^5
  • 若是是B 硬幣，那麼P(H₅T₅|B) = 0.5^5 * 0.5^5
  • 將 P(H₅T₅|A) 和 P(H₅T₅|B) 歸一化處理，可得：
  • P(H₅T₅|A) = 0.45，P(H₅T₅|B) = 0.55
• 第2輪：9 正 1 負，計算出現這種結果的機率：
  • 若是是A 硬幣，那麼P(H₉T₁|A) = 0.6^9 * 0.4^1
  • 若是是B 硬幣，那麼P(H₉T₁|B) = 0.5^9 * 0.5^1
  • 將 P(H₉T₁|A) 和 P(H₉T₁|B) 歸一化處理，可得：
  • P(H₉T₁|A) = 0.8，P(H₉T₁|B) = 0.2
• 第3輪：8 正 2 負，計算出現這種結果的機率：
  • 若是是A 硬幣，那麼P(H₈T₂|A) = 0.6^8 * 0.4^2
  • 若是是B 硬幣，那麼P(H₈T₂|B) = 0.5^8 * 0.5^2
  • 將 P(H₈T₂|A) 和 P(H₈T₂|B) 歸一化處理，可得：
  • P(H₈T₂|A) = 0.73，P(H₈T₂|B) = 0.27
• 第4輪：4 正 6 負，計算出現這種結果的機率：
  • 若是是A 硬幣，那麼P(H₄T₆|A) = 0.6^4 * 0.4^6
  • 若是是B 硬幣，那麼P(H₄T₆|B) = 0.5^4 * 0.5^6
  • 將 P(H₄T₆|A) 和 P(H₄T₆|B) 歸一化處理，可得：
  • P(H₄T₆|A) = 0.35，P(H₄T₆|B) = 0.65
• 第5輪：7 正 3 負，計算出現這種結果的機率：
  • 若是是A 硬幣，那麼P(H₇T₃|A) = 0.6^7 * 0.4^3
  • 若是是B 硬幣，那麼P(H₇T₃|B) = 0.5^7 * 0.5^3
  • 將 P(H₇T₃|A) 和 P(H₇T₃|B) 歸一化處理，可得：
  • P(H₇T₃|A) = 0.65，P(H₇T₃|B) = 0.35

而後，根據每一輪的 P(H_mT_n|A) 和 P(H_mT_n|B)，能夠計算出每一輪的正 / 負面次數。

m 爲正面次數，n 爲負面次數。

對於硬幣A，結果以下：

輪數	P(HmTn|A)	m	n	正面數	負面數
1	0.45	5	5	0.45*5=2.2	0.45*5=2.2
2	0.8	9	1	0.8*9=7.2	0.8*1=0.8
3	0.73	8	2	0.73*8=5.9	0.73*2=1.5
4	0.35	4	6	0.35*4=1.4	0.35*6=2.1
5	0.65	7	3	0.65*7=4.5	0.65*3=1.9
總計	-	-	-	21.3	8.6

對於硬幣B，結果以下：

輪數	P(HmTn|B)	m	n	正面數	負面數
1	0.55	5	5	0.55*5=2.8	0.55*5=2.8
2	0.2	9	1	0.2*9=1.8	0.2*1=0.2
3	0.27	8	2	0.27*8=2.1	0.27*2=0.5
4	0.65	4	6	0.65*4=2.6	0.65*6=3.9
5	0.35	7	3	0.35*7=2.5	0.35*3=1.1
總計	-	-	-	11.7	8.4

第三步

根據上面兩個表格，能夠得出（第1次迭代的結果） θ_A 和 θ_B：

根據這個估計值，再次回到第一步去計算。

如此往復第1、2、三步，通過10次迭代以後，θ_A 和 θ_B 的估計值爲：

最終，θ_A 和 θ_B 將收斂到一個幾乎不變的值，此時迭代結束。這樣咱們就求解出了θ_A 和 θ_B 的最大似然估計值。

咱們將上述過程當中，第一步稱爲初始化參數，第二步稱爲觀察預期，第三步稱爲從新估計參數。

第1、二步爲E 過程，第三步爲M 過程，這就是EM 算法的過程。

若是咱們有一個待聚類的數據集，咱們把潛在的類別當作隱變量，樣本當作觀察值，這樣就能夠把聚類問題轉化成參數估計問題。這就是EM 聚類的原理。

6，硬聚類與軟聚類

與 K 均值算法相比，K 均值算法是經過距離來區分樣本之間的差異，且每一個樣本在計算的時候只能屬於一個分類，咱們稱之爲硬聚類算法。

而 EM 聚類在求解的過程當中，實際上每一個樣本都有必定的機率和每一個聚類相關，這叫作軟聚類算法。

7，EM 聚類的缺點

EM 聚類算法存在兩個比較明顯的問題。

第一個問題是，EM 算法計算複雜，收斂較慢，不太適合大規模數據集和高維數據。

第二個問題是，EM 算法不必定能給出全局最優解：

• 當優化的目標函數是凸函數時，必定能夠獲得全局最優解。
• 當優化的目標函數不是凸函數時，可能會獲得局部最優解，而非全局最優解。

8，GMM 高斯混合模型

上文中介紹過，常見的EM 算法有GMM 高斯混合模型和HMM 隱馬爾科夫模型。這裏主要介紹GMM 高斯混合模型的實現。

sklearn 庫的mixture 模塊中的GaussianMixture 類是GMM 算法的實現。

先來看下 GaussianMixture 類的原型：

GaussianMixture(
  n_components=1,
  covariance_type='full', 
  tol=0.001, 
  reg_covar=1e-06, 
  max_iter=100, 
  n_init=1, 
  init_params='kmeans', 
  weights_init=None, 
  means_init=None, 
  precisions_init=None, 
  random_state=None, 
  warm_start=False, 
  verbose=0, 
  verbose_interval=10)

這裏介紹幾個重要的參數：

• n_components：表明高斯混合模型的個數，也就是咱們要聚類的個數，默認值爲 1。
• covariance_type：表明協方差類型。一個高斯混合模型的分佈是由均值向量和協方差矩陣決定的，因此協方差的類型也表明了不一樣的高斯混合模型的特徵。協方差類型有 4 種取值：
  • full，表明徹底協方差，也就是元素都不爲 0；
  • tied，表明相同的徹底協方差；
  • diag，表明對角協方差，也就是對角不爲 0，其他爲 0；
  • spherical，表明球面協方差，非對角爲 0，對角徹底相同，呈現球面的特性。
• max_iter：表明最大迭代次數，默認值爲 100。

9，對鳶尾花數據集聚類

在《決策樹算法-實戰篇-鳶尾花及波士頓房價預測》中咱們介紹過鳶尾花數據集。這裏咱們使用GMM 算法對該數據進行聚類。

首先加載數據集：

from sklearn.datasets import load_iris

iris = load_iris()   	# 加載數據集
features = iris.data	# 獲取特徵集
labels = iris.target    # 獲取目標集

在聚類算法中，只須要特徵數據 features，而不須要目標數據labels，但可使用 labels 對聚類的結果作驗證。

構造GMM聚類：

from sklearn.mixture import GaussianMixture

# 原數據中有 3 個分類，因此這裏咱們將 n_components 設置爲 3
gmm = GaussianMixture(n_components=3, covariance_type='full')

對數據集進行聚類：

prediction_labels = gmm.fit_predict(features)

查看原始分類：

>>> print(labels)
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2]

查看聚類結果：

>>> print(prediction_labels)
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 2 1
 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
 2 2]

對比原始分類和聚類結果，聚類結果中只有個別數據分類錯誤，我用紅圈標了出來：

10，評估聚類結果

咱們可使用 Calinski-Harabaz 指標對聚類結果進行評估。

sklearn 庫實現了該指標的計算，即 calinski_harabasz_score 方法，該方法會計算出一個分值，分數越高，表明聚類效果越好，也就是相同類中的差別性小，不一樣類之間的差別性大。

下面對鳶尾花數據集的聚類結果進行評估，傳入特徵數據和聚類結果：

>>> from sklearn.metrics import calinski_harabasz_score
>>> calinski_harabasz_score(features, prediction_labels)
481.78070899745234

咱們也能夠傳入特徵數據和原始結果：

>>> calinski_harabasz_score(features, labels)
487.33087637489984

能夠看到，對於原始結果計算出的分值是487.33，對於預測結果計算出的分值是481.78，相差並很少，說明預測結果仍是不錯。

通常狀況下，一個須要聚類的數據集並無目標數據，因此只能對預測結果進行評分。咱們須要人工對聚類的含義結果進行分析。

11，總結

本篇文章主要介紹了以下內容：

• EM 算法的過程及原理，介紹了一個投擲硬幣的例子。
• 硬聚類與軟聚類的區別。
• EM 聚類的缺點：
  • 計算複雜度較大。
  • 有可能得不到全局最優解。
• 使用GMM 算法對鳶尾花數據進行聚類。
• 使用 Calinski-Harabaz 指標對聚類結果進行評估。

（本節完。）

---

推薦閱讀：

K 均值算法-如何讓數據自動分組

Apriori 算法-如何進行關聯規則挖掘

PageRank 算法-Google 如何給網頁排名

數據變換-歸一化與標準化

如何使用Python 進行數據可視化

---

歡迎關注做者公衆號，獲取更多技術乾貨。