《统计学习方法》——EM算法及其推广(上)_统计方法em例题

引言

EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计。
理解EM算法需要很多概率论的知识，所以下面先贴出所需要的知识。便于对后文的理解。

补充知识

期望

对于离散型随机变量$X$的概率分布为 p i = p { X = x i } p_i=p\{X=x_i\} pi​=p{X=xi​}，期望$E(X)$为：
$$E(X)=\sum _i{x}_i{p}_i$$
其中$p_i$是概率，满足两个条件$1\le p_i\le 0$且${\sum }_i{p}_i=1$。

若连续型随机变量$X$的概率密度函数为$f(x)$，则期望$E(X)$为：
$$E(X)={\int }_{-\infty }^{+\infty }xf(x)dx$$
设$Y=g(X)$，若$X$是离散型随机变量，则：
$$E(Y)=\sum _{i}g(x_i)p_i$$
若$X$是连续型随机变量，则：
$$E(X)={\int }_{-\infty }^{+\infty }g(x)f(x)dx$$

$C$为常数，$X$和$Y$是两个随机变量，以下是期望的性质：

• $E(C)=C$
• $E(CX)=CE(X)$
• $E(X+Y)=E(X)+E(Y)$

联合概率与边缘概率

假设有离散型随机变量X与Y

• $P(X=a,Y=b)$用于表示$X=a$且$Y=b$的概率，这种包含多个条件且所有条件同时成立的概率称为联合概率。
• 与之对应地，$P(X=a)$或$P(Y=b)$这种只与单个随机变量有关的概率称为边缘概率。

联合概率一览表称为联合分布，边缘概率一览表称为边缘分布。

联合概率与边缘概率的关系如下：
$$P(X=a)=\sum _{b}P(X=a,Y=b)$$
上式相当于是固定$X$，$Y$取所有的取值($Y$跑一圈)。也可以固定$Y$，让$X$取所有的值：

$$P(Y=b)=\sum _{a}P(X=a,Y=b)$$
来看一下简单的题理解一下。

例题1.抛两枚硬币(出现正反面的概率皆为50%)，且随机变量X,Y分别表示正面及反面出现的次数，求X与Y的联合概率及X的边缘概率。

这是离散型的，我们可以写出所有可能：(正，正)、(正，反)、(反，正)、(反，反)

然后画一个表格，写出所有可能的概率。

这个就是X与Y的联合分布。
X的边缘分布则为：

其中$P(X=0)=P(X=0,Y=0)+P(X=0,Y=1)+P(X=0,Y=2)=0+0+\frac{1}{4}=\frac{1}{4}$

我们可以像上面这样通过联合分布计算边缘分布。
联合分布还有一条性质就是：
$$\sum _a\sum _{b}P(X=a,Y=b)=1$$

条件概率

条件概率的定义为：
$$P(Y=b|X=a)=\frac{P(X=a,Y=b)}{P(X=a)}$$

条件概率有如下性质：
在条件$Y=b$下$X$的条件分布也是一种"X的概率分布"，因此穷举$X$可取的值后，所有这些值对应的概率之和为$1$。即：
$$\sum _{a}P(X=a|Y=b)=1$$

要看清楚$X$和$Y$哪个是条件，如果是${\sum }_{a}P(Y=b|X=a)$的值，则不一定是1。

下面简单地证明下上式：
$$\sum _{a}P(X=a|Y=b)=\frac{{\sum }_{a}P(X=a,Y=b)}{P(Y=b)}=\frac{P(Y=b)}{P(Y=b)}=1$$

三个或更多的随机变量

含有三个或更多的随机变量时，条件概率的定义依然与之前相同。

• $Y=b$且$Z=c$时，$X=a$的条件概率

$$P(X=a|Y=b,Z=c)=\frac{P(X=a,Y=b,Z=c)}{P(Y=b,Z=c)}$$

• $Z=c$时，$X=a$且$Y=b$的条件概率

$$P(X=a,Y=b|Z=c)=\frac{P(X=a,Y=b,Z=c)}{P(Z=c)}$$

通常，我们可以像下面这样分解联合概率。
$$P(X=a,Y=b,Z=c)=P(X=a|Y=b,Z=c)P(Y=b|Z=c)P(Z=c)$$

☆条件联合分布的分解☆

标了☆表示这个是重点，上面说了这么多，就是为了引入这一小节。

$$P(X=a,Y=b|Z=c)=P(X=a|Y=b,Z=c)P(Y=b|Z=c)\text{\hspace{1em}(p1)}$$
乍一看好像两边不相等，但是看完这一小节后，保证你乍一看，两边就是相等的。

上式$(p1)$左侧为
$$\frac{P(X=a,Y=b,Z=c)}{P(Z=c)}$$
$(p1)$右侧可以写为(为右侧两项分别应用条件概率的定义)：
$$\frac{P(X=a,Y=b,Z=c)}{P(Y=b,Z=c)}\cdot \frac{P(Y=b,Z=c)}{P(Z=c)}=\frac{P(X=a,Y=b,Z=c)}{P(Z=c)}$$
这样就证明了等式$(p1)$的成立，也可以把等式两边同乘以$P(Z=c)$来证明。

还可以换一种方式来理解，$(p1)$式中每一个$P$都附有$Z=c$的条件。即条件$Z=c$是整个式子的大前提。因此，我们可以先声明这一前提，表示之后将基于这一大前提分解概率，而不必每次强调$Z=c$。

于是，就转化为了以下式子：
$$P(X=a,Y=b)=P(X=a|Y=b)P(Y=b)$$

好了，下面正式进入EM算法的理解。

Jensen不等式

Jensen不等式和凸函数的定义有关。首先看凸函数的定义。

设$C$是非空凸集，$f$是定义在$C$上的函数，如果对任意的$x_1,x_2\in C,\lambda \in (0,1)$，均有：
$$f(\lambda x_1+(1-\lambda )x_2)\le \lambda f(x_1)+(1-\lambda )f(x_2)\text{\hspace{1em}(p2)}$$
则称$f$为$C$上的凸函数；

凸函数示意图，来自convex.indigits.com

如上图所示($\theta$即$\lambda$)，即凸函数任意两点之间直线位于两点之间曲线的上方。

这也是Jensen不等式的两点形式，我们可以将其推广到一般的情况，即：
$$f(\sum _{i=1}^n{\lambda }_i{x}_i)\le \sum _{i=1}^n{\lambda }_{i}f(x_i))$$
其中${\sum }_{i=1}^n{\lambda }_i=1$。

如果为凹函数，则不等式符号取反，有：
$$f(\sum _{i=1}^n{\lambda }_i{x}_i)\ge \sum _{i=1}^n{\lambda }_{i}f(x_i))$$

下面我们来证明一下。对于$n=2$的情况，由凸函数的定义，已经成立了。

下面，我们通过数学归纳法的思想来证明。

假设在$n=k$的情况下，不等式成立，即
$$f(\sum _{i=1}^k{\lambda }_i{x}_i)\le \sum _{i=1}^k{\lambda }_{i}f(x_i)\text{\hspace{1em}(p3)}$$
成立。

下面我们来证明当$n=k+1$时，公式仍然成立，即
$$f(\sum _{i=1}^{k+1}{\lambda }_i{x}_i)\le \sum _{i=1}^{k+1}{\lambda }_{i}f(x_i))$$
首先把函数里面$k+1$个数分成两部分，
$$\begin{array}{rl}\sum _{i=1}^{k+1}{\lambda }_i{x}_i& =\sum _{i=1}^k{\lambda }_i{x}_i+{\lambda }_{k+1}{x}_{k+1}\\ & =(1-{\lambda }_{k+1})\sum _{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}{x}_i+{\lambda }_{k+1}{x}_{k+1}\end{array}$$
令$x^{\ast }={\sum }_{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}{x}_i$，则：
$$f\left(\sum _{i=1}^{k+1}{\lambda }_i{x}_i\right)=f((1-{\lambda }_{k+1})x^{\ast }+{\lambda }_{k+1}{x}_{k+1})$$
$x^{\ast }$相当于一个新的点，而$1-{\lambda }_{k+1}+{\lambda }_{k+1}=1$，变成了满足公式$(p2)$的形式，因此代入得：
$$\begin{array}{rl}f\left(\sum _{i=1}^{k+1}{\lambda }_i{x}_i\right)& =f((1-{\lambda }_{k+1})x^{\ast }+{\lambda }_{k+1}{x}_{k+1})\\ & \le (1-{\lambda }_{k+1})f(x^{\ast })+{\lambda }_{k+1}f(x_{k+1})\\ & =(1-{\lambda }_{k+1})f\left(\sum _{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}{x}_i\right)+{\lambda }_{k+1}f(x_{k+1})\end{array}\text{\hspace{1em}(p4)}$$
下面我们针对函数中求和这一项，我们假设$n=k+1$，因此有${\sum }_{i=1}^k{\lambda }_i+{\lambda }_{k+1}=1$，即有${\sum }_{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}=1$成立，根据公式$(p3)$，有：
$$f\left(\sum _{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}{x}_i\right)\le \sum _{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}f(x_i)$$
把该结果代入$(p4)$，则有：
$$\begin{array}{rl}f\left(\sum _{i=1}^{k+1}{\lambda }_i{x}_i\right)& \le (1-{\lambda }_{k+1})\left(\sum _{i=1}^k\frac{{\lambda }_i}{1-{\lambda }_{k+1}}f(x_i)\right)+{\lambda }_{k+1}f(x_{k+1})\\ & =\sum _{i=1}^k{\lambda }_{i}f(x_i)+{\lambda }_{k+1}f(x_{k+1})\\ & =\sum _{i=1}^{k+1}{\lambda }_{i}f(x_i)\end{array}$$
上面$\frac{1}{1-{\lambda }_{k+1}}$与求和中的$i$无关，可以提到求和符号外面，刚好与外面的$1-{\lambda }_{k+1}$相乘得1。证毕。

对于对数函数，即$f(\cdot )$为$\log (\cdot )$，注意$\log$函数为凹函数，因此我们可以写成对应的Jensen不等式：
$$\log (\sum _{i=1}^n{\lambda }_i{x}_i)\ge \sum _{i=1}^n{\lambda }_i\log (x_i)$$

EM算法的引入

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

EM算法

首先介绍一个使用EM算法的例子。
例9.1
假设有3 枚硬币，分别记作A,B,C。这些硬币正面出现的概率分别是$\pi$， $p$ 和$q$。进行如下掷硬币试验: 先掷硬币A ，根据其结果选出硬币B或硬币C ，正面选硬币B ，反面边硬币C; 然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0; 独立地重复n 次试验(这里， n= 10) ，观测结果如下:
$$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数（即$\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}\text{\hspace{1em}(9.1)}$$
这里，随机变量$y$是观测变量，表示一次试验观测的结果是1 或0; 随机变量$z$是隐变量，表示未观测到的掷硬币A 的结果;
$\theta =(\pi ,p,q)$是模型参数。这一模型是以上数据的生成模型。注意，随机变量$y$的数据可以观测，随机变量$z$ 的数据不可观测。

第一次看到这个式子，半天没明白，然后只能跑回去补概率论的知识。下面写出本人的理解，如有误，请斧正。
首先$P(y|\theta )={\sum }_{z}P(y,z|\theta )$。这个式子左右两边条件中都有一个$\theta$，我们可以去掉它，就是联合概率和边缘概率关系的式子。也可以进行简单的证明。
$$\begin{array}{rl}& P(y|\theta )=\sum _{z}P(y,z|\theta )\\ & \frac{P(y,\theta )}{P(\theta )}=\frac{{\sum }_{z}P(y,z,\theta )}{P(\theta )}\\ & \frac{P(y,\theta )}{P(\theta )}=\frac{P(y,\theta )}{P(\theta )}\end{array}$$
而对于${\sum }_{z}P(y,z|\theta )={\sum }_{z}P(z|\theta )P(y|z,\theta )$则是补充知识中条件联合分布的分解这一小节介绍的内容。条件中都有$\theta$，我们直接去掉它，就好理解了，这里就不证明了。

下面来看式$(9.1)$最后一个等式是怎么来的。

(上图参考了程序员的数学2 概率论)
我们可以画一个类似上面的正方形来理解，把概率转换为面积的思想来理解，整个正方形的面积为1。首先在该正方形中间画一条红线，分成了两个矩阵。
A硬币出现正面的概率为$\pi$(这只是一个符号，所有的概率都是$[0,1]$之间，不是圆周率那个$\pi$)，它就占了红线上方的那个矩形。A硬币出现反面的概率就是$1-\pi$，即线色下面那个矩形。

在A正的情况下，硬币B正面概率为$p$，就是上图深灰色区域。在A反的情况下，C正的概率为$q$，为上图浅灰色区域。
式$(9.1)$是一次试验的结果。$y$为观测结果，$y$要么为0要么为1。

• 当y=1时，这次实验的概率就是上图灰色面积之和：

$$\pi p+(1-\pi )q$$

• 当y=0时，这次实验的概率就是上图白色面积之和：

$$\pi (1-p)+(1-\pi )(1-q)$$

写到一起就是最终的等式。

也可以直接根据${\sum }_{z}P(z|\theta )P(y|z,\theta )$写出最终结果，这里$\theta$是固定的，可以先不看。
$P(z=1)$就是$\pi$，$P(z=0)=1-\pi$。
乘以$P(y|z)$，若$z=1$，则考虑硬币B。根据$y$的观测结果，则乘以硬币B相应的概率$p$或$1-p$。写到一起也可以得到$\pi p^y(1-p{)}^{1-y}$。

当y=1时，$\pi p^y(1-p{)}^{1-y}$就是$\pi p$；当y=0时，就是$\pi (1-p)$。虽然写到了一起，但是根据实际结果只会出现一个。

所以将观测数据表示为$Y=(Y_1,Y_2,\cdots ,Y_n{)}^T$ ，未观测数据表示为$Z=(Z_1,Z_2,\cdots ,Z_n{)}^T$则观测数据的似然函数为

$$P(Y|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\text{\hspace{1em}(9.2)}$$
因为每次试验结果都是独立的，出现$Y$的概率就可以写成每次试验概率连乘：
$$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}]\text{\hspace{1em}(9.3)}$$
考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即
$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )\text{\hspace{1em}(9.4)}$$
就是求得使得式$(9.3)$概率最大的参数$\theta$。加了对数使连乘变成连加。
EM算法就是可以用于求解这个问题的一种迭代算法。下面给出针对以上问题的EM算法。

这里先给出EM算法的一个计算过程，其推导可见下文。
EM算法首先选取参数的初始值，记作${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，然后通过下面的步骤迭代计算参数的估计值，直到收敛为止。

所谓迭代，就是第$i+1$次的参数值可以通过第$i$次来写出。第$i+1$次迭代过程如下：

• E步：计算在模型参数${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据$y_j$来自掷硬币B的概率

$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_j}}\text{\hspace{1em}(9.5)}$$

不要被这个式子吓到了，其实就是第$j$份观测数据$y_j$下，分母是公式$(9.1)$的式子。

再啰嗦一下，上标$(i)$表示第$i$次迭代得到的参数值。而分子是来自硬币B的概率，根据$y_j$的结果把正面和反面的式子写到了一起。再来看分母，分母就是来自硬币B的概率加上来自硬币C的概率。

• M步：计算模型参数的新估计量

$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}\text{\hspace{1em}(9.6)}$$

重点来了，我们说$\pi$是硬币A正面朝上的概率，即得到的结果$y$来自硬币B的概率，所以我们对E步计算出来的来自硬币B的概率求一个均值，就得到了更新后的$\pi$。
$$p^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}\text{\hspace{1em}(9.7)}$$
结合面积图的话，就是计算深灰色区域。这里更新的是B正面朝上的概率。因为 y ∈ { 0 , 1 } y \in \{0,1\} y∈{0,1}，所以可以这么写。B正面向上的概率除以B出现的总概率(加上B反面朝上的概率)。
$$q^{(i+1)}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})}\text{\hspace{1em}(9.8)}$$
这个就很容易理解了，$\mu$是来自硬币B的概率，那么$(1-\mu )$就是来自硬币C的概率。剩下的和公式$(9.7)$类似。

下面通过程序来实现上面的公式，并赋予不同的初值，看得到的结果会怎样。

import numpy as np
Y = np.array([1,1,0,1,0,0,1,0,1,1])
pi,p,q = 0.5,0.5,0.5 #初始值，可变
count = 2 #简单实现，先固定迭代次数
pis,ps,qs = [pi],[p],[q] #保存历史计算结果
u = np.zeros(len(Y),dtype=np.float)

for c in range(count):
    for i,y in enumerate(Y):
        B = pis[c] * ps[c]**y * (1-ps[c])**(1-y)
        C = (1-pis[c]) * qs[c]**y * (1-qs[c])**(1-y)
        u[i] = B / (B +C)
    pi = np.mean(u) #π
    p = np.sum(u * Y) / np.sum(u)
    q = np.sum((1-u) * Y) / np.sum(1-u)
    pis.append(pi)
    ps.append(p)
    qs.append(q)
    print('after iteration %d,pi=%2.4f,p=%2.4f,q=%2.4f' %(c+1,pi,p,q))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

可以看到得到的结果和书上的是一样的。在给定不同的初始值，看得到的结果。

可以看到，EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起成为完全数据，观测数据Y又称为不完全数据。假设给定观测数据Y，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数，那么不完全数据$Y$的似然函数是$P(Y|\theta )$，对数似然函数$L(\theta )=\log P(Y|\theta )$；假设Y和Z的联合概率分布是$P(Y,Z|\theta )$，那么完全数据的对数似然函数是$\log P(Y,Z|\theta )$。

EM算法通过迭代求$L(\theta )=\log P(Y|\theta )$的极大似然顾及。每次迭代包含两步：E步，求期望；M步，求极大化。下面来介绍EM算法。

算法9.1 (EM算法)

输入：观测变量数据Y，隐变量数据Z，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$；
输出：模型参数$\theta$。

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

(2) E步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ & =\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}\text{\hspace{1em}(9.9)}$$

这里，$P(Z|Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布；

$Q(\theta ,{\theta }^{(i)})$是条件分布$P(Z|Y,{\theta }^{(i)})$下的期望。

(3) M步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9.10)}$$
(4) 重复第(2)步和第(3)步，直到收敛。
式$(9.9)$的函数$Q(\theta ,{\theta }^{(i)})$是EM算法的核心，称为$Q$函数(Q function)。

定义 9.1(Q函数) 完全数据的对数似然函数$\log P(Y,Z|\theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为$Q$函数，即
$$Q(\theta ,{\theta }^{(i)})=E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\text{\hspace{1em}(9.11)}$$
对EM算法的几个步骤进行说明：

步骤(1) 参数的初值可以任意选择，但EM算法对初值是敏感的。

步骤(2) E步求$Q(\theta ,{\theta }^{(i)})$。$Q$函数式中$Z$是未观测数据，$Y$是观测数据。$Q(\theta ,{\theta }^{(i)})$中的第一个参数表示要极大化的参数，第二个参数是当前的估计值，即是一个确定的值。每次迭代实际在求$Q$函数及其极大。

步骤(3) M步求$Q(\theta ,{\theta }^{(i)})$的极大化，得到${\theta }^{(i+1)}$，完成一次${\theta }^{(i)}\to {\theta }^{(i+1)}$的迭代。

步骤(4) 给出停止迭代的条件，一般是对较小的正数${ϵ}_1,{ϵ}_2$，若满足
$$||{\theta }^{(i+1)}-{\theta }^{(i)}||\le {ϵ}_1\text{或}||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||\le {ϵ}_2$$
则停止迭代。

EM算法的导出

本节推导出式$(9.9)$的$Q$函数。

我们面对一个含有隐变量的概率模型，目标是极大化观测数据(不完全数据)$Y$关于参数$\theta$的对数似然函数，即极大化式$(9.2)$的似然函数，取对数：
$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Z|\theta )P(Y|Z,\theta )\text{\hspace{1em}(9.12)}$$
这一极大化的主要难点在于其中包含未观测数据($Z$)，并且包含和的对数$\log \sum$形式。

EM算法是通过迭代逐步近似极大化$L(\theta )$的。假设在第$i$次迭代后$\theta$的估计值是${\theta }^{(i)}$。我们希望新估计值$\theta$能使$L(\theta )>L({\theta }^{(i)})$，并逐步达到极大值。因此，可以考虑两者的差：
$$L(\theta )-L({\theta }^{(i)})=\log \left(\sum _{Z}P(Y|Z,\theta )P(Z|\theta )\right)-\log P(Y|{\theta }^{(i)})$$

上式中的第二项$\log P(Y|{\theta }^{(i)})$中的$Y$和${\theta }^{(i)}$都是确定的值，因此不展开。根据我们上面补充知识中介绍的Jessen不等式，对上式中的$\log \sum$转换为$\sum \log$的形式：

$$\begin{array}{rlr}L(\theta )-L({\theta }^{(i)})& =\log \left(\sum _{Z}P(Y|Z,\theta )P(Z|\theta )\right)-\log P(Y|{\theta }^{(i)})& \\ & =\log \left(\sum _Z\frac{P(Z|Y,{\theta }^{(i)})}{P(Z|Y,{\theta }^{(i)})}\cdot P(Y|Z,\theta )P(Z|\theta )\right)-\log P(Y|{\theta }^{(i)})& \\ & =\log \left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\cdot \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\right)-\log P(Y|{\theta }^{(i)})& \\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)})& \\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y|{\theta }^{(i)})& \text{... 因为}\sum _{Z}P(Z|Y,{\theta }^{(i)})=1\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}& \end{array}$$
这里引入$p(Z|Y,{\theta }^{(i)})$，给定$Y,{\theta }^{(i)}$下$Z$的条件概率分布，也是一个概率分布，所以${\sum }_{Z}p(Z|Y,{\theta }^{(i)})=1$。符合Jesson不等式的要求。

令
$$B(\theta ,{\theta }^{(i)})\widehat{=}L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\text{\hspace{1em}(9.13)}$$
相当于把前面的$L({\theta }^{(i)})$移到等式右边，则
$$L(\theta )\ge B(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9.14)}$$
这样就求得$L(\theta )$的一个下界，即$L(\theta )$最小值为$B(\theta ,{\theta }^{(i)})$。将${\theta }^{(i)}$代入式$(9.13)$，有
$$B({\theta }^{(i)},{\theta }^{(i)})\widehat{=}L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y,Z|{\theta }^{(i)})}{P(Y,Z|{\theta }^{(i)})}=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log 1=L({\theta }^{(i)})$$
即
$$L({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})\text{\hspace{1em}(9.15)}$$
说明$B(\theta ,{\theta }^{(i)})$与$L(\theta )$在$\theta ={\theta }^{(i)}$处相交，${\theta }^{(i)}$作为一个新的起点，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可以使$L(\theta )$增大。即为了$L(\theta )$尽可能大地增长，选择使$B(\theta ,{\theta }^{(i)})$达到极大${\theta }^{(i+1)}$，即
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9.16)}$$
现在求${\theta }^{(i+1)}$的表达式，省去对$\theta$极大化来说是常数的项，由式$(9.16)$、$(9.13)$及$(9.10)$有
$$\begin{array}{rl}{\theta }^{(i+1)}& =\arg \underset{\theta }{\max }\left(L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\right)\\ & =\underset{L({\theta }^{(i)})是常数}{\underbrace{\cancel{\arg \underset{\theta }{\max }L({\theta }^{(i)})}}}+\arg \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log (P(Y|Z,\theta )P(Z|\theta ))\right)-\underset{Y,{\theta }^{(i)}已知，且对所有Z进行求和后也是常数}{\underbrace{\cancel{\arg \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log (P(Z,Y|{\theta }^{(i)}))\right)}}}\\ & =\arg \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log (P(Y|Z,\theta )P(Z|\theta ))\right)\\ & =\arg \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )\right)\\ & =\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}\text{\hspace{1em}(9.17)}$$
这就得到了$Q$函数的表达式。上式等价于EM算法包含E步和M步的一次迭代，即求$Q$函数及其极大化。

可以看到，EM算法是通过不断求解下界的极大化来逼近求解对数似然函数极大化的算法。

图9.1给出了EM算法的直观解释，来自原书。图中加粗上方曲线为$L(\theta )$，下方曲线为$B(\theta ,{\theta }^{(i)})$。$B(\theta ,{\theta }^{(i)})$为对数似然函数$L(\theta )$的下界。由式$(9.15)$，两个函数在点$\theta ={\theta }^{(i)}$处相等。EM算法找到下一个点${\theta }^{(i+1)}$使函数$B(\theta ,{\theta }^{(i)})$极大化，也是函数$Q(\theta ,{\theta }^{(i)})$极大化。由于$L(\theta )\ge B(\theta ,{\theta }^{(i)})$，函数$B(\theta ,{\theta }^{(i)})$的增加，保证对数似然函数$L(\theta )$在每次迭代中也是增加的。

EM算法在点${\theta }^{(i+1)}$处重新计算$Q$函数值，进行下一次迭代。在此过程中，对数似然函数$L(\theta )$不断增大，从图可以看出EM算法不能保证找到全局最优值。

EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。那么EM算法得到的估计序列是否收敛？下面给出EM算法收敛性的两个定理。

定理 9.1 设$P(Y|\theta )$为观测数据的似然函数，${\theta }^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$P(Y|{\theta }^{(i)})(i=1,2,\cdots )$为对应的似然函数序列，则$P(Y|{\theta }^{(i)})$是单调递增的，即
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})\text{\hspace{1em}(9.18)}$$
证明 由于
$$P(Y|\theta )=\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$
取对数有
$$\log P(Y|\theta )=\log P(Y,Z|\theta )-\log P(Z|Y,\theta )$$
上述等式两边在条件分布$P(Z|Y,{\theta }^{(i)})$下求期望，得：
$$E_{Z|Y,{\theta }^{(i)}}[\log P(Y|\theta )]=E_{Z|Y,{\theta }^{(i)}}[\log P(Y,Z|\theta )-\log P(Z|Y,\theta )]$$
等式左边有：
$$\begin{array}{rl}{E}_{Z|Y,{\theta }^{(i)}}[\log P(Y|\theta )]& =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y|\theta )\\ & =\log P(Y|\theta )\sum _{Z}P(Z|Y,{\theta }^{(i)})\\ & =\log P(Y|\theta )\end{array}$$
等式右边有：
$$E_{Z|Y,{\theta }^{(i)}}[\log P(Y,Z|\theta )-\log P(Z|Y,\theta )]=\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Z|Y,\theta )$$

由式$(9.11)$
$$Q(\theta ,{\theta }^{(i)})=\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
可知上述等式右边第一项就是$Q(\theta ,{\theta }^{(i)})$。

令
$$H(\theta ,{\theta }^{(i)})=\sum _Z\log P(Z|Y,\theta )P(Z|Y,{\theta }^{(i)})\text{\hspace{1em}(9.19)}$$
上述等式左边等于右边，于是对数似然函数可以写成
$$\log P(Y|\theta )=Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9.20)}$$

为了证明式$(9.18)$，在式$(9.20)$中分别取$\theta$为${\theta }^{(i+1)}$和${\theta }^{(i)}$并相减，有
$$\begin{array}{rl}\log P& (Y|{\theta }^{(i+1)})-\log P(Y|{\theta }^{(i)})\\ & =Q({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i+1)},{\theta }^{(i)})-[Q({\theta }^{(i)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})]\\ & =[Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})]-[H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})]\end{array}\text{\hspace{1em}(9.21)}$$
我们只要证明式$(9.21)$两端是非负的，就能。式$(9.21)$右端的第1项，由于${\theta }^{(i+1)}$使$Q(\theta ,{\theta }^{(i)})$达到极大，所以有
$$Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})\ge 0\text{\hspace{1em}(9.22)}$$
其第2项，由式$(9.19)$可得：
$$\begin{array}{rl}H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})& =\sum _Z\log P(Z|Y,{\theta }^{(i+1)})P(Z|Y,{\theta }^{(i)})-\sum _Z\log P(Z|Y,{\theta }^{(i)})P(Z|Y,{\theta }^{(i)})\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\left(\log \frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}\right)\\ & \le \log \left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\cdot \frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}\right)\\ & =\log \left(\sum _{Z}P(Z|Y,{\theta }^{(i+1)}\right)\\ & =\log 1\\ & =0\end{array}\text{\hspace{1em}(9.23)}$$
这里再次利用了Jessen不等式，由式$(9.22)$和式$(9.23)$可知式$(9.21)$右端是非负的。

定理 9.2 设$L(\theta )=\log P(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$L({\theta }^{(i)})(i=1,2\cdots )$为对应的对数似然函数序列。

(1) 如果$P(Y|\theta )$有上界，则$L({\theta }^{(i)})=\log P(Y|{\theta }^{(i)})$收敛到某一值$L^{\ast }$；

(2) 在函数$Q(\theta ,{\theta }^{\prime })$与$L(\theta )$满足一定条件下，由EM算法得到的参数估计序列${\theta }^{(i)}$的收敛值${\theta }^{\ast }$是$L(\theta )$的稳定点。