EM算法 (第十五周周报12.6-12.12)

一.引例

1.引例1¹

设四种实验结果发生的概率依次为：$\frac{1}{2}-\frac{\theta }{4},\frac{1}{4}-\frac{\theta }{4},\frac{1}{4}+\frac{\theta }{4},\frac{\theta }{4}$;
发生的次数依次为：$y_1,y_2,y_3,y_4$次，求$\hat{\theta }$.

• 最大似然估计：
  $$\begin{array}{rl}L(\theta )& =(\frac{1}{2}-\frac{\theta }{4}{)}^{y_1}(\frac{1}{4}-\frac{\theta }{4}{)}^{y_2}(\frac{1}{4}+\frac{\theta }{4}{)}^{y_3}(\frac{\theta }{4}{)}^{y_4}\\ \ln L(\theta )& =y_1\ln (\frac{1}{2}-\frac{\theta }{4})+y_2\ln (\frac{1}{4}-\frac{\theta }{4})+y_3\ln (\frac{1}{4}+\frac{\theta }{4})+y_4\ln (\frac{\theta }{4})\\ \frac{\mathrm{d}\ln L(\theta )}{\mathrm{d}\theta }& =-\frac{y_1}{2-\theta }-\frac{y_2}{1-\theta }-\frac{y_3}{1+\theta }-\frac{y_4}{\theta }=0\end{array}$$
  由于导数是关于$\theta$的一元三次方程，求解困难，当然也可以用数值方法求解，但是如果涉及的实验结果不止4个，就显得麻烦，所以为了使方法更加通用化，则使用EM算法来求解.
• EM算法
  1）把事件拆分为两个和事件：
  假设发生结果概率为$\frac{1}{2}-\frac{\theta }{4}$的事件拆分为发生概率为$\frac{1}{4}-\frac{\theta }{4}$和$\frac{1}{4}$的两个事件，发生次数分为$z_1,y_1-z_1;$
  假设发生结果概率为$\frac{1}{4}+\frac{\theta }{4}$的事件拆分为发生概率为$\frac{\theta }{4}$和$\frac{1}{4}$的两个事件，发生次数分为$z_2,y_3-z_2;$
  2）由最大似然估计：
  $$\begin{array}{rl}L(\theta )& =(\frac{1}{4}{)}^{y_1-z_1}(\frac{1}{4}-\frac{\theta }{4}{)}^{y_2+z_1}(\frac{1}{4}{)}^{y_3-z_2}(\frac{\theta }{4}{)}^{y_4+z_2}\\ \ln L(\theta )& =(y_1-z_1)\ln \frac{1}{4}+(y_2+z_1)\ln (\frac{1}{4}-\frac{\theta }{4})+(y_3-z_2)\ln \frac{1}{4}+(y_4+z_2)\ln \frac{\theta }{4}\\ \frac{\mathrm{d}\ln L(\theta )}{\mathrm{d}\theta }& =-\frac{y_2+z_1}{1-\theta }+\frac{y_4+z_2}{\theta }=0\end{array}$$
  $$\hat{\theta }=\frac{z_2+y_4}{z_2+z_1+y_2+y_4}\text{\hspace{1em}(1)}$$
  3）由于$z_1,z_2$未知，但次数服从二项分布，故$z_1\sim B(y_1,\frac{\frac{1}{4}-\frac{\theta }{4}}{\frac{1}{2}-\frac{\theta }{4}}=\frac{1-\theta }{2-\theta })$,同理，$z_2\sim B(y_3,\frac{\theta }{1+\theta })$
  4）EM算法：
  • 第一步（E步）：求期望；目的：消去潜在变量$z_1,z_2$
    $$E(z_1)=y_1\frac{1-\theta }{2-\theta },E(z_2)=y_3\frac{\theta }{1+\theta }$$
    对$(1)$式两边求期望，带入即可.
  • 第二步（M步）：求最大
    带入后得：$$\hat{\theta }=\frac{y_3\frac{\theta }{1+\theta }+y_4}{y_3\frac{\theta }{1+\theta }+y_1\frac{1-\theta }{2-\theta }+y_2+y_4}$$
    迭代法求：
    $${\theta }^{i+1}=\frac{y_3\frac{{\theta }^i}{1+{\theta }^i}+y_4}{y_3\frac{{\theta }^i}{1+{\theta }^i}+y_1\frac{1-{\theta }^i}{2-{\theta }^i}+y_2+y_4}$$

2.引例2²

3枚硬币A，B，C出现正面的概率分别为：$\pi ,p,q$.现做如下实验：先抛硬币A，若出现正面选B，若出现反面选C，再投掷选出的硬币，出现正面记为1，反面记为0；独立重复$n$次实验,假设只能观测掷硬币的结果，不能观测掷硬币的过程，如何求出$\pi ,p,q$？
$$\begin{array}{rl}P(y|\theta )& =\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$
其中$y$为一次观测的结果，$z$为隐变量，即中间A的结果，$\theta =(\pi ,p,q)$.
将观测数据表示为：$Y=(y_1,y_2,\dots ,y_n{)}^{\mathrm{T}}$,隐变量表示为：$Z=(z_1,z_2,\dots ,z_n{)}^{\mathrm{T}}$,则：
$$P(Y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$
根据最大似然估计：
$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )$$
由于这个问题无解析解，故用数值方法迭代法求解，即EM法.
以下介绍EM算法的理论由来部分.

二.证明EM的收敛性³

1.单个高斯

当总体$X\sim N(\mu ,{\sigma }^2),x_i\stackrel{iid}{\sim }X,i=1,2,\dots ,n$,令$\theta =(\mu ,{\sigma }^2)$
$$\hat{\theta }=\arg \underset{\theta }{\max }\sum _{i=1}^n\log N(x_i|\mu ,{\sigma }^2)\text{\hspace{1em}(2)}$$
当总体服从单个高斯分布时，易根据最大似然估计法求得：$\hat{\mu }=\bar{x},\widehat{{\sigma }^2}=S^2$;
其中$L(\theta |x_1,x_2,\dots ,x_n)={\sum }_{i=1}^n\log N(x_i|\mu ,{\sigma }^2)$称为对数似然函数.

2.高斯混合模型

当总体服从混合高斯模型时，假设有$k$个高斯模型，样本$x_i,i=1,2,\dots ,n$ ,$\theta =({\mu }_1,{\mu }_2,\dots ,{\mu }_k,{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_k^2,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{k-1})$,则$x_i$出现的概率为$k$个高斯的叠加，即：
$$\begin{array}{rl}P(x_i|\theta )& =\sum _{j=1}^k{\lambda }_{j}N({\mu }_j,{\sigma }_j^2)\\ \mathrm{s}.\mathrm{t}.\sum _{j=1}^k{\lambda }_j& =1\end{array}$$
若使用最大似然估计，则得(即用$P(x_i|\theta )$代替$(2)$式的$N(x_i|\mu ,{\sigma }^2)$)：
$$\hat{\theta }=\arg \underset{\theta }{\max }\sum _{i=1}^n\log \sum _{j=1}^k{\lambda }_{j}N({\mu }_j,{\sigma }_j^2)$$
由于对每一个参数求导为零是一件很困难的事，所以使用迭代的方法(EM法)求解$\hat{\theta }$,迭代公式为：
$$\begin{array}{rl}{\theta }^{(g+1)}& =\arg \underset{\theta }{\max }{\int }_Z\log P(X,Z|\theta )P(Z|X,{\theta }^g)\mathrm{d}Z\\ \mathrm{s}.\mathrm{t}.\log P(X|{\theta }^{(g+1)})& \ge \log P(X|{\theta }^g)\end{array}\text{\hspace{1em}(4)}$$
其中$Z$为隐变量集合，$X$为数据集合

3.收敛性证明

即证明：$\log P(X|{\theta }^{(g+1)})\ge \log P(X|{\theta }^g)$
证明：
由：$$\log P(X|\theta )=\log P(X,Z|\theta )-\log P(Z|X,\theta )\text{\hspace{1em}(3)}$$

因为P(AB)=P(A)P(B|A),故$\log P(A)=\log P(AB)-\log P(B|A)$,两边同时加上$\theta$即可

对$(3)$式两边对分布$P(Z|X,{\theta }^g)$求期望：
$${\int }_Z\log P(X|\theta )P(Z|X,{\theta }^g)\mathrm{d}Z={\int }_Z\log P(X,Z|\theta )P(Z|X,{\theta }^g)\mathrm{d}Z-{\int }_Z\log P(Z|X,\theta )P(Z|X,{\theta }^g)\mathrm{d}Z$$
$$左端=\log P(X|\theta )$$
令$Q(\theta ,{\theta }^g)={\int }_Z\log P(X,Z|\theta )P(Z|X,{\theta }^g)\mathrm{d}Z$
$H(\theta ,{\theta }^g)={\int }_Z\log P(Z|X,\theta )P(Z|X,{\theta }^g)\mathrm{d}Z$
故：$$右端=Q(\theta ,{\theta }^g)-H(\theta ,{\theta }^g)$$
假设：$\forall \theta ,都有H({\theta }^g,{\theta }^g)\ge H(\theta ,{\theta }^g)$,得：$H({\theta }^g,{\theta }^g)\ge H({\theta }^{(g+1)},{\theta }^g)$;又由$(4)$式，得$Q({\theta }^g,{\theta }^g)\le Q({\theta }^{(g+1)},{\theta }^g)$故：
$$Q({\theta }^g,{\theta }^g)-H({\theta }^g,{\theta }^g)\le Q({\theta }^{(g+1)},{\theta }^g)-H({\theta }^{(g+1)},{\theta }^g)$$
由此可得$\log P(X|{\theta }^{(g+1)})\ge \log P(X|{\theta }^g)$.
现证满足假设：$\forall \theta ,都有H({\theta }^g,{\theta }^g)\ge H(\theta ,{\theta }^g)$
证明：
$$\begin{array}{rl}H({\theta }^g,{\theta }^g)-H(\theta ,{\theta }^g)& ={\int }_Z\log P(Z|X,{\theta }^g)P(Z|X,{\theta }^g)\mathrm{d}Z-{\int }_Z\log P(Z|X,\theta )P(Z|X,{\theta }^g)\mathrm{d}Z\\ & ={\int }_Z-\log \frac{P(Z|X,\theta )}{P(Z|X,{\theta }^g)}P(Z|X,{\theta }^g)\mathrm{d}Z\\ & \stackrel{\ast }{\ge }0\end{array}$$

（ *）步得由来：
$f(x)=-\log x$是一个凸函数，即满足定义域内$\forall x,y，\lambda \in [0,1]$，$$\mathrm{s}.\mathrm{t}.\lambda f(x)+(1-\lambda )f(y)\ge f(\lambda x+(1-\lambda )y)$$
即：两点连线在函数的上方.
还可以将式子两边视作期望：函数的期望大于等于期望的函数；
故函数$-\log \frac{P(Z|X,\theta )}{P(Z|X,{\theta }^g)}$的期望等于它期望的函数：
∫ Z − log ⁡ P ( Z ∣ X , θ ) P ( Z ∣ X , θ g ) P ( Z ∣ X , θ g ) d Z ≥ − log ⁡ { ∫ Z P ( Z ∣ X , θ ) P ( Z ∣ X , θ g ) P ( Z ∣ X , θ g ) d Z } ≥ − log ⁡ 1 ≥ 0 \begin{aligned}\int_Z -\log \frac{P(Z \vert X,\theta)}{P(Z \vert X,\theta^g)} P(Z|X,\theta^g)\mathrm{d}Z &\geq -\log\{\int_Z \frac{P(Z \vert X,\theta)}{P(Z \vert X,\theta^g)}P(Z|X,\theta^g)\mathrm{d}Z\} \\ &\geq -\log1 \\&\geq 0\end{aligned} ∫Z​−logP(Z∣X,θg)P(Z∣X,θ)​P(Z∣X,θg)dZ​≥−log{∫Z​P(Z∣X,θg)P(Z∣X,θ)​P(Z∣X,θg)dZ}≥−log1≥0​

三.EM算法的步骤

由：$$\begin{array}{rl}{\theta }^{(g+1)}& =\arg \underset{\theta }{\max }{\int }_Z\log P(X,Z|\theta )P(Z|X,{\theta }^g)\mathrm{d}Z\end{array}$$
我们只需要得知每个模型的$P(X,Z|\theta )$和$P(Z|X,{\theta }^g)$即可迭代求出$\hat{\theta }$

四.混合高斯举例

• 求$P(X,Z|\theta )$
  $$P(X,Z|\theta )=\prod _{i=1}^{n}P(x_i,z_i|\theta )=\prod _{i=1}^{n}P(x_i|z_i,\theta )P(z_i|\theta )=\prod _{i=1}^n{\lambda }_{z_i}N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)$$
• 求$P(Z|X,{\theta }^g)$
  $$\begin{array}{rl}P(Z|X,{\theta }^g)& =\prod _{i=1}^{n}P(z_i|x_i,{\theta }^g)\\ & \stackrel{\ast \ast }{=}\frac{{\lambda }_{z_i}N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)}{{\sum }_{z_i=1}^k{\lambda }_{z_i}N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)}\end{array}$$

(**)是由全概率公式：$$P(A|B)=\frac{P(B|A)P(A)}{{\sum }_{i=1}^{n}P(B|A)P(A)}$$
推导而来

• 带入$(4)$式：

  • E-step(即求期望步骤)：
    $$\begin{array}{rl}原式& =\sum _{z_1=1}^k\sum _{z_2=1}^k\cdots \sum _{z_n=1}^k[\sum _{i=1}^n(\log {\lambda }_{z_i}+\log N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2))\prod _{i=1}^{n}P(z_i|x_i,{\theta }^g)]\\ & \stackrel{\ast \ast \ast }{=}\sum _{i=1}^n\sum _{z_i=1}^k(\log {\lambda }_{z_i}+\log N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2))P(z_i|x_i,{\theta }^g)\end{array}$$

  (***)的由来：
  令$f_i(z_i)=\log {\lambda }_{z_i}+\log N(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2),P(z_1,z_2,\dots ,z_n)={\prod }_{i=1}^{n}P(z_i|x_i,{\theta }^g)$
  $$\begin{array}{rl}原式& =\sum _{z_1=1}^k\sum _{z_2=1}^k\cdots \sum _{z_n=1}^k(f_1(z_1)+f_2(z_2)+f_n(z_n))P(z_1,z_2,\dots ,z_n)\\ & =\sum _{z_1=1}^k\sum _{z_2=1}^k\cdots \sum _{z_n=1}^k{f}_1(z_1)P(z_1,z_2,\dots ,z_n)+\dots \\ & =\sum _{z_1=1}^k{f}_1(z_1)\sum _{z_2=1}^k\cdots \sum _{z_n=1}^{k}P(z_1,z_2,\dots ,z_n)+\dots \\ & =\sum _{z_1=1}^k{f}_1(z_1)P(z_1)+\dots \end{array}$$

  • M-step(argmax步骤)
    • 求${\lambda }_{z_i}$：
      $$\frac{\mathrm{d}\log {\lambda }_{z_i}P(z_i|x_i,{\theta }^g)}{\mathrm{d}{\lambda }_{z_i}}\stackrel{令}{=}0$$
      $$\mathrm{s}.\mathrm{t}.\sum _{z_i=1}^k=1$$
      用拉格朗日乘数法求解即可.
      解得：${\lambda }_{z_i}=\frac{1}{n}{\sum }_{i=1}^{n}P(z_i|x_i,\theta )$
      含义：所有的高斯的占比的和求平均.
    • 求${\mu }_i,{\sigma }_i^2$:
      用矩阵求导为零计算所得.
      综上所述：
      $$\begin{array}{rl}{\lambda }_l^{(g+1)}& =\frac{1}{n}\sum _{l=1}^{n}P(l|x_l,{\theta }^g)\\ {\mu }_l^{(g+1)}& =\frac{{\sum }_{l=1}^n{x}_{l}P(l|x_l,{\theta }^g)}{{\sum }_{l=1}^{n}P(z_l|x_l,{\theta }^g)}\\ {{\sigma }_l^2}^{(g+1)}& =\frac{{\sum }_{l=1}^n(x_l-{\mu }_l^{l+1})(x_l-{\mu }_l^{l+1}{)}^{\mathrm{T}}P(l|x_l,{\theta }^g)}{{\sum }_{l=1}^{n}P(l|x_l,{\theta }^g)}\end{array}$$

五.EM算法推广——GEM

…

---

1. 此引例来自于mathtsing ↩︎

2. 此引例来自于李航的《统计学习方法》 ↩︎

3. 以下内容来源于徐亦达的机器学习 ↩︎