EM算法详解

EM算法

EM算法是含有隐变量的概率模型参数的极大似然估计法。

用$Y$表示观测变量的数据，$Z$表示隐变量的数据，$\theta$表示要估计的参数，$Y$和$Z$连在一起称为完全数据，观测数据$Y$称为不完全数据，假设$Y$的概率分布是$P(Y\mid \theta )$，那么不完全数据$Y$的对数似然函数是$\log P(Y\mid \theta )$，假设$Y$和$Z$的联合概率分布是$P(Y,Z\mid \theta )$，那么完全数据的对数似然函数是$\log P(Y,Z\mid \theta )$。

对含有隐变量的概率模型，目标是极大化观测数据$Y$对于模型参数$\theta$的对数似然函数，即极大化
$$L(\theta )=\log P(Y\mid \theta )=\log (\sum _{z}P(Y\mid Z,\theta )P(Z\mid \theta ))$$
这个式子的困难在于存在未观测数据并且包含和的对数，EM算法通过迭代逐步近似极大化$L(\theta )$，假设在第$i$次迭代后$\theta$的估计值是${\theta }^{(i)}$，我们希望新的估计值$\theta$能使$L(\theta )$增加，即$L(\theta )>L({\theta }^{(i)})$，并逐步达到极大值，为此考虑两者的差
$$L(\theta )-L({\theta }_{(i)})=\log (\sum _{z}P(Y\mid Z,\theta )P(Z\mid \theta ))-\log P(Y\mid {\theta }^{(i)})$$
利用$Jensen$不等式
$$\begin{array}{rl}L(\theta )-L({\theta }_{(i)})& =\log (\sum _{z}P(Z\mid Y,{\theta }^{(i)})\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})})-\log P(Y\mid {\theta }^{(i)})\\ & \ge \sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})})-\log P(Y\mid {\theta }^{(i)})\\ & =\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})})-\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log P(Y\mid {\theta }^{(i)})\\ & =\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})P(Y\mid {\theta }^{(i)})})\end{array}$$
令
$$B(\theta ,{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})P(Y\mid {\theta }^{(i)})})$$

即
$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$
并且$L({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})$，因此任何能使$B(\theta ,{\theta }^{(i)})$增大的$\theta$也一定能使$L(\theta )$增大，为使$L(\theta )$增长尽可能的大，应选择${\theta }^{(i+1)}$使$B(\theta ,{\theta }^{(i)})$达到极大，即
$$\begin{array}{rl}{\theta }^{(i+1)}& =\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})\\ & =\arg \underset{\theta }{\max }L({\theta }^{(i)})+\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P(Z\mid Y,{\theta }^{(i)})P(Y\mid {\theta }^{(i)})})\\ & =\arg \underset{\theta }{\max }\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log (P(Y\mid Z,\theta )P(Z\mid \theta )\\ & =\arg \underset{\theta }{\max }\sum _{z}P(Z\mid Y,{\theta }^{(i)})\log P(Y,Z\mid \theta )\\ & =\arg \underset{\theta }{\max }{E}_z[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}]\end{array}$$
令
$$Q(\theta ,{\theta }^{(i)})=E_z[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}]$$
即完全数据的对数似然$log(Y,Z\mid \theta )$关于在给定观测数据和当前参数${\theta }^{(i)}$下对未观测数据$Z$的期望，因此
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
EM算法：

输入：观测数据变量$Y$，隐变量数据$Z$，联合分布$P(Y,Z\mid \theta )$，条件分布$P(Z\mid Y,\theta )$

输出：模型参数$\theta$

(1) 选择参数的初始值${\theta }^{(0)}$，开始迭代

(2) E步：计算$Q(\theta ,{\theta }^{(i)})={\sum }_{z}P(Z\mid Y,{\theta }^{(i)})\log P(Y,Z\mid \theta )$

(3) M步：计算${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$

(4) 重复执行第(2)步和第(3)步，直至收敛

高斯混合模型

高斯混合分布是具有如下形式的概率分布
$$P(x\mid \theta )=\sum _{k=1}^K{\alpha }_k\varphi (x\mid {\theta }_k)$$
其中${\alpha }_k$是系数，${\sum }_{k=1}^K{\alpha }_k=1$，${\alpha }_k\ge 0$，$\varphi (y\mid {\theta }_k)$是高斯概率密度函数，${\theta }_k=({\mu }_k,{\sigma }_k^2)$
$$\varphi (x\mid {\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp (-\frac{(x-{\mu }_k{)}^2}{2{\sigma }_k^2})$$
假设样本的生成过程由高斯混合分布给出：首先根据${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k$选择一个高斯混合成分，然后根据被选择的高斯混合成分生成观测数据。这是观测数据是已知的，观测数据来自哪个高斯分布是未知的，以隐变量${\gamma }_{jk}$表示，其定义如下：
$$\begin{array}{rl}{\gamma }_{jk}=& \{\begin{array}{l}1,第j个观测来自第k个分模型\\ 0,其他\end{array}\\ & \\ & j=1,2,\dots ,N;k=1,2,\dots ,k\end{array}$$
那么完全数据是
$$(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jk}),j=1,2,\dots ,N$$
于是完全数据的似然函数
$$\begin{array}{rl}P(\gamma ,y\mid \theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jk}\mid \theta )\\ & =\prod _{j=1}^N\prod _{k=1}^K[{\alpha }_k\varphi (x_j\mid {\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (x_j\mid {\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp (-\frac{(x_j-{\mu }_k{)}^2}{2{\sigma }_k^2}){]}^{{\gamma }_{jk}}\end{array}$$
其中$n_k={\sum }_{j=1}^N{\gamma }_{jk}$，${\sum }_{k=1}^K{n}_k=N$，那么完全数据的对数似然是
$$\log P(\gamma ,y\mid \theta )=\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{k=1}^K\sum _{j=1}^N[{\gamma }_{jk}(\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_k-\frac{(x_j-{\mu }_k{)}^2}{2{\sigma }_k^2})]$$
需要极大化的$Q$函数是
Q(θ,θ(i))=Eγ[log⁡P(γ,y∣θ)∣y,θ(i)]=Eγ{∑k=1Knklog⁡αk+∑k=1K∑j=1N[γjk(log⁡12π−log⁡σk−(xj−μk)22σk2)]}=∑k=1K{∑j=1N(Eγjk)log⁡αk+∑j=1N(Eγjk)[log⁡12π−log⁡σk−(xj−μk)22σk2]} \begin{aligned} Q(\theta,\theta^{(i)}) &=E_\gamma[\log P(\gamma,y\mid\theta)\mid y,\theta^{(i)}]\\ &=E_\gamma\{\sum_{k=1}^{K}n_k\log\alpha_k+\sum_{k=1}^{K}\sum_{j=1}^{N}[\gamma_{jk}(\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(x_j-\mu_k)^2}{2\sigma_k^2})]\}\\ &=\sum_{k=1}^{K}\{\sum_{j=1}^{N}(E\gamma_{jk})\log\alpha_k+\sum_{j=1}^{N}(E\gamma_{jk})[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(x_j-\mu_k)^2}{2\sigma_k^2}]\} \end{aligned} Q(θ,θ(i))​=Eγ​[logP(γ,y∣θ)∣y,θ(i)]=Eγ​{k=1∑K​nk​logαk​+k=1∑K​j=1∑N​[γjk​(log2π​1​−logσk​−2σk2​(xj​−μk​)2​)]}=k=1∑K​{j=1∑N​(Eγjk​)logαk​+j=1∑N​(Eγjk​)[log2π​1​−logσk​−2σk2​(xj​−μk​)2​]}​
这里需要计算$E{\gamma }_{jk}=E({\gamma }_{jk}\mid y,\theta )$
$$\begin{array}{rl}E({\gamma }_{jk}\mid y,\theta )& =P({\gamma }_{jk}=1\mid y,\theta )\\ & =\frac{P({\gamma }_{jk}=1,y_j\mid \theta )}{P(y_j\mid \theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j\mid \theta )}{{\sum }_{k=1}^{K}P({\gamma }_{jk}=1,y_j\mid \theta )}\\ & =\frac{P(y_j\mid {\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1\mid \theta )}{{\sum }_{k=1}^{K}P(y_j\mid {\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1\mid \theta )}\\ & =\frac{{\alpha }_k\varphi (y_j\mid {\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j\mid {\theta }_k)}\end{array}$$
$E({\gamma }_{jk}\mid y,\theta )$表示当前参数下第$j$个观测数据来自第$k$个混合成分的概率，记为${\hat{\gamma }}_{jk}$。综上所述
Q(θ,θ(i))=∑k=1K{nklog⁡αk+∑j=1Nγ^jk[log⁡12π−log⁡σk−(xj−μk)22σk2]} Q(\theta,\theta^{(i)})=\sum_{k=1}^{K}\{n_k\log\alpha_k+\sum_{j=1}^{N}\hat{\gamma}_{jk}[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(x_j-\mu_k)^2}{2\sigma_k^2}]\} Q(θ,θ(i))=k=1∑K​{nk​logαk​+j=1∑N​γ^​jk​[log2π​1​−logσk​−2σk2​(xj​−μk​)2​]}
其中$n_k={\sum }_{j=1}^N{\hat{\gamma }}_{jk}$，接下来需对$Q(\theta ,{\theta }^{(i)})$求极大，需要求$Q(\theta ,{\theta }^{(i)})$对每个参数的偏导。

$Q(\theta ,{\theta }^{(i)})$对${\mu }_k$的偏导：
$$\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\mu }_k}=\sum _{j=1}^N\frac{{\gamma }_{jk}(x_j-{\mu }_k)}{{\sigma }_k^2}$$
所以
$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\gamma }_{jk}{x}_j}{{\sum }_{j=1}^N{\gamma }_{jk}}$$
$Q(\theta ,{\theta }^{(i)})$对${\sigma }_k^2$的偏导：
$$\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\sigma }_k^2}=-\frac{1}{2}\sum _{j=1}^N{\hat{\gamma }}_{jk}(\frac{1}{{\sigma }_k^2}-\frac{(x_j-{\mu }_k{)}^2}{{\sigma }_k^4})$$
所以
$${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}^{jk}(x_j-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}$$
$Q(\theta ,{\theta }^{(i)})$对${\alpha }_k$的偏导：由于存在约束条件${\sum }_{k=1}^K{\alpha }_k=1$，所以考虑$Q(\theta ,{\theta }^{(i)})$的拉格朗日函数
$$L(\theta ,{\theta }^{(i)})=Q(\theta ,{\theta }^{(i)})+\lambda (\sum _{k=1}^K{\alpha }_k-1)$$
求偏导
$$\frac{\partial L(\theta ,{\theta }^{(i)})}{\partial {\alpha }_k}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{{\alpha }_k}+\lambda$$
令偏导等于零，即
$$\sum _{j=1}^N{\hat{\gamma }}_{jk}+\lambda {\alpha }_k=0$$
为求解$\lambda$，对所有分模型求和得
$$\sum _{k=1}^K\sum _{j=1}^N{\hat{\gamma }}_{jk}+\lambda \sum _{k=1}^K{\alpha }_k=0$$
解得$\lambda =-N$，所以
$${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N}$$

多元高斯分布与此类似，更新公式几乎是一模一样的，但计算却更为复杂，涉及到复杂的求导。

根据高斯混合模型的计算结果，可以得到一种聚类方法，即根据$\gamma$的取值判断每个样本属于哪一类，对于样本$x_j$，其所属的类别是
$$\arg \underset{k}{\max }{\gamma }_{jk}$$