EM（Expectation-Maximization）演算法是一種迭代式方法，主要應用於包含隱藏變數（latent variable）的參數估計，在無監督學習中有著廣泛的應用，EM演算法迭代包含兩步：

  • 利用估計的參數值來求對數似然期望（expectation）。
  • 通過最大化對數似然期望來更新參數。

上述是EM演算法的兩個基本迭代部分，在實際應用中EM演算法更多的看作是一種演算法思想，而不是特定的演算法步驟，所以接下來我將通過一些具體應用來進一步闡述EM演算法的主要思想。

1、K-Means

K-Means演算法是一種簡單好用的聚類方法，之所以把它提到最前面來是因為我覺得同樣作為無監督學習方法，K-Means在很多方面都深刻的體現著EM演算法的思想。

那麼首先我們來假設一個問題的形式：

input: 一組不帶label的數據： left{ x^{(1)},x^{(2)},...,x^{(m)} ight} ，現在要對其進行聚類。

output: 聚類模型

因為與監督學習不同，數據沒有帶label，所以對於之前介紹的邏輯回歸等分類方法就無法應用，而是需要我們從數據自身進行發掘。下面給出K-Means演算法的主要步驟：

1、隨機初始化聚類重心（cluster centroids）: mu _ { 1 } , mu _ { 2 } , ldots , mu _ { k } in mathbb { R } ^ { n }

2、Repeat until convergence：{ for every i, set: c ^ { ( i ) } : = arg min _ { j } left| x ^ { ( i ) } - mu _ { j } ight| ^ { 2 } for every j , set: mu _ { j } : = frac { sum _ { i = 1 } ^ { m } 1 left{ c ^ { ( i ) } = j ight} x ^ { ( i ) } } { sum _ { i = 1 } ^ { m } 1 left{ c ^ { ( i ) } = j ight} } }

可以看到，K-Means主要包含兩個迭代部分，首先是最小化重心與數據距離來更新參數，接著再對屬於同一類的數據求平均值來更新重心。當進行了一段時間後重心不發生變化，K-Means演算法也達到了收斂。

對應的，第一步更新參數類似於EM演算法中的maximization，而更新中心則對應於expectation步。雖然和真正的EM演算法有著很大的差別，但是這種迭代更新的思想卻是不謀而合，所以我們可以將其看作一種簡化版的EM來輔助理解。

2、高斯混合模型（GMM）

接下來要介紹的高斯混合模型（Gaussian Mixture Model）是EM演算法的一個重要應用，它應用廣泛，而很多時候EM演算法是對GMM進行參數估計的有效手段。

2.1單高斯分布

高斯分布又名正態分布，可以說是我們最常見的一種分布之一了，而且實驗表明許多受多個獨立因素影響的數據（如體重身高等）都滿足高斯分布。

如上圖所示就是一個典型的二維高斯分布

高斯分布有兩個重要參數：均值 mu 以及方差 sigma^{2} ，若隨機變數 X 服從高斯分布，則將其寫作： X sim N left( mu , sigma ^ { 2 } ight) 。

一維高斯分布的概率密度函數為： f ( x ) = frac { 1 } { sqrt { 2 pi } sigma } exp left( - frac { ( x - mu ) ^ { 2 } } { 2 sigma ^ { 2 } } ight)

多維高斯分布的概率密度函數為： f _ { mathbf { x } } left( x _ { 1 } , ldots , x _ { k } ight) = frac { 1 } { sqrt { ( 2 pi ) ^ { k } | mathbf { Sigma } | } } exp left( - frac { 1 } { 2 } ( mathbf { x } - oldsymbol { mu } ) ^ { mathrm { T } } mathbf { Sigma } ^ { - 1 } ( mathbf { x } - oldsymbol { mu } ) ight)

其中 | Sigma | 是協方差矩陣的行列式。

2.2高斯混合模型

顧名思義，高斯混合模型就是多個高斯分布線性組合在一起形成的新的分布，理論上來說，任何分布都可以由多個混合分布組合表示，但由於高斯混合模型地應用廣泛，這裡我們主要討論它。

如下圖所示就是兩個二維高斯分布的混合：

那麼如何來描述上述的分布呢？很顯然使用單個高斯模型是不夠準確的（從圖中可以直觀地看到兩個不同地聚類），所以由此引入了高斯混合模型地數學表達：

\p ( oldsymbol { x } ) = sum _ { k = 1 } ^ { K } alpha _ { k } mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } ) qquad(1)

其中 alpha_{k} 是權重，且sum _ { k = 1 } ^ { K } alpha _ { k }=1 ，而 mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , mathbf { Sigma } _ { k } ) 代表第k個分量模型。

如何理解GMM呢，首先我們先引入一個K維的二值變數 z ，其中某個特定的元素 z_{k}=1 ，其餘元素都為0，因此 z_{k} 滿足 z _ { k } in { 0,1 } 且 sum _ { k } z _ { k } = 1 。

根據上述定義， z 有K個不同的狀態，用來表示觀測數據 x_{i} 屬於哪個分量模型，因此每一個觀測數據 x_{i} 都對應一個 z^{i} ，譬如 z_{k}^{i}=1 表示 x_{i} 來自第k個分量模型，同時變數 z 的邊緣概率p(z_{k}^{i}=1)=alpha_{k} 則表示觀測數據 x_{i} 屬於第k個分量模型的概率。而由於 z 是「1 of K」的形式，因此變數 z 的分布也可以寫成如下形似：

\p ( oldsymbol { z } ) = prod _ { k = 1 } ^ { K }alpha _ { k } ^ { z _ { k } }

相應的，當對 z 指定一個特定的值時，那麼 x 的條件概率分布就是一個特定的高斯分布：

\p ( oldsymbol { x } | z _ { k } = 1 ) = mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , mathbf { Sigma } _ { k } )

也可以寫成：

\p ( oldsymbol { x } | oldsymbol { z } ) = prod _ { k = 1 } ^ { K } mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } ) ^ { z _ { k } }

那麼就可以得到聯合概率分布 p(x,z)=p(z)p(x|z) ，而變數 x 的邊緣概率分布就可以通過聯合概率對所有可能的 z 值進行求和得到，即：

\p ( oldsymbol { x } ) = sum _ { oldsymbol { z } } p ( oldsymbol { z } ) p ( oldsymbol { x } | oldsymbol { z } ) = sum _ { k = 1 } ^ { K } alpha _ { k } mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } )qquad(2)

關於公式（2）第二個等式，乍一看會有些跳躍，但仔細想一想就會知道，所有可能的 z 值其實就是考慮 x 屬於所有K個分量模型的可能性，所以體現的形式就是對所有模型加權求和。同時我們也發現經過一番推導我們得到了最初我們給出的高斯混合模型的表達式，所以我們所定義的變數 z 可以看作是一個隱含變數，用於決定觀測數據屬於哪個模型。

2.3使用EM方法估計GMM的參數

為了方便使用EM方法估計GMM的參數，我們首先定義一個新的量 gamma left( z _ { k } ight) ，用它來表示後驗概率 p left( z _ { k } = 1 | oldsymbol { x } ight) 。根據貝葉斯理論，我們通過先驗概率 p(oldsymbol { z}) 以及條件概率 p(oldsymbol { x }|oldsymbol { z }) 我們可以寫出：

\egin{aligned} gamma left( z _ { k } ight) equiv p left( z _ { k } = 1 | oldsymbol { x } ight) & = frac { p left( z _ { k } = 1 ight) p ( oldsymbol { x } | z _ { k } = 1 ) } { sum _ { j = 1 } ^ { K } p left( z _ { j } = 1 ight) p ( oldsymbol { x } | z _ { j } = 1 ) } \ & = frac { alpha _ { k } mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } ) } { sum _ { j = 1 } ^ { K } alpha _ { j } mathcal { N } ( oldsymbol { x } | oldsymbol { mu } _ { j } , oldsymbol { Sigma } _ { j } ) } end{aligned}qquad(3)

上面的概率公式給出的就是在已知觀測樣本的前提下估計隱藏變數的概率，這是我們在使用EM演算法中的一個重要的過渡參數。

2.3.1極大似然函數

為了估計GMM的參數，我們要估計 alpha,mu,Sigma 這三個參數，所以我們首先將公式（1）的GMM模型改寫成： p ( oldsymbol { x } | oldsymbol { pi } , oldsymbol { mu } , oldsymbol { Sigma } ) 。然後由於 x_{i}， i，2...，n 之間獨立同分布，因此：

\p(X|alpha,mu,Sigma)=prod_{i=1}^{N}p(oldsymbol { x }|alpha,mu,Sigma)qquad(4)

為了估計參數，我們要對其進行最大化似然函數操作，為了計算方便，我們首先給出公式（4）對數似然函數：

\ln p ( oldsymbol { X } | oldsymbol {alpha } , oldsymbol { mu } , mathbf { Sigma } ) = sum _ { n = 1 } ^ { N } ln left{ sum _ { k = 1 } ^ { K }alpha _ { k } mathcal { N } left( oldsymbol { x } _ { n } | oldsymbol { mu } _ { k } , mathbf { Sigma } _ { k } ight) ight}qquad(5)

在最大化似然函數的前提下，我們首先來估計參數 mu ，很自然的做法就是讓公式（5）對變數 mu_{k} 求導，然後令其導數為0，即：

\0 = sum _ { n = 1 } ^ { K } underbrace { frac { alpha _ { k } mathcal { N } left( oldsymbol { x } _ { n } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } ight) } { sum _ { j }alpha_ { j } mathcal { N } left( oldsymbol { x } _ { n } | oldsymbol { mu } _ { j } , oldsymbol { Sigma } _ { j } ight) } } _ { gamma left( z _ { n k } ight) } mathbf { Sigma } _ { k } ^ { - 1 } left( oldsymbol { x } _ { n } - oldsymbol { mu } _ { k } ight)qquad(6)

接著對公式（6）兩邊同乘Sigma_{k} 後化簡可得：

\oldsymbol { mu } _ { k } = frac { 1 } { N _ { k } } sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) oldsymbol { x } _ { n }qquad(7)

其中 N _ { k } = sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) ，（化簡借用公式（3））。

同樣我們對 oldsymbol { Sigma } _ { k } 求導並令其導數為0，可得：

\oldsymbol { Sigma } _ { k } = frac { 1 } { N _ { k } } sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) left( mathbf { x } _ { n } - oldsymbol { mu } _ { k } ight) left( mathbf { x } _ { n } - oldsymbol { mu } _ { k } ight) ^ { mathrm { T } }qquad(8)

最後我們要估計的參數是 alpha_{k} ，但這裡需要注意的是，在一開始我們對權值 alpha 有過限制，即： sum _ { k = 1 } ^ { K } alpha _ { k }=1 ，因此在極大化似然函數時要考慮這一條件。所以我們使用拉格朗日法乘數法來求下面式子的極大值：

\ln p ( mathbf { X } | alpha , mu , Sigma ) + lambda left( sum _ { k = 1 } ^ { K } alpha _ { k } - 1 ight)qquad(9)

然後可得：

\0 = sum _ { n = 1 } ^ { N } frac { mathcal { N } left( mathbf { x } _ { n } | oldsymbol { mu } _ { k } , mathbf { Sigma } _ { k } ight) } { sum _ { j } alpha _ { j } mathcal { N } left( mathbf { x } _ { n } | oldsymbol { mu } _ { j } , mathbf { Sigma } _ { j } ight) } + lambdaqquad(10)

對公式（10）前面分子分母約分後可得： 0=sum_{n=1}^{N}{frac{1}{sum_{j}^{}{alpha_{j}}}}+lambda ，同時又由於 sum _ { k = 1 } ^ { K } alpha _ { k }=1 ，所以可得 lambda=-N 。接著我們對公式（10）兩邊同乘 alpha_{k} ，進行化簡後我們可以得到：

\pi _ { k } = frac { N _ { k } } { N }qquad(11)

2.3.2 EM演算法應用於GMM的演算法步驟

根據（7）、（8）、（11）三式我們已經得到了GMM重要參數的更新迭代公式，這就是EM演算法應用於高斯混合模型的實例，其過程包含E步（計算 gamma left( z _ { k } ight) ）和M步（更新參數 alpha,mu,Sigma ）。接下來我們把具體的演算法過程寫下來：

Input：觀測數據 X ： left{ x_{1},x_{2},...,x_{N} ight}

（1）初始化GMM參數：均值 mu_{k} ，協方差 Sigma_{k} ，權值 alpha_{k}

（2）E步：根據當前模型的參數，計算 分模型k對觀測數據x_{j} 的響應度（responsibility） \gamma left( z _ { n k } ight) = frac { pi _ { k } mathcal { N } left( mathrm { x } _ { n } | oldsymbol { mu } _ { k } , oldsymbol { Sigma } _ { k } ight) } { sum _ { j = 1 } ^ { K } pi _ { j } mathcal { N } left( mathbf { x } _ { n } | oldsymbol { mu } _ { j } , oldsymbol { Sigma } _ { j } ight) } （3）M步: 迭代更新模型的參數值：\egin{aligned} oldsymbol { mu } _ { k } ^ { mathrm { new } } & = frac { 1 } { N _ { k } } sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) mathbf { x } _ { n } \ mathbf { Sigma } _ { k } ^ { mathrm { new } } & = frac { 1 } { N _ { k } } sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) left( mathbf { x } _ { n } - oldsymbol { mu } _ { k } ^ { mathrm { new } } ight) left( mathbf { x } _ { n } - oldsymbol { mu } _ { k } ^ { mathrm { new } } ight) ^ { mathrm { T } } \ pi _ { k } ^ { mathrm { new } } & = frac { N _ { k } } { N } end{aligned} 其中：\N _ { k } = sum _ { n = 1 } ^ { N } gamma left( z _ { n k } ight) （4）：估算對數似然函數的值：\ln p ( mathbf { X } | oldsymbol { mu } , mathbf { Sigma } , oldsymbol { pi } ) = sum _ { n = 1 } ^ { N } ln left{ sum _ { k = 1 } ^ { K } pi _ { k } mathcal { N } left( mathbf { x } _ { n } | oldsymbol { mu } _ { k } , mathbf { Sigma } _ { k } ight) ight} 重複（2）、（3）步直到演算法收斂（似然函數值變化或參數變化小於設定值）

2.3.3 Example

下面我們使用一組真實數據來模擬使用EM演算法求解GMM的過程。

數據：老忠實泉數據 （一組泉水噴發與時間的數據）

我們首先表示出數據的分布：

可以很明顯地看出數據分為兩類，因此我們可以使用包含兩個高斯分布的GMM對其建模。

在實際操作時為了有利計算，往往會把數據進行歸一化處理。

首先是對參數的初始化：

def init_params(shape, K): # K為模型中高斯分布數量
N, D = shape # 樣本數量和數據維度
mu = np.random.rand(K, D)
cov = np.array([np.eye(D)] * K)
alpha = np.array([1.0 / K] * K)

return mu, cov, alpha

然後是E步：

def getExpectation(Y, mu, cov, alpha):
N = Y.shape[0]
K = alpha.shape[0]

# 計算響應度矩陣
gamma = np.zeros((N, K))

for k in range(K):
gamma[:, k] = phi(Y, mu[k], cov[k])
gamma = np.mat(gamma)

# 計算每個模型對每個樣本的響應度
for k in range(K):
gamma[:, k] = alpha[k] * gamma[:, k]
for i in range(N):
gamma[i, :] /= np.sum(gamma[i, :])
return gamma

接著是M步：

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
print(mu)
# 更新 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

然後只需要重複迭代E步和M步達到收斂即可（可以自行設定收斂條件）：

def GMM(Y, K, times):
Y = scale_data(Y)
mu, cov, alpha = init_params(Y.shape, K)
# 收斂條件
delta = 0.000001
for i in range(times):
print("times: {}".format(i))
gamma = Expectation(Y, mu, cov, alpha)
last_mu = mu
last_cov = cov
last_alpha = alpha
mu, cov, alpha = maximize(Y, gamma)
mu_delta = np.sum(np.abs(last_mu - mu))
cov_delta = np.sum(np.abs(last_cov - cov))
alpha_delta = np.sum(np.abs(last_alpha - alpha))
# 判斷是否收斂
if mu_delta < delta and cov_delta < delta and alpha_delta <delta:
break

return mu, cov, alpha

下面是迭代之後得到的結果：

從圖中可以看到原始分布被有效的劃分為了兩個部分。

3、Conclusion

總結來說，EM演算法的思路就是先求對數似然函數的期望值，然後再最大化似然函數以此來更新參數，通過不斷地迭代來找的模型的參數。

對於使用EM演算法比較重要的一點就是要找好隱含變數，通過隱含變數把觀測數據與


推薦閱讀：
相关文章