EM算法与GMM高斯混合模型

参考知乎上的这篇文章：EM算法详解 by Microstrong，感谢大佬写得非常详细。

EM算法

预备知识：MLE，Jensen不等式。（不展开介绍）
EM算法可以理解为，从MLE出发，为了解决包含隐变量（或缺失数据）的参数估计问题而设计的一种迭代算法，该算法的迭代步骤分为E步和M步，E步在有了初始化参数值的情况下，需要求解隐变量的条件分布，进而求得一个包含隐变量的函数的期望（目标函数），M步需要最大化目标函数，求得更新后的参数值。

数学推导

MLE中的对数似然函数，可以简单写为，
$$\ln L=\sum _i\ln P(x=x_i)$$
或者考虑隐变量 $z$（或缺失数据），以及模型参数时，公式写为，
$$\ln L=\sum _i\ln P(x=x_i,z;\theta )$$
这里用了一个很神奇的思路，就是利用不等式和等号成立的条件进行恒等变换，以方便设计迭代算法。
这里就用了Jensen不等式，以上凸函数（凹函数）为例（$g(x)=\ln x$)，
$$g(E(X))\ge E(g(X))$$
假设$X$服从一个分布$P(X)$，上式可以写为，
$$g(\sum _i{x}_{i}p(x=x_i))\ge \sum _{i}g(x_i)p(x=x_i)$$
或连续情况下，
$$g(\int xf(x)dx)\ge \int g(x_i)f(x)dx$$
考虑对前面的对数似然函数用Jensen不等式，假设$x=x_i$条件下，$z$服从某个分布$Q_i(z)$，可以写为，
ln ⁡ L = ∑ i ln ⁡ P ( x = x i , z ; θ ) = ∑ i ln ⁡ { ∑ z Q i ( z ) P ( x = x i , z ; θ ) Q i ( z ) } ≥ ∑ i ∑ z Q i ( z ) ln ⁡ P ( x = x i , z ; θ ) Q i ( z ) \begin{aligned} \ln{L}&=\sum_i{\ln{P(x=x_i,z;\theta)}} \\ &= \sum_i{\ln\{\sum_zQ_i(z)\frac{{P(x=x_i,z;\theta)}}{Q_i(z)}\}} \\ &\ge \sum_i{\sum_zQ_i(z)\ln\frac{{P(x=x_i,z;\theta)}}{Q_i(z)}} \\ \end{aligned} lnL​=i∑​lnP(x=xi​,z;θ)=i∑​ln{z∑​Qi​(z)Qi​(z)P(x=xi​,z;θ)​}≥i∑​z∑​Qi​(z)lnQi​(z)P(x=xi​,z;θ)​​
这里不等式等号成立的条件为，严格凸函数上取的坐标点都是同一点，也就是
$$\frac{P(x=x_i,z;\theta )}{Q_i(z)}=Const$$
再根据分布满足的条件，
$$\sum _z{Q}_i(z)=1$$
联立求得，
$$Q_i(z)=\frac{P(x=x_i,z;\theta )}{{\sum }_{z}P(x=x_i,z;\theta )}$$
若$z$的取值连续，也可用上面这个公式表示，而这个式子就是贝叶斯公式（后验概率），即，
$$Q_i(z)=P(z|x_i;\theta )$$

E步和M步

E步首先需要求出$P(z|x_i;\theta )$，至于具体怎么求，就是由模型决定的了，用贝叶斯公式总能算出来。注意，因为求这个后验概率时，用的是上一步（或初始化）得到的$\theta$，这里可以记为，
$$P(z|x_i;{\theta }_t)$$
有了这个之后，就可以写出目标函数，也就是MLE的对数似然函数再用Jensen不等式放缩得到的目标函数，
$$l(\theta ,{\theta }_t)=\sum _i\sum _{z}P(z|x_i;{\theta }_t)\ln \frac{P(x=x_i,z;\theta )}{P(z|x_i;{\theta }_t)}$$
之后进入M步，最大化这个目标函数以求得更新后的$\theta$，
$${\theta }_{t+1}=\underset{\theta }{\min }l(\theta ,{\theta }_t)$$
如果把E步和M步合起来，我们可以说EM算法就是下面这个更新公式。
$${\theta }_{t+1}=\underset{\theta }{\min }\sum _i\sum _{z}P(z|x_i;{\theta }_t)\ln \frac{P(x=x_i,z;\theta )}{P(z|x_i;{\theta }_t)}$$

GMM（Gaussian Mixture Model）

GMM是一种聚类算法，认为所有样本服从若干个高斯分布（正态分布）的叠加，模型为，
p ( x = x i ) = ∑ j p ( x = x i , μ = μ j ) = ∑ j p ( μ = μ j ) p ( x = x i ∣ μ = μ j ) = ∑ j p ( μ = μ j ) 1 ( 2 π ) p / 2 ∣ Σ j ∣ 1 / 2 e x p { − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) } \begin{aligned} p(x=x_i)&=\sum_j{p(x=x_i,\mu=\mu_j)} \\ &= \sum_j{p(\mu=\mu_j)p(x=x_i|\mu=\mu_j)} \\ &= \sum_j{p(\mu=\mu_j)\frac{1}{{(2\pi)}^{p/2}{|\Sigma_j|}^{1/2}}exp\{-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}\}} \\ \end{aligned} p(x=xi​)​=j∑​p(x=xi​,μ=μj​)=j∑​p(μ=μj​)p(x=xi​∣μ=μj​)=j∑​p(μ=μj​)(2π)p/2∣Σj​∣1/21​exp{−21​(xi​−μj​)TΣj−1​(xi​−μj​)}​
注意这个模型中的参数除了${\mu }_j$、${\Sigma }_j$外，还有$p(\mu ={\mu }_j)$，这个其实是一个定义每个分布权重或高度的，和为1的系数，也可以写作$w_j$。

而聚类问题其实就是缺失标签值的分类问题，因此标签值 $\mu$ 其实就相当于一个离散的隐变量，可以用EM算法解决。
此外，我们可以定义，
$$z_{ij}=\{\begin{array}{ll}1& \text{if }{x}_i\in {\mu }_j\\ 0& \text{if }{x}_i\notin {\mu }_j\end{array}$$
对$z_{ij}$求一下期望，也就是，
$$E(z_{ij})=P(x_i\in {\mu }_j)=p(\mu ={\mu }_j|x=x_i)$$
发现，$E(z_{ij})$就是把$\mu$作为隐变量时，要求的隐变量的后验概率。

E步

在初始化参数之后，首先求隐变量的条件分布，也就是$E(z_{i\cdot })$或$p(\mu |x=x_i)$，
E ( z i j ) = p ( μ = μ j ∣ x = x i ) = p ( μ = μ j ) p ( x = x i ∣ μ = μ j ) ∑ s p ( μ = μ s ) p ( x = x i ∣ μ = μ s ) = p ( μ = μ j ) 1 ( 2 π ) p / 2 ∣ Σ j ∣ 1 / 2 e x p { − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) } ∑ s p ( μ = μ s ) 1 ( 2 π ) p / 2 ∣ Σ s ∣ 1 / 2 e x p { − 1 2 ( x i − μ s ) T Σ s − 1 ( x i − μ s ) } \begin{aligned} E(z_{ij})&=p(\mu=\mu_j|x=x_i) \\ &= \frac{p(\mu=\mu_j)p(x=x_i|\mu=\mu_j)}{\sum_s{p(\mu=\mu_s)p(x=x_i|\mu=\mu_s)}} \\ &= \frac{p(\mu=\mu_j)\frac{1}{{(2\pi)}^{p/2}{|\Sigma_j|}^{1/2}}exp\{-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}\}}{\sum_s{p(\mu=\mu_s)\frac{1}{{(2\pi)}^{p/2}{|\Sigma_s|}^{1/2}}exp\{-\frac{1}{2}{(x_i-\mu_s)}^T\Sigma_s^{-1}{(x_i-\mu_s)}\}}} \\ \end{aligned} E(zij​)​=p(μ=μj​∣x=xi​)=∑s​p(μ=μs​)p(x=xi​∣μ=μs​)p(μ=μj​)p(x=xi​∣μ=μj​)​=∑s​p(μ=μs​)(2π)p/2∣Σs​∣1/21​exp{−21​(xi​−μs​)TΣs−1​(xi​−μs​)}p(μ=μj​)(2π)p/2∣Σj​∣1/21​exp{−21​(xi​−μj​)TΣj−1​(xi​−μj​)}​​
这里$E(z_{ij})$求出后就是个常数了，代入EM算法的目标函数中，
l ( μ , Σ , w ) = ∑ i ∑ j E ( z i j ) ln ⁡ P ( x = x i , μ = μ j ; μ , Σ , w ) E ( z i j ) = ∑ i ∑ j E ( z i j ) ln ⁡ w j 1 ( 2 π ) p / 2 ∣ Σ j ∣ 1 / 2 e x p { − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) } E ( z i j ) \begin{aligned} l(\mu,\Sigma,w)&=\sum_i\sum_jE(z_{ij})\ln\frac{{P(x=x_i,\mu=\mu_j;\mu,\Sigma,w)}}{E(z_{ij})} \\ &=\sum_i\sum_jE(z_{ij})\ln\frac{w_j\frac{1}{{(2\pi)}^{p/2}{|\Sigma_j|}^{1/2}}exp\{-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}\}}{E(z_{ij})} \end{aligned} l(μ,Σ,w)​=i∑​j∑​E(zij​)lnE(zij​)P(x=xi​,μ=μj​;μ,Σ,w)​=i∑​j∑​E(zij​)lnE(zij​)wj​(2π)p/2∣Σj​∣1/21​exp{−21​(xi​−μj​)TΣj−1​(xi​−μj​)}​​

M步

把对数部分展开化简，
l ( μ , Σ , w ) = ∑ i ∑ j E ( z i j ) { ln ⁡ w j − p 2 ln ⁡ 2 π − 1 2 ln ⁡ ∣ Σ j ∣ − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) − ln ⁡ E ( z i j ) } l(\mu,\Sigma,w)=\sum_i\sum_jE(z_{ij})\{\ln{w_j}-\frac{p}{2}\ln{2\pi}-\frac{1}{2}\ln{|\Sigma_j|}-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}-\ln{E(z_{ij})}\} l(μ,Σ,w)=i∑​j∑​E(zij​){lnwj​−2p​ln2π−21​ln∣Σj​∣−21​(xi​−μj​)TΣj−1​(xi​−μj​)−lnE(zij​)}
去除常数项后，
l ∗ ( μ , Σ , w ) = ∑ i ∑ j E ( z i j ) { ln ⁡ w j − 1 2 ln ⁡ ∣ Σ j ∣ − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) } l^*(\mu,\Sigma,w)=\sum_i\sum_jE(z_{ij})\{\ln{w_j}-\frac{1}{2}\ln{|\Sigma_j|}-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}\} l∗(μ,Σ,w)=i∑​j∑​E(zij​){lnwj​−21​ln∣Σj​∣−21​(xi​−μj​)TΣj−1​(xi​−μj​)}
注意这里有一个约束条件是，
$$\sum _j{w}_j=1$$
我们可以把包含$w_j$的部分单独拿出来研究，
$$\begin{array}{rl}\underset{w}{\max }& \sum _j[\ln w_j\sum _{i}E(z_{ij})]\\ s.t.& \sum _j{w}_j=1\end{array}$$
可以直接写出拉格朗日函数求解，
$$L_1(w_j,\alpha )=\sum _j[\ln w_j\sum _{i}E(z_{ij})]-\alpha (\sum _j{w}_j-1)$$
对$w_j$、$\alpha$分别求导，
$$\begin{array}{rl}\frac{\partial L_1}{\partial w_j}& =\frac{{\sum }_{i}E(z_{ij})}{w_j}-\alpha =0\\ \frac{\partial L_1}{\partial \alpha }& =\sum _j{w}_j-1=0\end{array}$$
又根据，
$$\sum _{j}E(z_{ij})=\sum _{j}p(\mu ={\mu }_j|x=x_i)=1$$
很容易求得，（式子中的$n$为样本数量）
$$\begin{gathered} \alpha =\sum _j\sum _{i}E(z_{ij})=n \\ {w}_j=\frac{1}{\alpha }\sum _{j}E(z_{ij})=\frac{1}{n}\sum _{j}E(z_{ij}) \end{gathered}$$
再看目标函数$l^{\ast }(\cdot )$去除$w_j$后的剩余部分为，
L 2 ( μ j , Σ j ) = ∑ i ∑ j E ( z i j ) { − 1 2 ln ⁡ ∣ Σ j ∣ − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) } L_2(\mu_j,\Sigma_j)=\sum_i\sum_jE(z_{ij})\{-\frac{1}{2}\ln{|\Sigma_j|}-\frac{1}{2}{(x_i-\mu_j)}^T\Sigma_j^{-1}{(x_i-\mu_j)}\} L2​(μj​,Σj​)=i∑​j∑​E(zij​){−21​ln∣Σj​∣−21​(xi​−μj​)TΣj−1​(xi​−μj​)}
对${\mu }_j$求导，
$$\frac{\partial L_2}{\partial {\mu }_j}=\sum _{i}E(z_{ij}){\Sigma }_j^{-1}(x_i-{\mu }_j)=0$$
很容易求得，
$${\mu }_j=\frac{{\sum }_{i}E(z_{ij})x_i}{{\sum }_{i}E(z_{ij})}$$
对${\Sigma }_j$求导（这里不会可以见我写的多元正态分布的最大似然估计），
∂ L 2 ∂ Σ j = ∑ i E ( z i j ) { − Σ j ∗ 2 ∣ Σ j ∣ + 1 2 Σ j − 2 ( x i − μ j ) ( x i − μ j ) T } = 0 \frac{\partial L_2}{\partial \Sigma_j}=\sum_i E(z_{ij})\{-\frac{\Sigma_j^*}{2|\Sigma_j|}+\frac{1}{2}\Sigma_j^{-2}{(x_i-\mu_j){(x_i-\mu_j)}^T}\}=0 ∂Σj​∂L2​​=i∑​E(zij​){−2∣Σj​∣Σj∗​​+21​Σj−2​(xi​−μj​)(xi​−μj​)T}=0
化简，
$${\Sigma }_j^{-1}\sum _{i}E(z_{ij})={\Sigma }_j^{-2}\sum _{i}E(z_{ij})(x_i-{\mu }_j){(x_i-{\mu }_j)}^T$$
得到，
$${\Sigma }_j=\frac{{\sum }_{i}E(z_{ij})(x_i-{\mu }_j){(x_i-{\mu }_j)}^T}{{\sum }_{i}E(z_{ij})}$$

这样，我们就得到了GMM用EM算法求解时的参数更新公式，最后整理一遍，如下，
$$w_j^{t+1}=\frac{1}{n}\sum _j{E}^t(z_{ij})$$$${\mu }_j^{t+1}=\frac{{\sum }_i{E}^t(z_{ij})x_i}{{\sum }_i{E}^t(z_{ij})}$$$${\Sigma }_j^{t+1}=\frac{{\sum }_i{E}^t(z_{ij})(x_i-{\mu }_j^{t+1}){(x_i-{\mu }_j^{t+1})}^T}{{\sum }_i{E}^t(z_{ij})}$$
其中，$E^t(z_{ij})$由E步，根据迭代前的参数值$w_j^t$、${\mu }_j^t$、${\Sigma }_j^t$求得。