極大似然估計
考慮一個高斯分布\(p(\mathbf{x}\mid{\theta})\)，其中\(\theta=(\mu,\Sigma)\)。樣本集\(X=\{x_1,...,x_N\}\)中每個樣本都是獨立的從該高斯分布中抽取得到的，滿足獨立同分布假設。
因此，取到這個樣本集的概率為：
\[\begin{aligned} p(X\mid{\theta}) &= \prod_{i=1}^Np(x_i\mid\theta) \end{aligned}\]
我們要估計模型參數向量\(\theta\)的值，由於現在樣本集\(X\)已知，所以可以把\(p(X\mid{\theta})\)看作是參數向量\(\theta\)的函數，稱之為樣本集\(X\)下的似然函數\(L(\theta\mid{X})\)。
\[\begin{aligned} L(\theta\mid{X}) &= \prod_{i=1}^Np(x_i\mid\theta) \end{aligned}\]
極大似然估計要做的就是最大化該似然函數，找到能夠使得\(p(X\mid{\theta})\)最大的參數向量\(\theta\)，換句話說，就是找到最符合當前觀測樣本集\(X\)的模型的參數向量\(\theta\)
為方便運算，通常我們使用似然函數的對數形式（可以將乘積運算轉化為求和形式），稱之為“對數似然函數”。由於底數大於1的對數函數總是單調遞增的，所以使對數似然函數達到最大值的參數向量也能使似然函數本身達到最大值。故，對數似然函數為：
\[\begin{aligned} L(\theta\mid{X}) &= \log\Big(\prod_{i=1}^Np(x_i\mid\theta)\Big) \\ &= \sum_{i=1}^N\log{p(x_i\mid\theta)} \end{aligned}\]
參數的估計值 \(\hat{\theta}=\arg\underset{\theta}{\max}L(\theta\mid{X})\)
混合高斯模型及其求解困境
混合高斯模型是指由\(k\)個高斯模型疊加形成的一種概率分布模型，形式化表示為：
\[p(\mathbf{x}\mid{\Theta}) = \sum_{l=1}^{k}\alpha_lp_l(\mathbf{x}\mid{\theta_l}) \]
其中，參數 \(\Theta=(\alpha_1,...,\alpha_k,\theta_1,...,\theta_k)\)，並且有\(\Sigma_{l=1}^{k}\alpha_l=1\)，\(\alpha_l\)代表單個高斯分布在混合高斯模型中的權重。
現在我們假設觀測樣本集\(X=(x_1,...x_N)\) 來自於該混合高斯模型，滿足獨立同分布假設。為了估計出該混合高斯模型的參數\(\Theta\)，我們寫出這n個數據的對數似然函數：
\[\begin{aligned} L(\Theta\mid{X}) &= \log\Big(p(X\mid{\Theta})\Big) \\ &= \log\Big(\prod_{i=1}^{N}p(x_i\mid{\Theta})\Big) \\ &= \sum_{i=1}^{N}\log\Big(p(x_i\mid{\Theta})\Big) \\ &= \sum_{i=1}^{N}\log\Big(\sum_{l=1}^{k}\alpha_lp_l(x_i\mid{\theta_l})\Big) \end{aligned}\]
我們的目標就是要通過最大化該似然函數從而估計出混合高斯模型的參數 \(\Theta\)
觀察該式，由於對數函數中還包含求和式，因此如果仿照極大似然估計單個高斯分布參數的方法來求解這個問題，是無法得到解析解的。所以我們要尋求更好的方式來解決這個問題。
EM算法
基本過程
考慮一個樣本集\(X\)是從某種分布（參數為\(\Theta\)）中觀測得到的，我們稱之為不完全數據(incomplete data)。我們引入一個無法直接觀測得到的隨機變量集合\(Z\)，叫作隱變量。X和Z連在一起稱作完全數據。我們可以得到它們的聯合概率分布為：
\[p(X,Z\mid{\Theta})=p(Z\mid{X,\Theta})p(X\mid{\Theta}) \]
我們定義一個新的似然函數叫作完全數據似然：
\[L(\Theta\mid{X,Z})=p(X,Z\mid{\Theta}) \]
我們可以通過極大化這樣一個似然函數來估計參數Theta。然而Z是隱變量，其分布未知上式無法直接求解，我們通過計算完全數據的對數似然函數關於Z的期望，來最大化已觀測數據的邊際似然。我們定義這樣一個期望值為Q函數：
\[\begin{aligned} Q(\Theta,\Theta^{g}) &= \mathbb{E}_Z\Big[\log{p(X,Z\mid\Theta)}\mid{X,\Theta^g}\Big]\\ &= \sum_Z\log{p(X,Z\mid{\Theta})p(Z\mid{X,\Theta^g})} \end{aligned}\]
其中，\(\Theta^g\)表示當前的參數估計值，X是觀測數據，都作為常量。\(\Theta\)是我們要極大化的參數，\(Z\)是來自於分布\(p(Z\mid{X,\Theta^g})\)的隨機變量。\(p(Z\mid{X,\Theta^g})\)是未觀測變量的邊緣分布，並且依賴於觀測數據\(X\)和當前模型參數\(\Theta^g\)
EM算法中的E-setp就是指上面計算期望值的過程。
M-step則是極大化這樣一個期望，從而得到能夠最大化期望的參數\(\Theta\)
\[\begin{aligned} \Theta^{g+1} &= \arg\underset{\theta}{\max}Q(\Theta,\Theta^{g}) \\ &= \arg\underset{\theta}{\max}\sum_Z\log{p(X,Z\mid{\Theta})p(Z\mid{X,\Theta^g})} \end{aligned}\]
重復執行E-step和M-step直至收斂。
收斂性證明
(略)
EM算法應用於高斯混合模型
回到高斯混合模型的參數估計問題，我們將樣本集\(X\)看作不完全數據，考慮一個無法觀測到的隨機變量\(Z=\{z_i\}_{i=1}^{N}\)，其中\(z_i\in\{1,...,k\}\)，用來指示每一個數據點是由哪一個高斯分布產生的，（可以類比成一個類別標簽，這個標簽我們是無法直接觀測得到的）。比如，\(z_i=k\) 表示第\(i\)個樣本點是由混合高斯模型中的第\(k\)個分量模型生成的。
那么完全數據\((X,Z)\)的概率分布為：
\[\begin{aligned} p(X,Z\mid{\Theta}) &= \prod_{i=1}^{N}p(x_i,z_i\mid{\Theta}) \\ &= \prod_{i=1}^{N}p(z_i\mid{\Theta})p(x_i\mid{z_i,\Theta}) \\ &= \prod_{i=1}^{N}\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}}) \end{aligned}\]
在混合高斯模型問題中，\(p(z_i\mid\Theta)\) 是指第\(z_i\)個模型的先驗概率，也就是其權重\(\alpha_{z_i}\)。給定樣本點來源的高斯分布類別后，\(p(x_i\mid{z_i,\Theta})\)可以寫成對應的高斯分布下的概率密度形式，即\(p_{z_i}(x_i\mid{\theta_{z_i}})\)
根據貝葉斯公式，又可以得到隱變量的條件概率分布：
\[\begin{aligned} p(Z\mid{X,\Theta}) &= \prod_{i=1}^Np(z_i\mid{x_i,\Theta}) \\ &= \prod_{i=1}^N\frac{p(x_i,z_i,\Theta)}{p(x_i,\Theta)} \\ &= \prod_{i=1}^N\frac{\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})}{\sum_{z_i=1}^k\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})} \end{aligned}\]
因此，我們可以寫出相應的\(Q\)函數：
\[\begin{aligned} Q(\Theta,\Theta^{g}) &= \sum_Z\log{p(X,Z\mid{\Theta})p(Z\mid{X,\Theta^g})} \\ &= \sum_Z\log{\prod_{i=1}^{N}\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})}\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &= \sum_Z\sum_{i=1}^{N}\log\Big(\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})\Big)\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &= \sum_{z_1=1}^k\sum_{z_2=1}^k...\sum_{z_N=1}^k\sum_{i=1}^{N}\log\Big(\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})\Big)\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &= \sum_{z_1=1}^k\sum_{z_2=1}^k...\sum_{z_N=1}^k\log\Big(\alpha_{z_1}p_{z_1}(x_1\mid{\theta_{z_1}})\Big)\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &\quad+ \sum_{z_1=1}^k\sum_{z_2=1}^k...\sum_{z_N=1}^k\sum_{i=2}^{N}\log\Big(\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})\Big)\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &= \mathcal{A} + \mathcal{B} \end{aligned}\]
其中，
\[\begin{aligned} \mathcal{A} &= \sum_{z_1=1}^k\sum_{z_2=1}^k...\sum_{z_N=1}^k\log\Big(\alpha_{z_1}p_{z_1}(x_1\mid{\theta_{z_1}})\Big)\prod_{i=1}^Np(z_i\mid{x_i,\Theta^g}) \\ &= \sum_{z_1=1}^k\log\Big(\alpha_{z_1}p_{z_1}(x_1\mid{\theta_{z_1}})\Big)p(z_1\mid{x_1,\Theta^g})\underbrace{\sum_{z_2=1}^k...\sum_{z_N=1}^k\prod_{i=2}^Np(z_i\mid{x_i,\Theta^g})}_{result=1} \\ &= \sum_{z_1=1}^k\log\Big(\alpha_{z_1}p_{z_1}(x_1\mid{\theta_{z_1}})\Big)p(z_1\mid{x_1,\Theta^g}) \\ \end{aligned}\]
\(\mathcal{B}\)式可以按照同樣的操作技巧進行分解，故不贅述。
並且，我們用\(l\)代替\(z_i\)來簡化我們的表達式。
因此，\(Q\)函數可以化簡為：
\[\begin{aligned} Q(\Theta,\Theta^{g}) &= \sum_{i=1}^N\sum_{z_i=1}^k\log\Big(\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})\Big)p(z_i\mid{x_i,\Theta^g}) \\ &= \sum_{i=1}^N\sum_{l=1}^k\log\Big(\alpha_{l}p_{l}(x_i\mid{\theta_{l}})\Big)p(l\mid{x_i,\Theta^g}) \\ &= \sum_{l=1}^k\sum_{i=1}^N\log\Big(\alpha_{l}p_{l}(x_i\mid{\theta_{l}})\Big)p(l\mid{x_i,\Theta^g}) \\ &= \sum_{l=1}^k\sum_{i=1}^N\log(\alpha_l)p(l\mid{x_i,\Theta^g})+\sum_{l=1}^k\sum_{i=1}^N\log\Big(p_l(x_i\mid{\theta_{l}})\Big)p(l\mid{x_i,\Theta^g}) \\ \end{aligned}\]
這樣，我們可以對包含參數\(\alpha_l\) 和參數\(\theta_l\) 的項分別進行最大化從而得到各自的估計值。
估計參數\(\alpha_l\)
這個問題可以表示為下面的約束最優化問題：
\[\begin{aligned} \underset{\alpha_l}{\max}&\quad \sum_{l=1}^k\sum_{i=1}^N\log(\alpha_l)p(l\mid{x_i,\Theta^g}) \\ s.t.& \quad\sum_{l=1}^k{\alpha_l}=1 \end{aligned}\]
引入拉格朗日乘子\(\lambda\)，構建拉格朗日函數：
\[\begin{aligned} \mathcal{L}(\alpha_1,...,\alpha_2,\lambda) &= \sum_{l=1}^k\sum_{i=1}^N\log(\alpha_l)p(l\mid{x_i,\Theta^g})-\lambda\Big(\sum_{l=1}^k{\alpha_l}-1\Big) \\ &= \sum_{l=1}^k\log(\alpha_l)\sum_{i=1}^Np(l\mid{x_i,\Theta^g})-\lambda\Big(\sum_{l=1}^k{\alpha_l}-1\Big) \end{aligned}\]
對參數\(\alpha_l\)求偏導並令其為0：
\[\frac{\partial\mathcal{L}}{\partial\alpha_l}=\frac{1}{\alpha_l}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})-\lambda=0 \]
得到：
\[\alpha_l=\frac{1}{\lambda}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g}) \]
代回約束條件，有：
\[1-\frac{1}{\lambda}\sum_{l=1}^{k}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})=0 \]
在之前的推導中，我們得到過這樣一個公式，即隱變量\(z\)的條件分布：
\[\begin{aligned} p(z_i\mid{x_i,\Theta}) &= \frac{p(x_i,z_i,\Theta)}{p(x_i,\Theta)} \\ &= \frac{\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})}{\sum_{z_i=1}^k\alpha_{z_i}p_{z_i}(x_i\mid{\theta_{z_i}})} \end{aligned}\]
同樣用\(l\)代替\(z_i\)來簡化我們的表達式，得到：
\[\begin{aligned} p(l\mid{x_i,\Theta}) &= \frac{p(x_i,l,\Theta)}{p(x_i,\Theta)} \\ &= \frac{\alpha_{l}p_{l}(x_i\mid{\theta_{l}})}{\sum_{l=1}^k\alpha_{l}p_{l}(x_i\mid{\theta_{l}})} \end{aligned}\]
故將其代回之前的式子，得到：
\[\begin{aligned} 1-\frac{1}{\lambda}\sum_{l=1}^{k}\sum_{i=1}^{N}\frac{\alpha_{l}p_{l}(x_i\mid{\theta_{l}})}{\sum_{l=1}^k\alpha_{l}p_{l}(x_i\mid{\theta_{l}})}&=0 \\ 1-\frac{1}{\lambda}\sum_{i=1}^{N}\sum_{l=1}^{k}\frac{\alpha_{l}p_{l}(x_i\mid{\theta_{l}})}{\sum_{l=1}^k\alpha_{l}p_{l}(x_i\mid{\theta_{l}})}&=0 \\ 1-\frac{1}{\lambda}\sum_{i=1}^{N}(1)&=0 \\ 1-\frac{N}{\lambda}&=0 \\ \lambda&=N \end{aligned}\]
故，
\[\alpha_l=\frac{1}{N}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g}) \]
在估計剩下的參數之前，先補充一下一會要用到的線性代數知識。
【知識點~Matrix Algebra】
- 矩陣的跡等於矩陣的主對角線上元素之和，且有以下性質
- \(|A|\)表示矩陣\(A\)的行列式，當\(A\)可逆時，有以下性質
- \(a_{i,j}\)表示矩陣\(A\)的第\(i\)行，\(j\)列的元素，給出矩陣函數求導的一些公式：
-
\[\frac{\partial{x^TAx}}{\partial{x}}=(A+A^T)x \]
-
\[\frac{\partial\log|A|}{\partial{A}}=2A^{-1}-diag(A^{-1}) \]
-
\[\frac{\partial{tr(AB)}}{\partial{A}}=B+B^T-diag(B) \]
估計參數\(\theta_l\)
對於\(d\)維高斯分布來說，參數\(\theta=(\mu,\Sigma)\)，且：
\[p_l(x\mid{\mu_l,\Sigma_l})=\frac{1}{(2\pi)^{d/2}\vert\Sigma_l\vert^{1/2}}\exp\big[-\frac{1}{2}(x-\mu_l)^T\Sigma_l^{-1}(x-\mu_l)\big] \]
估計參數\(\mu_l\)
\[\begin{aligned} &\quad\sum_{l=1}^k\sum_{i=1}^N\log\Big(p_l(x_i\mid{\mu_l,\Sigma_l})\Big)p(l\mid{x_i,\Theta^g}) \\ &= \sum_{l=1}^k\sum_{i=1}^N\Big(\log(2\pi^{-d/2})+\log|\Sigma_l|^{-1/2}-\frac{1}{2}(x_i-\mu_l)^T\Sigma_l^{-1}(x_i-\mu_l)\Big)p(l\mid{x_i,\Theta^g}) \\ \end{aligned}\]
忽略其中的常數項（因為常數項在求導之后為0），上式化簡為：
\[\sum_{l=1}^k\sum_{i=1}^N\Big(-\frac{1}{2}\log|\Sigma_l|-\frac{1}{2}(x_i-\mu_l)^T\Sigma_l^{-1}(x_i-\mu_l)\Big)p(l\mid{x_i,\Theta^g}) \]
關於\(\mu_l\)求導，得到：
\[\sum_{l=1}^k\sum_{i=1}^N\Bigg(-\frac{1}{2}\Big(\Sigma_l^{-1}+(\Sigma_l^{-1})^T\Big)\Big(x_i-\mu_l\Big)\Big(-1\Big)\Bigg)p(l\mid{x_i,\Theta^g}) \]
因為協方差矩陣Sigma為對稱陣，故\(\frac{1}{2}\Big(\Sigma_l^{-1}+(\Sigma_l^{-1})^T\Big)=\Sigma_l^{-1}\)，因此，上式繼續化簡為：
\[\sum_{l=1}^k\sum_{i=1}^N\Sigma_l^{-1}(x_i-\mu_l)p(l\mid{x_i,\Theta^g}) \]
令該式為0，得到：
\[\begin{aligned} \sum_{l=1}^k\Sigma_l^{-1}\sum_{i=1}^N{\mu_l}p(l\mid{x_i,\Theta^g})&=\sum_{l=1}^k\Sigma_l^{-1}\sum_{i=1}^N{x_i}p(l\mid{x_i,\Theta^g}) \\ \sum_{i=1}^N{\mu_l}p(l\mid{x_i,\Theta^g})&=\sum_{i=1}^N{x_i}p(l\mid{x_i,\Theta^g}) \\ \mu_l\sum_{i=1}^Np(l\mid{x_i,\Theta^g})&=\sum_{i=1}^N{x_i}p(l\mid{x_i,\Theta^g}) \\ \mu_l&=\frac{\sum_{i=1}^N{x_i}p(l\mid{x_i,\Theta^g})}{\sum_{i=1}^Np(l\mid{x_i,\Theta^g})} \\ \end{aligned}\]
估計參數\(\Sigma_l\)
\[\begin{aligned} &\quad\sum_{l=1}^k\sum_{i=1}^N\log\Big(p_l(x_i\mid{\mu_l,\Sigma_l})\Big)p(l\mid{x_i,\Theta^g}) \\ &= \sum_{l=1}^k\sum_{i=1}^N\Big(\log(2\pi^{-d/2})+\log|\Sigma_l|^{-1/2}-\frac{1}{2}(x_i-\mu_l)^T\Sigma_l^{-1}(x_i-\mu_l)\Big)p(l\mid{x_i,\Theta^g}) \\ \end{aligned}\]
忽略其中的常數項（因為常數項在求導之后為0），上式化簡為：
\[\begin{aligned} &\quad\sum_{l=1}^k\sum_{i=1}^N\Big(-\frac{1}{2}\log|\Sigma_l|-\frac{1}{2}(x_i-\mu_l)^T\Sigma_l^{-1}(x_i-\mu_l)\Big)p(l\mid{x_i,\Theta^g}) \\ &= \sum_{l=1}^k\sum_{i=1}^N\Big(\frac{1}{2}\log|\Sigma_l^{-1}|p(l\mid{x_i,\Theta^g})-\frac{1}{2}(x_i-\mu_l)^T\Sigma_l^{-1}(x_i-\mu_l)p(l\mid{x_i,\Theta^g})\Big) \\ &= \sum_{l=1}^k\Big(\frac{1}{2}\log|\Sigma_l^{-1}|\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})-\frac{1}{2}\Sigma_l^{-1}\sum_{i=1}^{N}(x_i-\mu_l)^T(I)(x_i-\mu_l)p(l\mid{x_i,\Theta^g})\Big) \\ &= \sum_{l=1}^k\Bigg(\frac{1}{2}\log|\Sigma_l^{-1}|\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})-\frac{1}{2}tr\Big(\Sigma_l^{-1}\sum_{i=1}^{N}(x_i-\mu_l)(x_i-\mu_l)^Tp(l\mid{x_i,\Theta^g})\Big)\Bigg) \\ \end{aligned}\]
考慮方程\(S(\mu_l,\Sigma_{l}^{-1})\)為上式中\(\sum_{l=1}^k\)內部的式子，即：
\[S(\mu_l,\Sigma_{l}^{-1})= \frac{1}{2}\log|\Sigma_l^{-1}|\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})-\frac{1}{2}tr\Big(\Sigma_l^{-1}\underbrace{\sum_{i=1}^{N}(x_i-\mu_l)(x_i-\mu_l)^Tp(l\mid{x_i,\Theta^g})}_{M_l}\Big) \]
\(S\)關於\(\Sigma_l^{-1}\)求導數，得到：
\[\begin{aligned} \frac{\partial{S(\mu_l,\Sigma_l^{-1})}}{\partial{\Sigma_l^{-1}}} &= \frac{1}{2}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})\Big(2\Sigma_l-diag(\Sigma_l)\Big)-\frac{1}{2}\frac{\partial{tr(\Sigma_l^{-1}M_l)}}{\partial\Sigma_l^{-1}} \\ &= \frac{1}{2}\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})\Big(2\Sigma_l-diag(\Sigma_l)\Big)-\frac{1}{2}\Big(2M_l-diag(M_l)\Big) \\ \end{aligned}\]
令該導數值為0，得到：
\[\begin{aligned} \sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})\Big(2\Sigma_l-diag(\Sigma_l)\Big) &= \Big(2M_l-diag(M_l)\Big) \\ 2\Big(\underbrace{\sum_{i=1}^{N}\Sigma_lp(l\mid{x_i,\Theta^g})-M_l}_{\mathcal{A}}\Big) &= diag\Big(\underbrace{\sum_{i=1}^{N}\Sigma_lp(l\mid{x_i,\Theta^g})-M_l}_{\mathcal{A}}\Big) \\ 2\mathcal{A} &= diag(\mathcal{A}) \end{aligned}\]
解得\(\mathcal{A}=0\)，即：
\[\begin{aligned} \sum_{i=1}^{N}\Sigma_lp(l\mid{x_i,\Theta^g}) &= M_l \\ \sum_{i=1}^{N}\Sigma_lp(l\mid{x_i,\Theta^g}) &= \sum_{i=1}^{N}(x_i-\mu_l)(x_i-\mu_l)^Tp(l\mid{x_i,\Theta^g}) \\ \Sigma_l\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g}) &= \sum_{i=1}^{N}(x_i-\mu_l)(x_i-\mu_l)^Tp(l\mid{x_i,\Theta^g}) \\ \Sigma_l &= \frac{\sum_{i=1}^{N}(x_i-\mu_l)(x_i-\mu_l)^Tp(l\mid{x_i,\Theta^g})}{\sum_{i=1}^{N}p(l\mid{x_i,\Theta^g})} \end{aligned}\]
