項目和代碼下載地址：文章末尾。
最後結果：
高斯混合模型 EM 算法的 Python 實現
2017年11月24日 Wray Zheng704
文章目錄
- 高斯混合模型簡介
- GMM 的 EM 算法
- E 步
- M 步
- Python 實現
- 測試結果
- 測試數據
- 參考資料
- 相關文章
高斯混合模型簡介
首先簡單介紹一下，高斯混合模型（GMM, Gaussian Mixture Model）是多個高斯模型的線性疊加，高斯混合模型的概率分佈可以表示如下：
P(x)=∑
k=1
K
α
k
ϕ(x;μ
k
,Σ
k
)
P(x)=∑k=1Kαkϕ(x;μk,Σk)
其中，K 表示模型的個數，α
k
αk 是第 k 個模型的係數，表示出現該模型的概率，ϕ(x;μ
k
,Σ
k
)
ϕ(x;μk,Σk) 是第 k 個高斯模型的概率分佈。
這裡討論的是多個隨機變量的情況，即多元高斯分佈，因此高斯分佈中的參數不再是方差 σ
k
σk，而是協方差矩陣 Σ
k
Σk 。
我們的目標是給定一堆沒有標籤的樣本和模型個數 K，以此求得混合模型的參數，然後就可以用這個模型來對樣本進行聚類。
GMM 的 EM 算法
我們知道，EM 算法是通過不斷迭代來求得最佳參數的。在執行該算法之前，需要先給出一個初始化的模型參數。
我們讓每個模型的 μ
μ 為隨機值，Σ
Σ 為單位矩陣，α
α 為 1
K
1K，即每個模型初始時都是等概率出現的。
EM 算法可以分為 E 步和 M 步。
E 步
直觀理解就是我們已經知道了樣本 x
i
xi，那麼它是由哪個模型產生的呢？我們這裡求的就是：樣本 x
i
xi 來自於第 k 個模型的概率，我們把這個概率稱為模型 k 對樣本 x
i
xi 的“責任”，也叫“響應度”，記作 γ
ik
γik，計算公式如下：
γ
ik
=α
k
ϕ(x
i
;μ
k
,Σ
k
)
∑
K
k=1
α
i
ϕ(x
i
;μ
k
,Σ
k
)
γik=αkϕ(xi;μk,Σk)∑k=1Kαiϕ(xi;μk,Σk)
M 步
根據樣本和當前 γ
γ 矩陣重新估計參數，注意這裡 x 為列向量：
μ
k
=1
N
k
∑
i=1
N
γ
ik
x
i
μk=1Nk∑i=1Nγikxi
Σ
k
=1
N
k
∑
i=1
N
γ
ik
(x
i
−μ
k
)(x
i
−μ
k
)
T
Σk=1Nk∑i=1Nγik(xi−μk)(xi−μk)T
α
k
=N
k
N
αk=NkN
Python 實現
在給出代碼前，先作一些說明。
- 在對樣本應用高斯混合模型的 EM 算法前，需要先進行數據預處理，即把所有樣本值都縮放到 0 和 1 之間。
- 初始化模型參數時，要確保任意兩個模型之間參數沒有完全相同，否則迭代到最後，兩個模型的參數也將完全相同，相當於一個模型。
- 模型的個數必須大於 1。當 K 等於 1 時相當於將樣本聚成一類，沒有任何意義。
gmm.py 文件：
# -*- coding: utf-8 -*-# ----------------------------------------------------# Copyright (c) 2017, Wray Zheng. All Rights Reserved.# Distributed under the BSD License.# ----------------------------------------------------import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import multivariate_normalDEBUG = True####################################################### 調試輸出函數# 由全局變量 DEBUG 控制輸出######################################################def debug(*args, **kwargs): global DEBUG if DEBUG: print(*args, **kwargs)####################################################### 第 k 個模型的高斯分佈密度函數# 每 i 行表示第 i 個樣本在各模型中的出現概率# 返回一維列表######################################################def phi(Y, mu_k, cov_k): norm = multivariate_normal(mean=mu_k, cov=cov_k) return norm.pdf(Y)####################################################### E 步：計算每個模型對樣本的響應度# Y 為樣本矩陣，每個樣本一行，只有一個特徵時為列向量# mu 為均值多維數組，每行表示一個樣本各個特徵的均值# cov 為協方差矩陣的數組，alpha 為模型響應度數組######################################################def getExpectation(Y, mu, cov, alpha): # 樣本數 N = Y.shape[0] # 模型數 K = alpha.shape[0] # 為避免使用單個高斯模型或樣本，導致返回結果的類型不一致 # 因此要求樣本數和模型個數必須大於1 assert N > 1, "There must be more than one sample!" assert K > 1, "There must be more than one gaussian model!" # 響應度矩陣，行對應樣本，列對應響應度 gamma = np.mat(np.zeros((N, K))) # 計算各模型中所有樣本出現的概率，行對應樣本，列對應模型 prob = np.zeros((N, K)) for k in range(K): prob[:, k] = phi(Y, mu[k], cov[k]) prob = np.mat(prob) # 計算每個模型對每個樣本的響應度 for k in range(K): gamma[:, k] = alpha[k] * prob[:, k] for i in range(N): gamma[i, :] /= np.sum(gamma[i, :]) return gamma####################################################### M 步：迭代模型參數# Y 為樣本矩陣，gamma 為響應度矩陣######################################################def maximize(Y, gamma): # 樣本數和特徵數 N, D = Y.shape # 模型數 K = gamma.shape[1] #初始化參數值 mu = np.zeros((K, D)) cov = [] alpha = np.zeros(K) # 更新每個模型的參數 for k in range(K): # 第 k 個模型對所有樣本的響應度之和 Nk = np.sum(gamma[:, k]) # 更新 mu # 對每個特徵求均值 for d in range(D): mu[k, d] = np.sum(np.multiply(gamma[:, k], Y[:, d])) / Nk # 更新 cov cov_k = np.mat(np.zeros((D, D))) for i in range(N): cov_k += gamma[i, k] * (Y[i] - mu[k]).T * (Y[i] - mu[k]) / Nk cov.append(cov_k) # 更新 alpha alpha[k] = Nk / N cov = np.array(cov) return mu, cov, alpha####################################################### 數據預處理# 將所有數據都縮放到 0 和 1 之間######################################################def scale_data(Y): # 對每一維特徵分別進行縮放 for i in range(Y.shape[1]): max_ = Y[:, i].max() min_ = Y[:, i].min() Y[:, i] = (Y[:, i] - min_) / (max_ - min_) debug("Data scaled.") return Y####################################################### 初始化模型參數# shape 是表示樣本規模的二元組，(樣本數, 特徵數)# K 表示模型個數######################################################def init_params(shape, K): N, D = shape mu = np.random.rand(K, D) cov = np.array([np.eye(D)] * K) alpha = np.array([1.0 / K] * K) debug("Parameters initialized.") debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n") return mu, cov, alpha####################################################### 高斯混合模型 EM 算法# 給定樣本矩陣 Y，計算模型參數# K 為模型個數# times 為迭代次數######################################################def GMM_EM(Y, K, times): Y = scale_data(Y) mu, cov, alpha = init_params(Y.shape, K) for i in range(times): gamma = getExpectation(Y, mu, cov, alpha) mu, cov, alpha = maximize(Y, gamma) debug("{sep} Result {sep}".format(sep="-" * 20)) debug("mu:", mu, "cov:", cov, "alpha:", alpha, sep="\n") return mu, cov, alpha main.py 文件：
# -*- coding: utf-8 -*-# ----------------------------------------------------# Copyright (c) 2017, Wray Zheng. All Rights Reserved.# Distributed under the BSD License.# ----------------------------------------------------import matplotlib.pyplot as pltfrom gmm import *# 設置調試模式DEBUG = True# 載入數據Y = np.loadtxt("gmm.data")matY = np.matrix(Y, copy=True)# 模型個數，即聚類的類別個數K = 2# 計算 GMM 模型參數mu, cov, alpha = GMM_EM(matY, K, 100)# 根據 GMM 模型，對樣本數據進行聚類，一個模型對應一個類別N = Y.shape[0]# 求當前模型參數下，各模型對樣本的響應度矩陣gamma = getExpectation(matY, mu, cov, alpha)# 對每個樣本，求響應度最大的模型下標，作為其類別標識category = gamma.argmax(axis=1).flatten().tolist()[0]# 將每個樣本放入對應類別的列表中class1 = np.array([Y[i] for i in range(N) if category[i] == 0])class2 = np.array([Y[i] for i in range(N) if category[i] == 1])# 繪製聚類結果plt.plot(class1[:, 0], class1[:, 1], 'rs', label="class1")plt.plot(class2[:, 0], class2[:, 1], 'bo', label="class2")plt.legend(loc="best")plt.title("GMM Clustering By EM Algorithm")plt.show()測試結果
樣本 1 原始類別：
樣本 1 聚類結果：
樣本 2 原始類別：
樣本 2 聚類結果：
測試數據
或者通過下列代碼生成測試數據：
import numpy as npimport matplotlib.pyplot as pltcov1 = np.mat("0.3 0;0 0.1")cov2 = np.mat("0.2 0;0 0.3")mu1 = np.array([0, 1])mu2 = np.array([2, 1])sample = np.zeros((100, 2))sample[:30, :] = np.random.multivariate_normal(mean=mu1, cov=cov1, size=30)sample[30:, :] = np.random.multivariate_normal(mean=mu2, cov=cov2, size=70)np.savetxt("sample.data", sample)plt.plot(sample[:30, 0], sample[:30, 1], "bo")plt.plot(sample[30:, 0], sample[30:, 1], "rs")plt.show()完整項目下載地址：
https://github.com/wrayzheng/gmm-em-clustering.git
分享到:
閱讀更多 Python樂園 的文章
