注：本文是對《統計學習方法》EM算法的一個簡單總結。
1. 什么是EM算法？
引用書上的話：
概率模型有時既含有觀測變量，又含有隱變量或者潛在變量。如果概率模型的變量都是觀測變量，可以直接使用極大似然估計法或者貝葉斯的方法進行估計模型參數，但是當模型含有隱藏變量時，就不能簡單使用這些方法了。EM算法就是含有隱變量的概率模型參數的極大似然估計法，或者極大似然后驗概率估計法。
2. EM 算法的一個小例子:三硬幣模型
假設有3枚硬幣，記作A,B,C。這些硬幣的正面出現的概率分別為\(\pi\)、\(p\)、\(q\)。進行如下的試驗：先擲硬幣A，根據A的結果選擇B和C，如果擲A得到正面，則選擇B；如果擲A得到反面，則選擇C。接着擲出選出的硬幣。記錄下這次擲硬幣的結果，如果是正面，則記作1，反面則記作0。獨立重復做了n次試驗(這里取n=10)，得到結果如下：1,1,0,1,0,0,1,0,1,1。假設只能觀測到拋硬幣的結果，不能觀測到拋硬幣的過程，那么我們該如何估計三硬幣的參數\(\pi\)、\(p\)、\(q\)呢？（也就是估計三枚硬幣正面向上的概率）
EM算法分為E步和M步。EM 算法首先選取了參數的初始值，記作\(\theta^{(0)}\)=(\(\pi^{(0)}\),\(p^{(0)}\),\(q^{(0)}\))。然后通過下面的步驟迭代計算參數的估計值，直到收斂為止，第\(i\)次迭代的參數的估計值記作\(\theta^{(i)}\)=(\(\pi^{(i)}\),\(p^{(i)}\),\(q^{(i)}\)),則EM算法的第\(i+1\)次迭代為：
E步：計算模型在參數\(\pi^{(i)}\)、\(p^{(i)}\)、\(q^{(i)}\)下觀測數據\(y_j\)來自擲硬幣B的概率為
\(\mu_j^{(i+1)}\) = \(\frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j} + (1-\pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}\)
M步：計算模型新的參數的估計值：
\(\pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^{n}\mu_j^{(i+1)}\)
\(p^{(i+1)}=\frac{\sum_{j=1}^{n}\mu_j^{(i+1)}y_j}{\sum_{j=1}^{n}\mu_j^{(i+1)}}\)
\(q^{(i+1)}=\frac{\sum_{j=1}^{n}(1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^{n}(1-\mu_j^{(i+1)})}\)
下面帶入具體的數字計算一下。如果模型的初始參數取值為:\(\pi{(0)}=0.5,p^{(0)}=0.5,q^{(0)}=0.5\)
那么帶入上面的三個公式就可以計算出:\(\pi{(1)}=0.5,p^{(1)}=0.6,q^{(1)}=0.6\),繼續迭代可以得到 \(\pi{(2)}=0.5,p^{(2)}=0.6,q^{(2)}=0.6\)。於是我們可以得到原來參數\(\theta\)的極大似然估計為:\(\hat{\pi}=0.5,\hat{p}=0.6,\hat{q}=0.6\)。
如果取初值改為 \(\pi{(0)}=0.4,p^{(0)}=0.6,q^{(0)}=0.7\)，那么我們可以計算出:\(\mu_1^{(1)}=\frac{0.4\times0.6}{0.4\times0.6+0.6\times0.7}=0.364\)
:\(\mu_0^{(1)}=\frac{0.4\times0.4}{0.4\times0.4+0.6\times0.3}=0.47\)
那么由於觀測值中1出現了6次，0出現了4次，於是我們容易計算出：
\(\pi^{(1)}=\frac{0.364\times6+0.47\times4}{10}=0.4064\)
\(p^{(1)}=\frac{0.364\times1\times6+0.47\times0\times4}{0.364\times6+0.47\times4}=0.537\)
\(q^{(1)}=\frac{(1-0.364)\times6}{(1-0.364)\times6+(1-0.47)\times4}=0.643\)
顯然，我們由於初始值選擇的不同，導致了模型參數最終估計的不同，這說明了EM算法對於參數的初始值是非常敏感的。
3 EM 算法的步驟
輸入:觀測變量數據\(Y\)，隱變量數據\(Z\),聯合分布\(P(Y,Z|\theta)\)，條件分布\(P(Z|Y,\theta)\);
輸出：模型參數\(\theta\)。
(1) 選擇模型參數的初值\(\theta^{(0)}\)，開始迭代；
(2)\(E\)步:記\(\theta^{(i)}\)為第\(i\)次迭代參數\(\theta\)的估計值，在第\(i+1\)次迭代的E步，計算
\(Q(\theta,\theta^{(i)}=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)
這里\(P(Z|Y,\theta^{(i)})\)是在給定觀測數據Y和當前參數估計的\(\theta^{(i)}\)下隱變量數據\(Z\)的條件概率分布；
(3)\(M\)步:求使得\(Q(\theta,\theta^{(i)})\)極大化的\(\theta\)，確定第\(i+1\)次迭代的參數的估計值\(\theta^{(i+1)}\)
\(\theta^{(i+1)}=arg max_{\theta}Q(\theta,\theta^{(i)})\)
(4)重復第(2)步和第(3)步，直到算法收斂。
上面定義的函數\(Q(\theta,\theta^{(i)})\)是EM算法的核心，稱為Q函數(Q function)。下面給出Q函數的定義：
(Q函數) 完全數據的對數似然函數\(logP(Y,Z|\theta)\)關於給定觀測數據\(Y\)和當前參數\(\theta^{(i)}\)下對未觀測數據\(Z\)的條件概率分布\(P(Z|Y,\theta^{(i)})\)的期望稱為\(Q\)函數，即：
\(Q(\theta,\theta^{(i)}=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]\)
4. EM算法的收斂性
1. 設\(P(Y|\theta)\)為觀測數據的似然函數,\(\theta^{(i)}(i=1,2,...)\)為EM算法得到的參數估計序列，\(P(Y|\theta^{(i)}\)為對應的似然函數序列，則\(P(Y|\theta^{(i)}\)是單調遞增的，即：
\(P(Y|\theta^{(i+1)} \ge P(Y|\theta^{(i)})\)。
2. 設 \(L(\theta)=logP(Y|\theta)\)為觀測數據的對數似然函數，\(\theta^{(i)}\)為EM算法得到的參數估計序列,\(L(\theta^{(i)})\)為對應的對數似然函數序列，則
(1) 如果\(P(Y|\theta)\)有上界，則\(L(\theta^{(i)})=logP(Y|\theta^{(i)})\)收斂到某一值\(L^{*}\);
(2) 在函數 \(Q(\theta,\theta^{'})\)與\(L(\theta)\)滿足一定的條件下，由EM算法得到的參數估計序列\(\theta^{(i)}\)的收斂值\(\theta^{*}\)是\(L(\theta)\)的穩定點。
5. 高斯混合模型(GMM)
1. 高斯混合模型的定義:高斯混合模型是指具有如下形式的概率分布模型：
\(P(Y|\theta)=\sum_{k=1}^{K}\alpha_k\phi(y|\theta_k)\)
其中\(\alpha_k\)是系數，\(\alpha_k \ge 0\),\(\sum_{k=1}^{K}\alpha_k=1\);\(\phi(y|\theta_k)\)是高斯分布(正態分布)密度函數，\(\theta_k=(\mu_k,\sigma_k^{2})\),
\(\phi(y|\theta_k)=\frac{1}{\sqrt(2\pi)\sigma_k}exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})\) 稱為第\(k\)個分模型。
2. 使用EM算法估計高斯混合模型
具體推導過程從略，可以參見《統計學習方法》。這里直接給出結果：
高斯混合模型的EM估計算法
輸入：觀測數據\(y_1,y_2,...,y_N\)，高斯混合模型；
輸出：高斯混合模型的參數。
(1)取參數的初始值迭代
(2)E步，依據當前模型的參數，計算分模型k對觀測數據\(y_j\)的響應度，
\(\hat\gamma_{jk}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k\phi(y_j|\theta_k)},j=1,2,...,N,k=1,2,...K\)
(3)計算新一輪迭代的模型參數
\(\hat\mu_k=\frac{\sum_{j=1}^{N}\hat\gamma_{jk}y_j}{\sum_{j=1}^{N}\hat\gamma_{jk}},k=1,2,...,K\)
\(\hat\sigma_k^2=\frac{\sum_{j=1}^{N}\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^{N}\hat\gamma_{jk}},k=1,2,...,K\)
\(\hat\alpha_k=\frac{\sum_{j=1}^{N}\hat\gamma_{jk}}{N},k=1,2,...,K\)
(4)重復(2)和(3)直到收斂。
6. python 實現 EM算法求解高斯混合模型
求解的步驟參考5，這里直接給出代碼如下：
from __future__ import print_function
import numpy as np
def generateData(k,mu,sigma,dataNum):
'''
產生混合高斯模型的數據
:param k: 比例系數
:param mu: 均值
:param sigma: 標准差
:param dataNum:數據個數
:return: 生成的數據
'''
# 初始化數據
dataArray = np.zeros(dataNum,dtype=np.float32)
# 逐個依據概率產生數據
# 高斯分布個數
n = len(k)
for i in range(dataNum):
# 產生[0,1]之間的隨機數
rand = np.random.random()
Sum = 0
index = 0
while(index < n):
Sum += k[index]
if(rand < Sum):
dataArray[i] = np.random.normal(mu[index],sigma[index])
break
else:
index += 1
return dataArray
def normPdf(x,mu,sigma):
'''
計算均值為mu，標准差為sigma的正態分布函數的密度函數值
:param x: x值
:param mu: 均值
:param sigma: 標准差
:return: x處的密度函數值
'''
return (1./np.sqrt(2*np.pi))*(np.exp(-(x-mu)**2/(2*sigma**2)))
def em(dataArray,k,mu,sigma,step = 10):
'''
em算法估計高斯混合模型
:param dataNum: 已知數據個數
:param k: 每個高斯分布的估計系數
:param mu: 每個高斯分布的估計均值
:param sigma: 每個高斯分布的估計標准差
:param step:迭代次數
:return: em 估計迭代結束估計的參數值[k,mu,sigma]
'''
# 高斯分布個數
n = len(k)
# 數據個數
dataNum = dataArray.size
# 初始化gama數組
gamaArray = np.zeros((n,dataNum))
for s in range(step):
for i in range(n):
for j in range(dataNum):
Sum = sum([k[t]*normPdf(dataArray[j],mu[t],sigma[t]) for t in range(n)])
gamaArray[i][j] = k[i]*normPdf(dataArray[j],mu[i],sigma[i])/float(Sum)
# 更新 mu
for i in range(n):
mu[i] = np.sum(gamaArray[i]*dataArray)/np.sum(gamaArray[i])
# 更新 sigma
for i in range(n):
sigma[i] = np.sqrt(np.sum(gamaArray[i]*(dataArray - mu[i])**2)/np.sum(gamaArray[i]))
# 更新系數k
for i in range(n):
k[i] = np.sum(gamaArray[i])/dataNum
return [k,mu,sigma]
if __name__ == '__main__':
# 參數的准確值
k = [0.3,0.4,0.3]
mu = [2,4,3]
sigma = [1,1,4]
# 樣本數
dataNum = 5000
# 產生數據
dataArray = generateData(k,mu,sigma,dataNum)
# 參數的初始值
# 注意em算法對於參數的初始值是十分敏感的
k0 = [0.3,0.3,0.4]
mu0 = [1,2,2]
sigma0 = [1,1,1]
step = 6
# 使用em算法估計參數
k1,mu1,sigma1 = em(dataArray,k0,mu0,sigma0,step)
# 輸出參數的值
print("參數實際值:")
print("k:",k)
print("mu:",mu)
print("sigma:",sigma)
print("參數估計值:")
print("k1:",k1)
print("mu1:",mu1)
print("sigma1:",sigma1)
運行結果如下：
如果改變初始參數為:
k0=[0.1,0.4,0.5]
mu0=[3,3,2]
sigma0=[2,1,1.5]
那么得到的結果為：
對比可以發現，兩次參數估計的結果相差還是很大的。這也說明了EM算法對初始參數的敏感性。在實際應用的過程中，我們可以通過多次試驗驗證來選擇估計結果最好的那個作為初始參數。