最近接觸了pLSA模型，該模型需要使用期望最大化（Expectation Maximization）算法求解。
本文簡述了以下內容：
為什么需要EM算法
EM算法的推導與流程
EM算法的收斂性定理
使用EM算法求解三硬幣模型
為什么需要EM算法
數理統計的基本問題就是根據樣本所提供的信息，對總體的分布或者分布的數字特征作出統計推斷。所謂總體，就是一個具有確定分布的隨機變量，來自總體的每一個iid樣本都是一個與總體有相同分布的隨機變量。
參數估計是指這樣一類問題——總體所服從的分布類型已知，但某些參數未知：設 $Y_1,...,Y_N$ 是來自總體 $\textbf Y$ 的iid樣本，記 $Y=(y_1,...,y_N)$ 是樣本觀測值，如果隨機變量 $Y_1,...,Y_N$ 是可觀測的，那么直接用極大似然估計法就可以估計參數 $\theta$ 。
但是，如果里面含有不可觀測的隱變量，使用MLE就沒那么容易。EM算法正是服務於求解帶有隱變量的參數估計問題。
EM算法的推導與流程
下面考慮帶有隱變量 $\boldsymbol Z$ (觀測值為 $z$ )的參數估計問題。將觀測數據（亦稱不完全數據）記為 $Y=(y_1,...,y_N)$ ，不可觀測數據記為 $Z=(z_1,...,z_N)$ ， $Y$ 、$Z$ 合在一起稱為完全數據。那么觀測數據的似然函數為$$l(\theta)=\prod_{j=1}^NP(y_j|\theta)=\prod_{j=1}^N\sum_zP(z|\theta)P(y_j|z,\theta)$$
其中求和號表示對 $z$ 的所有可能取值求和。另外，當使用頻率學派的方法估計參數時，參數 $\theta$ 並不是一個隨機變量，所以概率的記號里凡是 $\theta$ 出現在“條件”的地方（比如 $P(y_j|\theta)$ ），我覺得應該是用分號隔開（ $P(y_j;\theta)$ ）。這里為了和參考資料里描述的一致，就不改了。
為了省事，表述成這個形式：$$l(\theta)=P(Y|\theta)=\sum_zP(z|\theta)P(Y|z,\theta)$$
不然的話后面所有的推導都要套一個求和號（因為會取對數似然，求積變求和）。
對數似然：$$L(\theta)=\ln P(Y|\theta)=\ln \sum_zP(z|\theta)P(Y|z,\theta)$$
EM算法是一種迭代算法，通過迭代的方式求取目標函數 $L(\theta)=\ln P(Y|\theta)$ 的極大值。因此，希望每一步迭代之后的目標函數值會比上一步迭代結束之后的值大。設當前第 $n$ 次迭代后參數取值為 $\theta_n$ ，我們的目的是使 $L(\theta_{n+1})>L(\theta_n)$ 。那么考慮：
$$L(\theta)-L(\theta_n)=\ln (\sum_zP(z|\theta)P(Y|z,\theta))-\ln P(Y|\theta_n)$$
使用琴生不等式（Jensen inequality）：
$$\ln\sum_j\lambda_jy_j\geq \sum_j\lambda_j\log y_j,\quad \lambda_j\ge 0,\sum_j\lambda_j=1$$
因為 $\sum_zP(z|Y,\theta_n)=1$ ，所以 $L(\theta)-L(\theta_n)$ 的第一項有
$$\begin{aligned} \ln (\sum_zP(z|\theta)P(Y|z,\theta))&=\ln (\sum_zP(z|Y,\theta_n)\frac{P(z|\theta)P(Y|z,\theta)}{P(z|Y,\theta_n)})\\&\ge \sum_zP(z|Y,\theta_n)\ln \frac{P(z|\theta)P(Y|z,\theta)}{P(z|Y,\theta_n)}\end{aligned}$$
第二項有
$$ -\ln P(Y|\theta_n)=-\sum_zP(z|Y,\theta_n)\ln P(Y|\theta_n) $$
則 $L(\theta)-L(\theta_n)$ 的下界為
$$\begin{aligned} L(\theta)-L(\theta_n)&\ge\sum_zP(z|Y,\theta_n)\ln\frac{P(z|\theta)P(Y|z,\theta)}{P(z|Y,\theta_n)}-\sum_zP(z|Y,\theta_n)\ln P(Y|\theta_n)\\&=\sum_z[P(z|Y,\theta_n)\ln\frac{P(z|\theta)P(Y|z,\theta)}{P(z|Y,\theta_n)}-P(z|Y,\theta_n)\ln P(Y|\theta_n)]\\&=\sum_zP(z|Y,\theta_n)\ln\frac{P(z|\theta)P(Y|z,\theta)}{P(Y|\theta_n)P(z|Y,\theta_n)} \end{aligned}$$
定義一個函數 $l(\theta|\theta_n)$ ：
$$l(\theta|\theta_n)\triangleq L(\theta_n)+\sum_zP(z|Y,\theta_n)\ln\frac{P(z|\theta)P(Y|z,\theta)}{P(Y|\theta_n)P(z|Y,\theta_n)}$$
從而有 $L(\theta)\ge l(\theta|\theta_n)$ ，也就是說第 $n$ 次迭代結束后， $L(\theta)$ 的一個下界為 $l(\theta|\theta_n)$ 。另外，有等式 $L(\theta_n)=l(\theta_n|\theta_n)$ 成立。
我們的目的是使下一次迭代后得到的目標函數值能夠大於當前的值： $L(\theta_{n+1})>L(\theta_n)$ ，即 $L(\theta_{n+1})>l(\theta_n|\theta_n)$ 。
而在當前， $L(\theta)$ 的下界為 $l(\theta|\theta_n)$ ，因此，任何能讓 $l(\theta|\theta_n)$ 增大的 $\theta$ ，也可以讓 $L(\theta)$ 增大。
也就是說，能滿足 $l(\theta_{n+1}|\theta_n)>l(\theta_n|\theta_n)$ 的 $\theta_{n+1}$ ，一定更能滿足$L(\theta_{n+1})>l(\theta_n|\theta_n)=L(\theta_n)$ 。
通過下圖（來源：參考資料[1]，自己做了點注釋）可以解釋：
需要注意的是，下界的曲線當然是隨着迭代的進行而變化的：在第 $i$ 次迭代結束后，總是有不等式 $L(\theta)\ge l(\theta|\theta_i)$ 和等式 $L(\theta_i)=l(\theta_i|\theta_i)$ 成立。
換句話說，EM算法通過優化對數似然在當前的下界，來間接優化對數似然。
ok，那么現在問題就從找滿足 $L(\theta_{n+1})>L(\theta_n)$ 的 $\theta_{n+1}$ ，
轉變成了找滿足 $l(\theta_{n+1}|\theta_n)>l(\theta_n|\theta_n)$ 的 $\theta_{n+1}$ 。如何找到這樣一個 $\theta_{n+1}$ ？直接找 $l(\theta|\theta_n)$ 的最優解唄：
$$\theta_{n+1}=\arg\max_\theta l(\theta|\theta_n)$$
把 $l(\theta|\theta_n)$ 中的幾個與 $\theta$ 無關的量去掉，從而有
$$\begin{aligned} \theta_{n+1}&=\arg\max_\theta \sum_zP(z|Y,\theta_n)\ln [P(z|\theta)P(Y|z,\theta)]\\&=\arg\max_\theta \sum_zP(z|Y,\theta_n)\ln P(Y,z|\theta) \end{aligned}$$
回顧一下隨機變量的期望的表達式：
$$\mathbb E[\boldsymbol Z]=\sum_kP(\boldsymbol Z=z_k)z_k$$
$$\mathbb E[g(\boldsymbol Z)]=\sum_kP(\boldsymbol Z=z_k)g(z_k)$$
$$\mathbb E[\boldsymbol Z|\boldsymbol Y=y]=\sum_kP(\boldsymbol Z=z_k|\boldsymbol Y=y)z_k$$
所以：
$$\begin{aligned} \theta_{n+1}&=\arg\max_\theta \mathbb E_{\boldsymbol Z|\boldsymbol Y,\theta_n}[\ln P(Y,z|\theta)]\\&=\arg\max_\theta Q(\theta|\theta_n) \end{aligned}$$
上式定義了一個函數 $Q(\theta|\theta_n)$ ，稱為 $Q$ 函數。
上式完整表明了EM算法中的一步迭代中所需要的兩個步驟：E-step，求期望；M-step，求極大值。有了上面的鋪墊，下面介紹EM算法的流程：
輸入：觀測數據 $Y$ ，不可觀測數據 $Z$；
輸出：參數 $\theta$ ；
步驟：1. 給出參數初始化值 $\theta_0$ ；
2. E步：記 $\theta_n$ 為第 $n$ 次迭代后參數的估計值。在第 $n+1$ 次迭代的E步，求期望（ $Q$ 函數）
$$Q(\theta|\theta_n)=\mathbb E_{\boldsymbol Z|\boldsymbol Y,\theta_n}[\ln P(Y,z|\theta)]=\sum_zP(z|Y,\theta_n)\ln P(Y,z|\theta)$$
3. M步：求 $Q$ 函數的極大值點，來作為第 $n+1$ 次迭代所得到的參數估計值 $\theta_{n+1}$
$$\theta_{n+1}=\arg\max_\theta Q(\theta|\theta_n)$$
4. 重復上面兩步，直至達到停機條件。
EM算法的收斂性定理
定理1：觀測數據的似然函數 $P(Y|\theta)$ 通過EM算法得到的序列 $P(Y|\theta_n)(n=1,2,...)$ 單調遞增：$P(Y|\theta_{n+1})\ge P(Y|\theta_n)$ 。
定理2：(1) 如果 $P(Y|\theta)$ 有上界，則 $L(\theta_n)=\ln P(Y|\theta_n)$ 收斂到某一值 $L^*$ ；
(2) 在 $Q$ 函數與 $L(\theta)$ 滿足一定條件下，由EM算法得到的參數估計序列 $\theta_n$ 的收斂值 $\theta^*$ 是 $L(\theta)$ 的穩定點。
定理2中第二點的“條件”在大多數情況下都滿足。只能保證收斂到穩定點，不能保證收斂到極大值點，因此EM算法受初值的影響較大。
使用EM算法求解三硬幣模型
參考資料[2]給出了三硬幣模型的描述：
假設有三枚硬幣A、B、C，這些硬幣正面出現的概率分別是 $\pi$ 、$p$ 、$q$ 。進行如下擲硬幣試驗： 先擲A，如果A是正面則再擲B，如果A是反面則再擲C。對於B或C的結果，如果是正面則記為1，如果是反面則記為0。進行 $N$ 次獨立重復實驗，得到結果。現在只能觀測到結果，不能觀測到擲硬幣的過程，估計模型參數 $\theta=(\pi,p,q)$ 。
在這個問題中，實驗結果是可觀測數據 $Y=(y_1,...,y_N)$ ，硬幣A的結果是不可觀測數據 $Z=(z_1,...,z_N)$ 且 $z$ 只有兩種可能取值1和0。
對於第 $j$ 次試驗，
$$\begin{aligned} P(y_j|\theta)&= \sum_zP(y_j,z|\theta)\\&=\sum_zP(z|\theta)P(y_j|z,\theta)\\&=P(z=1|\theta)P(y_j|z=1,\theta)+P(z=0|\theta)P(y_j|z=0,\theta) \\&=\begin{cases} \pi p+(1-\pi )q, & \text{if }y_j=1;\\ \pi (1-p)+(1-\pi )(1-q), & \text{if }y_j=0. \end{cases} \\&=\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi )q^{y_j}(1-q)^{1-y_j} \end{aligned}$$
所以有
$$P(Y|\theta)=\prod_{j=1}^NP(y_j|\theta)=\prod_{j=1}^N (\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi )q^{y_j}(1-q)^{1-y_j})$$
1. E-step，求期望（Q函數）：
$$\begin{aligned} Q(\theta|\theta_n)&=\sum_zP(z|Y,\theta_n)\ln P(Y,z|\theta)\\&=\sum_{j=1}^N\{ \sum_zP(z|y_j,\theta_n)\ln P(y_j,z|\theta) \}\\&=\sum_{j=1}^N\{ P(z=1|y_j,\theta_n)\ln P(y_j,z=1|\theta) + P(z=0|y_j,\theta_n)\ln P(y_j,z=0|\theta) \} \end{aligned}$$
先求 $P(z|y_j,\theta_n)$ ，
$$P(z|y_j,\theta_n)=\begin{cases} \frac{\pi_n p_n^{y_j}(1-p_n)^{1-y_j}}{\pi_n p_n^{y_j}(1-p_n)^{1-y_j}+(1-\pi_n )q_n^{y_j}(1-q_n)^{1-y_j}}=\mu_{j,n} & \text{if }z=1; \\1-\mu_{j,n} & \text{if }z=0. \end{cases}$$
再求 $P(y_j,z|\theta)=P(z|\theta)P(y_j|z,\theta)$ ，
$$P(y_j,z|\theta)=\begin{cases} \pi p^{y_j}(1-p)^{1-y_j} &\text{if }z=1;\\ (1-\pi )q^{y_j}(1-q)^{1-y_j} &\text{if }z=0. \end{cases}$$
因此，$Q$ 函數表達式為：
$$Q(\theta|\theta_n)=\sum_{j=1}^N\{\mu_{j,n}\ln [\pi p^{y_j}(1-p)^{1-y_j}]+(1-\mu_{j,n})\ln [(1-\pi )q^{y_j}(1-q)^{1-y_j}] \}$$
2. M-step，求 $Q$ 函數的極大值：
令 $Q$ 函數對參數求導數，並等於零。
$$\begin{aligned} \frac{\partial Q(\theta|\theta_n)}{\partial \pi}&=\sum_{j=1}^N\{\frac{\mu_{j,n}\ln [\pi p^{y_j}(1-p)^{1-y_j}]+(1-\mu_{j,n})\ln [(1-\pi )q^{y_j}(1-q)^{1-y_j}] }{\partial \pi}\}\\&=\sum_{j=1}^N\{ \mu_{j,n}\frac{p^{y_j}(1-p)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j}}+(1-\mu_{j,n})\frac{-q^{y_j}(1-q)^{1-y_j}}{(1-\pi )q^{y_j}(1-q)^{1-y_j}} \}\\&=\sum_{j=1}^N\{ \frac{\mu_{j,n}-\pi }{\pi (1-\pi)}\}\\&=\frac{(\sum_{j=1}^N\mu_{j,n})-n\pi }{\pi (1-\pi)} \end{aligned}$$
$$\begin{aligned}\frac{\partial Q(\theta_|\theta_n)}{\partial \pi}=0 &\implies \pi =\frac 1n\sum_{j=1}^N\mu_{j,n}\\ \therefore \pi_{n+1}&=\frac 1n\sum_{j=1}^N\mu_{j,n} \end{aligned}$$
$$\begin{aligned} \frac{\partial Q(\theta|\theta_n)}{\partial p}&=\sum_{j=1}^N\{\frac{\mu_{j,n}\ln [\pi p^{y_j}(1-p)^{1-y_j}]+(1-\mu_{j,n})\ln [(1-\pi )q^{y_j}(1-q)^{1-y_j}] }{\partial p}\}\\&=\sum_{j=1}^N\{ \mu_{j,n}\frac{\pi (y_jp^{y_j-1}(1-p)^{1-y_j}+p^{y_j}(-1)(1-y_j)(1-p)^{1-y_j-1})}{\pi p^{y_j}(1-p)^{1-y_j}}+0 \}\\&=\sum_{j=1}^N\{ \frac{\mu_{j,n}(y_j-p) }{p(1-p)}\}\\&=\frac{(\sum_{j=1}^N\mu_{j,n}y_j)-(p\sum_{j=1}^N\mu_{j,n}) }{p(1-p)} \end{aligned}$$
$$\begin{aligned}\frac{\partial Q(\theta_|\theta_n)}{\partial p}=0 &\implies p =\frac{\sum_{j=1}^N\mu_{j,n}y_j}{\sum_{j=1}^N\mu_{j,n}}\\ \therefore p_{n+1}&=\frac{\sum_{j=1}^N\mu_{j,n}y_j}{\sum_{j=1}^N\mu_{j,n}}\\q_{n+1}&=\frac{\sum_{j=1}^N(1-\mu_{j,n})y_j}{\sum_{j=1}^N(1-\mu_{j,n})} \end{aligned}$$
既然已經得到了三個參數的迭代式，便可給定初值，迭代求解了。
參考資料：
[1] The Expectation Maximization Algorithm: A short tutorial - Sean Borman
[2] 《統計學習方法》，李航
[3] 梯度下降和EM算法：系出同源，一脈相承 這是最近掃到的一篇，還沒有看，感覺不錯，有時間看一下～