1. EM算法-數學基礎
2. EM算法-原理詳解
3. EM算法-高斯混合模型GMM
4. EM算法-高斯混合模型GMM詳細代碼實現
5. EM算法-高斯混合模型GMM+Lasso
1. 前言
EM的前3篇博文分別從數學基礎、EM通用算法原理、EM的高斯混合模型的角度介紹了EM算法。按照慣例，本文要對EM算法進行更進一步的探究。就是動手去實踐她。
2. GMM實現
我的實現邏輯基本按照GMM算法流程中的方式實現。需要全部可運行代碼，請移步我的github。
輸入：觀測數據\(x_1,x_2,x_3,...,x_N\)
對輸入數據進行歸一化處理
#數據預處理
def scale_data(self):
for d in range(self.D):
max_ = self.X[:, d].max()
min_ = self.X[:, d].min()
self.X[:, d] = (self.X[:, d] - min_) / (max_ - min_)
self.xj_mean = np.mean(self.X, axis=0)
self.xj_s = np.sqrt(np.var(self.X, axis=0))
輸出：GMM的參數
- 初始化參數
#初始化參數
def init_params(self):
self.mu = np.random.rand(self.K, self.D)
self.cov = np.array([np.eye(self.D)] * self.K) * 0.1
self.alpha = np.array([1.0 / self.K] * self.K)
- E步：根據當前模型，計算模型\(k\)對\(x_i\)的影響
\[\gamma_{ik}=\frac{\pi_k\mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\sum_{k=1}^K\pi_k\mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \]
#e步，估計gamma
def e_step(self, data):
gamma_log_prob = np.mat(np.zeros((self.N, self.K)))
for k in range(self.K):
gamma_log_prob[:, k] = log_weight_prob(data, self.alpha[k], self.mu[k], self.cov[k])
log_prob_norm = logsumexp(gamma_log_prob, axis=1)
log_gamma = gamma_log_prob - log_prob_norm[:, np.newaxis]
return log_prob_norm, np.exp(log_gamma)
- M步：計算\(\mu_{k+1},\Sigma_{k+1}^2,\pi_{k+1}\)。
\[n_k=\sum_{i=1}^N\gamma_{ik} \]
\[\mu_{k+1}=\frac{1}{n_k}\sum_{i=1}^N\gamma_{ik}x_i \]
\[\Sigma_{k+1}^2=\frac{1}{n_k}\sum_{i=1}^N\gamma_{ik}(x_i-\mu_k)^2 \]
\[\pi_{k+1}=\frac{n_k}{N} \]
#m步，最大化loglikelihood
def m_step(self):
newmu = np.zeros([self.K, self.D])
newcov = []
newalpha = np.zeros(self.K)
for k in range(self.K):
Nk = np.sum(self.gamma[:, k])
newmu[k, :] = np.dot(self.gamma[:, k].T, self.X) / Nk
cov_k = self.compute_cov(k, Nk)
newcov.append(cov_k)
newalpha[k] = Nk / self.N
newcov = np.array(newcov)
return newmu, newcov, newalpha
- 重復2，3兩步直到收斂
最后加上loglikelihood的計算方法。
- 基本的計算方法按照公式定義。
\[L(\theta) = \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}Q_i(z^{(i)})P(x^{(i)},z^{(i)}|\theta)\;\;\;s.t.\sum\limits_{z}Q_i(z^{(i)}) =1 \]
實現如下
def loglikelihood(self):
P = np.zeros([self.N, self.K])
for k in range(self.K):
P[:,k] = prob(self.X, self.mu[k], self.cov[k])
return np.sum(np.log(P.dot(self.alpha)))
但是這樣的實現會有2個問題。
- 非矩陣運算，速度慢。
- 非常容易underflow，因為\(P.dot(self.alpha)\)非常容易是一個很小的數，系統把它當作0處理。
- 使用以下\(LogSumExp\)公式進行改進，並且令\(a_h = log(Q_i(z^{(i)}))+log(P(x^{(i)},z^{(i)}|\theta))\)，具體實現看github：
\[log(\sum_hexp(a_h)) = m + log(\sum_hexp(a_h - m))\;\;\;m=max(a_h) \]
3. 總結
首先gmm算法會很容易出現underflow和overflow，所以處理的時候有點麻煩。但是\(LogSumExp\)能解決大部分這個問題。還有就是我的實現方式是需要協方差矩陣一定要是正定矩陣，所以我的代碼中也做了處理。我們好想還不能夠滿足於最基礎的GMM算法，所以在下一篇文章中我們要對GMM加入一個懲罰項，並且用對角矩陣的方式代替協方差矩陣。