EM(最大期望)演算法推導、GMM的應用與程式碼實現

頎周發表於2020-06-21

  EM演算法是一種迭代演算法,用於含有隱變數的概率模型引數的極大似然估計。

使用EM演算法的原因

  首先舉李航老師《統計學習方法》中的例子來說明為什麼要用EM演算法估計含有隱變數的概率模型引數。

  假設有三枚硬幣,分別記作A, B, C。這些硬幣正面出現的概率分別是$\pi,p,q$。進行如下擲硬幣試驗:先擲硬幣A,根據其結果選出硬幣B或C,正面選硬幣B,反面邊硬幣C;然後擲選出的硬幣,擲硬幣的結果出現正面記作1,反面記作0;獨立地重複$n$次試驗,觀測結果為$\{y_1,y_2,...,y_n\}$。問三硬幣出現正面的概率。

  三硬幣模型(也就是第二枚硬幣正反面的概率)可以寫作

$ \begin{aligned} &P(y|\pi,p,q) \\ =&\sum\limits_z P(y,z|\pi,p,q)\\ =&\sum\limits_z P(y|z,\pi,p,q)P(z|\pi,p,q)\\ =&\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} \end{aligned} $

  其中$z$表示硬幣A的結果,也就是前面說的隱變數。通常我們直接使用極大似然估計,即最大化似然函式

$ \begin{aligned} &\max\limits_{\pi,p,q}\prod\limits_{i=1}^n P(y_i|\pi,p,q) \\ =&\max\limits_{\pi,p,q}\prod\limits_{i=1}^n[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}\sum\limits_{i=1}^n\log[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}L(\pi,p,q) \end{aligned} $

  分別對$\pi,p,q$求偏導並等於0,求解線性方程組來估計這三個引數。但是,由於它是帶有隱變數的,在獲取最終的隨機變數之前有一個分支選擇的過程,導致這個$\log$的內部是加和的形式,計算導數十分困難,而待求解的方程組不是線性方程組。當複雜度一高,解這種方程組幾乎成為不可能的事。以下推導EM演算法,它以迭代的方式來求解這些引數,應該也算一種貪心吧。

演算法匯出與理解

  對於引數為$\theta$且含有隱變數$Z$的概率模型,進行$n$次抽樣。假設隨機變數$Y$的觀察值為$\mathcal{Y} = \{y_1,y_2,...,y_n\}$,隱變數$Z$的$m$個可能的取值為$\mathcal{Z}=\{z_1,z_2,...,z_m\}$。

  寫出似然函式:

$ \begin{aligned} L(\theta) &= \sum\limits_{Y\in\mathcal{Y}}\log P(Y|\theta)\\ &=\sum\limits_{Y\in\mathcal{Y}}\log \sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ \end{aligned} $

  EM演算法首先初始化引數$\theta = \theta^0$,然後每一步迭代都會使似然函式增大,即$L(\theta^{k+1})\ge L(\theta^k)$。如何做到不斷變大呢?考慮迭代前的似然函式(為了方便不用$\theta^{k+1}$):

$ \begin{gather} \begin{aligned} L(\theta)=&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ =&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)}\\ \end{aligned} \label{} \end{gather} $

  至於上式的第二個等式為什麼取出$P(Z|Y,\theta^k)$而不是別的,正向的原因我想不出來,馬後炮原因在後面記錄。

  考慮其中的求和

$ \sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)=1$

  且由於$\log$函式是凹函式,因此由Jenson不等式得

$ \begin{gather} \begin{aligned} L(\theta) \ge&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)}\\ =&B(\theta,\theta^k) \end{aligned}\label{} \end{gather} $

  當$\theta = \theta^k$時,有

$ \begin{gather} \begin{aligned} L(\theta^k) \ge& B(\theta^k,\theta^k)\\ =&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta^k)}{P(Z|Y,\theta^k)}\\ =&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log P(Y|\theta^k)\\ =&\sum\limits_{Y\in \mathcal{Y}}\log P(Y|\theta^k)\\ =&L(\theta^k)\\ \end{aligned} \label{} \end{gather} $

  也就是在這時,$(2)$式取等,即$L(\theta^k) = B(\theta^k,\theta^k)$。取

$ \begin{gather} \theta^*=\text{arg}\max\limits_{\theta}B(\theta,\theta^k)\label{} \end{gather} $

  可得不等式

$L(\theta^*)\ge B(\theta^*,\theta^k)\ge B(\theta^k,\theta^k) = L(\theta^k)$

  所以,我們只要優化$(4)$式,讓$\theta^{k+1} = \theta^*$,即可保證每次迭代的非遞減勢頭,有$L(\theta^{k+1})\ge L(\theta^k)$。而由於似然函式是概率乘積的對數,一定有$L(\theta) < 0$,所以迭代有上界並且會收斂。以下是《統計學習方法》中EM演算法一次迭代的示意圖:

  進一步簡化$(4)$式,去掉優化無關項:

$ \begin{aligned} \theta^*=&\text{arg}\max\limits_{\theta}B(\theta,\theta^k) \\ =&\text{arg}\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)} \\ =&\text{arg}\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log P(Y,Z|\theta) \\ =&\text{arg}\max\limits_{\theta}Q(\theta,\theta^k) \\ \end{aligned} $

  $Q$函式使用導數求極值的方程與沒有隱變數的方程類似,容易求解。

  綜上,EM演算法的流程為:

  1. 設定$\theta^0$的初值。EM演算法對初值是敏感的,不同初值迭代出來的結果可能不同。

  2. 更新$\theta^k = \text{arg}\max\limits_{\theta}Q(\theta,\theta^{k-1})$。理解上來說,通常將這一步分為計算$Q$與極大化$Q$兩步,即求期望E與求極大M,但在程式碼中並不會將它們分出來,因此這裡濃縮為一步。另外,如果這個優化很難計算的話,因為有不等式的保證,直接取$\theta^k$為某個$\hat{\theta}$,只要有$Q(\hat{\theta},\theta^{k-1})\ge Q(\theta^{k-1},\theta^{k-1})$即可。

  3. 比較$\theta^k$與$\theta^{k-1}$的差異,比如求它們的差的二範數,若小於一定閾值就結束迭代,否則重複步驟2。

  下面記錄一下我對$(1)$式取出$P(Z|Y,\theta^k)$而不取別的$P$的理解:

  經過以上的推導,我認為這是為了給不等式取等創造條件。如果不能確定$L(\theta^k)$與$Q(\theta^k,\theta^k)$能否取等,那麼取$Q$的最大值$Q(\theta^*,\theta^k)$時,儘管有$Q(\theta^*,\theta^k)\ge Q(\theta^k,\theta^k)$,但並不能保證$L(\theta^*)\ge L(\theta^k)$,迭代的不減性質就就沒了。

  我這裡暫且把它看做一種巧合,是研究EM演算法的大佬,碰巧想用Jenson不等式來迭代而構造出來的一種做法。本人段位還太弱,無法正向理解其中的緣故,只能以這種方式來揣度大佬的思路了。知乎大佬發的EM演算法九層理解(點選連結),我當前只能到第3層,有時間一定要拜讀一下深度學習之父的著作。

高斯混合模型的應用

迭代式推導

  假設高斯混合模型混合了$m$個高斯分佈,引數為$\theta = (\alpha_1,\theta_1,\alpha_2,\theta_2,...,\alpha_m,\theta_m),\theta_i=(\mu_i,\sigma_i)$則整個概率分佈為:

$\displaystyle P(y|\theta) = \sum\limits_{i=1}^m\alpha_i \phi(y|\theta_i) =  \sum\limits_{i=1}^m\frac{\alpha_i }{\sqrt{2\pi}\sigma_i}\exp\left(-\frac{(y-\mu_i)^2}{2\sigma_i^2}\right),\;\text{where}\;\sum\limits_{j=1}^m\alpha_j = 1$

  對混合分佈抽樣$n$次得到$\{y_1,...,y_n\}$,則在第$k+1$次迭代,待優化式為:

$\begin{gather}\begin{aligned} &\max\limits_{\theta}Q(\theta,\theta^k) \\ =&\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log P(Y,Z|\theta) \\ =&\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} \frac{P(Z,Y|\theta^k)}{P(Y|\theta^k)}\log P(Y,Z|\theta) \\ =&\max\limits_{\theta}\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)} \log \left[\alpha_j\phi(y_i|\theta_j)\right] \\ =&\max\limits_{\theta}\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)} \log \left[ \frac{\alpha_j}{\sqrt{2\pi}\sigma_j}\exp\left(-\frac{(y_i-\mu_j)^2}{2\sigma_j^2}\right) \right]\\ =&\max\limits_{\theta}\sum\limits_{j=1}^m \sum\limits_{i=1}^n \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)} \left[ \log \alpha_j - \log \sigma_j-\frac{(y_i-\mu_j)^2}{2\sigma_j^2} \right]\\  \end{aligned} \label{}\end{gather}$

計算α

  定義

$\displaystyle n_j = \sum\limits_{i=1}^n \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)}$

  則對於$\alpha$,優化式為

$\begin{gather} \begin{aligned} \max\limits_{\alpha}\sum\limits_{j=1}^m n_j \log \alpha_j \end{aligned} \label{}\end{gather}$

  又因為$\sum\limits_{j=1}^m \alpha_j=1$,所以只需優化$m-1$個引數,上式變為:

$ \max\limits_\alpha \left[ \begin{matrix} n_1&n_2&\cdots &n_{m-1}&n_{m}\\ \end{matrix} \right] \cdot \left[ \begin{matrix} \log\alpha_1\\ \log\alpha_2\\ \vdots\\ \log\alpha_{m-1}\\ \log(1-\alpha_1-\cdots-\alpha_{m-1})\\ \end{matrix} \right] $

  對每個$\alpha_j$求導並等於0,得到線性方程組:

$\left[\begin{matrix}n_1+n_m&n_1&n_1&\cdots&n_1\\n_2&n_2+n_m&n_2&\cdots&n_2\\n_3&n_3&n_3+n_m&\cdots&n_3\\&&&\vdots&\\n_{m-1}&n_{m-1}&n_{m-1}&\cdots&n_{m-1}+n_m\\\end{matrix}\right]\cdot\left[\begin{matrix}\alpha_1\\\alpha_2\\\alpha_3\\\vdots\\\alpha_{m-1}\\\end{matrix}\right]=\left[\begin{matrix}n_1\\n_2\\n_3\\\vdots\\n_{m-1}\\\end{matrix}\right]$

  求解這個爪形線性方程組,得到

$\left[\begin{matrix}\sum_{j=1}^mn_j/n_1&0&0&\cdots&0\\-n_2/n_1&1&0&\cdots&0\\-n_3/n_1&0&1&\cdots&0\\&&&\vdots&\\-n_{m-1}/n_1&0&0&\cdots&1\\\end{matrix}\right]\cdot\left[\begin{matrix}\alpha_1\\\alpha_2\\\alpha_3\\\vdots\\\alpha_{m-1}\\\end{matrix}\right]=\left[\begin{matrix}1\\0\\0\\\vdots\\0\\\end{matrix}\right]$

  因為

$\displaystyle \sum\limits_{j=1}^m n_j =   \sum\limits_{j=1}^m\sum\limits_{i=1}^n \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)}=\sum\limits_{i=1}^n \sum\limits_{j=1}^m \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)} =\sum\limits_{i=1}^n 1 =  n$

  解得

$\displaystyle\alpha_j = \frac{n_j}{n} = \frac{1}{n}\sum\limits_{i=1}^n \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)}$

計算σ與μ

  與$\alpha$不同,它的方程組是所有$\alpha_j$之間聯立的;而$\sigma,\mu$的方程組則是$\sigma_j$與$\mu_j$之間聯立的。定義

$\displaystyle p_{ji} = \frac{\alpha_j^k\phi(y_i|\theta_j^k)} {\sum\limits_{l=1}^m \alpha_l^k\phi(y_i|\theta_l^k)}$

  則對於$\sigma_j,\mu_j$,優化式為(比較$(6),(7)$式的區別)

$\begin{gather}\displaystyle\min\limits_{\sigma_j,\mu_j}\sum\limits_{i=1}^n p_{ji} \left(\log \sigma_j+\frac{(y_i-\mu_j)^2}{2\sigma_j^2} \right)\label{}\end{gather}$

  對上式求導等於0,解得

$ \begin{aligned} &\mu_j = \frac{\sum\limits_{i=1}^np_{ji}y_i}{\sum\limits_{i=1}^np_{ji}} = \frac{\sum\limits_{i=1}^np_{ji}y_i}{n_j} = \frac{\sum\limits_{i=1}^np_{ji}y_i}{n\alpha_j}\\ &\sigma^2_j = \frac{\sum\limits_{i=1}^np_{ji}(y_i-\mu_j)^2}{\sum\limits_{i=1}^np_{ji}} = \frac{\sum\limits_{i=1}^np_{ji}(y_i-\mu_j)^2}{n_j} = \frac{\sum\limits_{i=1}^np_{ji}(y_i-\mu_j)^2}{n\alpha_j} \end{aligned} $

程式碼實現

  對於概率密度為$P(x) = −2x+2,x\in (0,1)$的隨機變數,以下程式碼實現GMM對這一概率密度的的擬合。共10000個抽樣,GMM混合了100個高斯分佈。

#%%定義引數、函式、抽樣
import numpy as np
import matplotlib.pyplot as plt

dis_num = 100 #用於擬合的分佈數量
sample_num = 10000 #用於擬合的分佈數量
alphas = np.random.rand(dis_num) 
alphas /= np.sum(alphas)  
mus = np.random.rand(dis_num)
sigmas = np.random.rand(dis_num)**2#方差,不是標準差
samples = 1-(1-np.random.rand(sample_num))**0.5 #樣本
C_pi = (2*np.pi)**0.5

dis_val = np.zeros([sample_num,dis_num])    #每個樣本在每個分佈成員上都有值,形成一個sample_num*dis_num的矩陣
pij = np.zeros([sample_num,dis_num])        #pij矩陣
def calc_dis_val(sample,alpha,mu,sigma,c_pi):
    return alpha*np.exp(-(sample[:,np.newaxis]-mu)**2/(2*sigma))/(c_pi*sigma**0.5) 
def calc_pij(dis_v):  
    return dis_v / dis_v.sum(axis = 1)[:,np.newaxis]      
#%%優化 
for i in range(1000):
    print(i)
    dis_val = calc_dis_val(samples,alphas,mus,sigmas,C_pi)
    pij = calc_pij(dis_val)  
    nj = pij.sum(axis = 0)
    alphas_before = alphas
    alphas = nj / sample_num
    mus = (pij*samples[:,np.newaxis]).sum(axis=0)/nj
    sigmas = (pij*(samples[:,np.newaxis] - mus)**2 ).sum(axis=0)/nj
    a = np.linalg.norm(alphas_before - alphas)
    print(a)
    if  a< 0.001:
        break

#%%繪圖 
plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標籤
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負號
def get_dis_val(x,alpha,sigma,mu,c_pi):
    y = np.zeros([len(x)]) 
    for a,s,m in zip(alpha,sigma,mu):   
        y += a*np.exp(-(x-m)**2/(2*s))/(c_pi*s**0.5)   
    return y
def paint(alpha,sigma,mu,c_pi,samples):
    x = np.linspace(-1,2,500)
    y = get_dis_val(x,alpha,sigma,mu,c_pi) 
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.hist(samples,density = True,label = '抽樣分佈') 
    ax.plot(x,y,label = "擬合的概率密度")
    ax.legend(loc = 'best')
    plt.show()
paint(alphas,sigmas,mus,C_pi,samples)

  以下是擬合結果圖,有點像是核函式估計,但是完全不同:

EM演算法的推廣

  EM演算法的推廣是對EM演算法的另一種解釋,最終的結論是一樣的,它可以使我們對EM演算法的理解更加深入。它也解釋了我在$(1)$式下方提出的疑問:為什麼取出$P(Z|Y,\theta^k)$而不是別的。

  定義$F$函式,即所謂Free energy自由能(自由能具體是啥先不研究了):

$ \begin{aligned} F(\tilde{P},\theta) &= E_{\tilde{P}}(\log P(Y,Z|\theta)) + H(\tilde{P})\\ &= \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log P(Y,Z|\theta) - \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log \tilde{P}(Z)\\ \end{aligned} $

  其中$\tilde{P}$是$Z$的某個概率分佈(不一定是單獨的分佈,可能是在某個條件下的分佈),$E_{\tilde{P}}$表示分佈$\tilde{P}$下的期望,$H$表示資訊熵。

  我們計算一下,對於固定的$\theta$,什麼樣的$\tilde{P}$會使$F(\tilde{P},\theta) $最大。也就是找到一個函式$\tilde{P}_{\theta}$,使$F$極大,寫成優化的形式就是(這裡是找函式而不是找引數哦,理解上可能要用到泛函分析的內容):

$ \begin{aligned} &\max\limits_{\tilde{P}} \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log P(Y,Z|\theta) - \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log \tilde{P}(Z)\\ &\;\text{s.t.}\; \sum\limits_{Z\in \mathcal{Z}}\tilde{P}(Z) = 1 \end{aligned} $

  拉格朗日函式(拉格朗日對偶性,點選連結)為:

$ \begin{aligned} L =  \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log P(Y,Z|\theta) - \sum\limits_{Z\in \mathcal{Z}} \tilde{P}(Z)\log \tilde{P}(Z)+ \lambda\left(1-\sum\limits_{Z\in \mathcal{Z}}\tilde{P}(Z)\right) \end{aligned} $

  因為每個$\tilde{P}(Z)$之間都是求和,沒有其它其它諸如乘積的操作,所以可以直接令$L$對某個$\tilde{P}(Z)$求導等於$0$來計算極值:

$ \begin{aligned} \frac{\partial L}{\partial \tilde{P}(Z)} = \log P(Y,Z|\theta) - \log \tilde{P}(Z) -1 -\lambda = 0 \end{aligned} $

  於是可以推出:

$ \begin{aligned} P(Y,Z|\theta) = e^{1+\lambda}\tilde{P}(Z) \end{aligned} $

  又由約束$\sum\limits_{Z\in \mathcal{Z}}\tilde{P}(Z) = 1$:

$P(Y|\theta) = e^{1+\lambda}$

  於是得到

$\begin{gather}\tilde{P}_{\theta}(Z) = P(Z|Y,\theta)\label{}\end{gather}$

  代回$F(\tilde{P},\theta)$,得到

$ \begin{aligned} F(\tilde{P}_\theta,\theta) &= \sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta)\log P(Y,Z|\theta) - \sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta)\log P(Z|Y,\theta)\\ &= \sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta)\log \frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}\\ &= \log P(Y|\theta)\\ \end{aligned} $

  也就是說,對$F$關於$\tilde{P}$進行最大化後,$F$就是待求分佈的對數似然;然後再關於$\theta$最大化,也就算得了最終要估計的引數$\hat{\theta}$。所以,EM演算法也可以解釋為$F$的極大-極大演算法。優化結果$(8)$式也解釋了我之前在$(1)$式下方的提問。

  那麼,怎麼使用$F$函式進行估計呢?還是要用迭代來算,迭代方式是和前面介紹的一樣的(懶得記錄了,統計學習方法上直接看吧)。實際上,$F$函式的方法只是提供了EM演算法的另一種解釋,具體方法上並沒有提升之處。

相關文章