统计计算二|EM算法（Expectation-Maximization Algorithm，期望最大化算法）

本文主要是介绍统计计算二|EM算法（Expectation-Maximization Algorithm，期望最大化算法），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

系列文章目录

统计计算一|非线性方程的求解

---

讲解EM方法的定义，使用场景，Q函数上升证明，实例运用，以及特殊情况下EM的一些变种。

一、基本概念

（一）极大似然估计和EM算法

• 频率学派通过极大化观测数据的对数似然函数，得到参数的估计。但是在有些情况下，得到的数据不是完整数据，还有一部分数据缺失。这导致无法直接极大化对数似然函数得到感兴趣的参数。
• EM 算法就是用来解决隐含数据存在情况下的参数估计问题。它是一种迭代优化策略, 它是受缺失思想，以及考虑给定已知项下缺失项的条件分布而激发产生的。

（二）EM算法的基本思想

EM算法是一种用于含有隐变量（或潜在变量）的概率模型参数估计的迭代方法。通过反复执行期望（E步）和最大化（M步）两个步骤，使模型参数的对数似然函数逐步逼近最大值。EM算法广泛应用于数据缺失情况下的最大似然估计和贝叶斯推断。

• E步：基于观测数据和上一步估计出的参数值，根据设定的给定已知项下缺失项的条件分布，预测完整数据的对数似然函数
• M步：极大化预测的完整数据的对数似然函数，得到感兴趣参数的估计
• E 步和 M 步反复迭代，直到估计参数基本无变化，算法收敛，得到最终的参数估计

（三）定义

对于给定的观测数据$x$，我们的目标是最大化观测数据的似然函数$L(\theta ,x)$，这个似然函数通常会难以处理，而采用$Y|\theta$和$Z|(x,\theta )$的密度则较容易处理。

1、缺失数据, 边际化和符号

定义$Y=(X,Z)$为完全数据，其中$X$为观测数据，$Z$为未观测到的缺失数据。设$f_X(x|\theta )$和$f_Y(y|\theta )$分别为观测数据和完全数据的密度。在给定观测数据下，缺失数据的条件密度为$f_{Z|X}(z|x;\theta )=f_Y(y|\theta )/f_X(x|\theta )$.

2、Q函数

设${\theta }^{(t)}$表示在迭代t时估计的最大值点（$t=0,1,...$），定义$Q(\theta |{\theta }^{(t)})$为观测数据$X=x$和当前迭代估计出的参数值${\theta }^{(t)}$条件下，完全数据的对数似然的对隐变量的条件期望，即：
Q ( θ ∣ θ ( t ) ) = E { l o g L ( θ ∣ Y ) ∣ x , θ ( t ) } = E { l o g f Y ( y ∣ θ ) ∣ x , θ ( t ) } = ∫ [ l o g f Y ( y ∣ θ ) ] f Z ∣ X ( z ∣ x , θ ( t ) ) d z \begin{equation*} \begin{aligned} Q(\theta|\theta^{(t)}) &= E\{log L(\theta|Y)|x,\theta^{(t)}\} \\ &= E\{log f_Y(y|\theta)|x,\theta^{(t)}\} \\ &= \int [log f_Y(y|\theta)]f_{Z|X}(z|x,\theta^{(t)})dz \\ \end{aligned} \end{equation*} Q(θ∣θ(t))​=E{logL(θ∣Y)∣x,θ(t)}=E{logfY​(y∣θ)∣x,θ(t)}=∫[logfY​(y∣θ)]fZ∣X​(z∣x,θ(t))dz​​

EM从${\theta }^0$开始，在E步和M步之间交替：

• 定义完全数据的似然函数
• E步：计算条件期望 Q ( θ ∣ θ ( t ) ) = E z { l o g f Y ( y ∣ θ ) ∣ x , θ ( t ) } Q(\theta|\theta^{(t)})=E_z\{log f_Y(y|\theta)|x,\theta^{(t)}\} Q(θ∣θ(t))=Ez​{logfY​(y∣θ)∣x,θ(t)}
• M步：寻求$\theta$来最大化$Q(\theta |{\theta }^{(t)})$，设${\theta }^{(t+1)}$为找到的点
• 返回E步，直到满足某停止规则为止

EM算法的停止规则通常依赖于$({\theta }^{(t+1)}-{\theta }^{(t)}{)}^T({\theta }^{(t+1)}-{\theta }^{(t)})$或$|Q({\theta }^{(t+1)}|{\theta }^{(t)})-Q({\theta }^{(t)}|{\theta }^{(t)})|$

3、混合高斯模型（Gaussian Mixture Model，简称GMM）

混合高斯模型是一种常用的概率模型，用于表示由多个高斯分布（正态分布）混合而成的数据分布。GMM广泛应用于聚类分析、密度估计和模式识别等领域。GMM的核心思想是用多个高斯分布来描述数据的整体分布，每个高斯分布代表数据的一个子群。
$$p(x)=\sum _{k=1}^K{\pi }_{k}N(x|{\mu }_k,{\Sigma }_k)$$

要训练混合高斯模型，即估计参数${\pi }_k$、${\mu }_k$和${\Sigma }_k$，通常使用期望最大化（EM）算法。

问题背景：（定义出完全数据的似然函数）

E步：计算条件期望 Q ( θ ∣ θ ( t ) ) = E z { l o g f Y ( y ∣ θ ) ∣ x , θ ( t ) } Q(\theta|\theta^{(t)})=E_z\{log f_Y(y|\theta)|x,\theta^{(t)}\} Q(θ∣θ(t))=Ez​{logfY​(y∣θ)∣x,θ(t)}

M步：寻求$\theta$来最大化$Q(\theta |{\theta }^{(t)})$，设${\theta }^{(t+1)}$为找到的点

4、一般混合模型

如果一个数据集是由$K$个不同总体组成，密度函数具有如下混合密度形式：（其中：${\sum }_{k=1}^K{\pi }_j=1$，模型参数包括$(\Pi ,\Theta )$）
$$p(x;\theta )=\sum _{k=1}^K{\pi }_k{p}_k(x;{\theta }_j)$$

假定观测到的数据为 X = { x 1 , x 2 , . . . , x n } X=\{x_1,x_2,...,x_n\} X={x1​,x2​,...,xn​}，未观测到的数据为 z = { z 1 , z 2 , . . . , z n } z=\{z_1,z_2,...,z_n\} z={z1​,z2​,...,zn​}，其中 z i ∈ { 1 , . . . , K } z_i\in\{1,...,K\} zi​∈{1,...,K}代表数据的来源，也就是：
$$p(x_i|z_i=k;\Theta )=p_k(x;{\theta }_k),p(z_i=k)={\pi }_j$$
则完全数据的似然函数可写成：
$$L(\Pi ,\Theta |X,Z)=\prod _{i=1}^n\prod _{k=1}^K({\pi }_k{p}_k(x_i;{\theta }_k){)}^{I(z_i=k)}$$

对数似然函数为：

$$l(\Pi ,\Theta |X,Z)=\sum _{i=1}^n\sum _{k=1}^{K}I(z_i=k)[(log{\pi }_k+log{p}_k(x_i;{\theta }_k))]$$

然后求Q函数：
Q ( θ ∣ θ ( t ) ) = E z { l ( Π , Θ ∣ X , Z ) ∣ X , Π ( t ) , Θ ( t ) } = ∑ i = 1 n ∑ k = 1 K E ( I ( z i = k ) ∣ X , Π ( t ) , Θ ( t ) ) [ ( l o g π k + l o g p k ( x i ; θ k ) ) ] \begin{equation*} \begin{aligned} Q(\theta|\theta^{(t)}) &= E_z\{l(\Pi,\Theta|X,Z)|X,\Pi^{(t)},\Theta^{(t)}\} \\ &= \sum_{i=1}^n\sum_{k=1}^KE(I(z_i=k)|X,\Pi^{(t)},\Theta^{(t)})[(log\pi_k+logp_k(x_i;\theta_k))] \\ \end{aligned} \end{equation*} Q(θ∣θ(t))​=Ez​{l(Π,Θ∣X,Z)∣X,Π(t),Θ(t)}=i=1∑n​k=1∑K​E(I(zi​=k)∣X,Π(t),Θ(t))[(logπk​+logpk​(xi​;θk​))]​​

对于其中的条件期望，有：
$$\begin{array}{c}\begin{array}{rl}E(I(z_i=k)|X,{\Pi }^{(t)},{\Theta }^{(t)})& =p(z_i=k|X,{\Pi }^{(t)},{\Theta }^{(t)})\\ & =\frac{p(x_i,z_i=k|{\theta }_k^{(t)})}{{\sum }_{l=1}^{K}p(x_i,z_i=l|{\theta }_l^{(t)})}\\ & =\frac{{\pi }_k^{(t)}p(x_i|{\theta }_k^{(t)})}{{\sum }_{l=1}^K{\pi }_l^{(t)}{p}_l(x_i|{\theta }_l^{(t)})}\end{array}\end{array}$$
代入到Q函数后，求导即可获得更新公式。

三、收敛性

（一）琴生不等式（Jensen’s inequality）

琴生不等式（Jensen’s inequality）是凸函数和期望值理论中的一个基本结果。它给出积分的凸函数值和凸函数的积分值间的关系，在此不等式最简单形式中，阐明了对一平均做凸函数变换，会小于等于先做凸函数变换再平均。

• 若将琴生不等式应用在二点上，就回到了凸函数的基本性质：过一个凸函数上任意两点所作割线一定在这两点间的函数图象的上方
  $$tf(x_1)+(1-t)f(x_2)\ge f(t{x}_1+(1-t)x_2),0\le t\le 1$$
• 如果$X$是随机变量且$\varphi$是凸函数，则：
  $$E[\varphi (X)]\ge \varphi (E[X])$$

（二） EM 算法的收敛性质

①EM算法保证了每次迭代后对数似然函数的值不会减少
观测数据密度的对数可重新表达为：
$$log{f}_X(x|\theta )=log{f}_Y(y|\theta )-log{f}_{Z|X}(z|x;\theta )$$
所以有：
E { l o g f X ( x ∣ θ ) ∣ x , θ ( t ) } = E { l o g f Y ( y ∣ θ ) ∣ x , θ ( t ) } − E { l o g f Z ∣ X ( z ∣ x , θ ) ∣ x , θ ( t ) } E\{log f_X(x|\theta)|x,\theta^{(t)}\}=E\{log f_Y(y|\theta)|x,\theta^{(t)}\}-E\{log f_{Z|X}(z|x,\theta)|x,\theta^{(t)}\} E{logfX​(x∣θ)∣x,θ(t)}=E{logfY​(y∣θ)∣x,θ(t)}−E{logfZ∣X​(z∣x,θ)∣x,θ(t)}
其中期望是关于$Z|(x,{\theta }^{(t)})$求取的，于是：
$$log{f}_X(x|\theta )=Q(\theta |{\theta }^{(t)})-H(\theta |{\theta }^{(t)})$$
其中 H ( θ ∣ θ ( t ) ) = E { l o g f Z ∣ X ( Z ∣ x , θ ) ∣ x , θ ( t ) } H(\theta|\theta^{(t)}) =E\{log f_{Z|X}(Z|x,\theta)|x,\theta^{(t)}\} H(θ∣θ(t))=E{logfZ∣X​(Z∣x,θ)∣x,θ(t)}

通过琴生不等式可以证明$H(\theta |{\theta }^{(t)})$是不断减小的：

由于每一步$Q$增大而$H$减小，所以观测数据的对数似然不断上升。

• 标准的 EM 算法：每次迭代中选择 ${\theta }^{(t+1)}$ 来最大化 $Q(\theta |{\theta }^{(t)})$
• 广义 EM（GEM）：简单选取任一个使得$Q({\theta }^{(t+1)}|{\theta }^{(t)})>Q({\theta }^{(t)}|{\theta }^{(t)})$的 ${\theta }^{(t+1)}$

若 $-\ell (\hat{\theta }|x)$ 是正定的，则 EM 算法有线性收敛，且收敛速度跟缺失信息比例有关。当缺失信息比例比较大，收敛会非常慢。

②EM 策略将优化问题由观测数据的对数似然函数 $l$ 转换到替代函数 $G$

EM 策略将优化问题由观测数据的对数似然函数 $l$ 转换到替代函数 $G$(有效地到 Q，因为后两项与 θ 无关，相当于常数)，这更便于最大化。$G$的最大值点保证了在$l$ 值上的增加.

$$\begin{gathered} H({\theta }^{(t)}|{\theta }^{(t)})-H(\theta |{\theta }^{(t)})\ge 0 \\ l(\theta |x)\ge Q(\theta |{\theta }^{(t)})+l({\theta }^{(t)}|x)-Q({\theta }^{(t)}|{\theta }^{(t)})=G(\theta |{\theta }^{(t)}) \end{gathered}$$
由于$Q(\theta |{\theta }^{(t)})$和$l({\theta }^{(t)}|x)$独立于$\theta$，函数$Q$和$G$在相同的$\theta$达到最大。此外，$G$在${\theta }^{(t)}$与$l$相切，并且在任一处低于$l$，称$G$是$l$的一个劣化函数。

（三）MM 算法 (Minorize-Maximization)

MM算法在每次迭代时找到一个目标函数的下界函数（替代函数），求下界函数的最大值。

• 每个 E 步等同于构造劣化函数 G
• 每个 M 步等同于最大化该函数以给出一个上升的路径.
  

四、方差估计

EM 的本质是极大化观测似然函数。因此, EM 参数估计后的协方差阵即观测数据的期望 Fisher 信息。一种方式是用观测数据的观测 Fisher 信息$-{\ell }^{\prime \prime }(\hat{\theta }|x)$ 来近似, 其中 ${\ell }^{\prime \prime }$是 $logL(\theta |x)$ 的二阶导数，即 Hessian 矩阵。

Hessian矩阵是一个函数的二阶偏导数构成的方阵，它描述了这个函数在参数空间中的曲率和形状。

在有些情形, 这个 Hessian 阵可以解析计算出来，而在其他情形, 要得到或编码 Hessian 阵会很困难。在这些场合, 可用多种其他方法来简化协方差阵的估计。

（一）Louis’s Method

1、提出原因

Louis’s Method 是一种用于估计参数估计协方差矩阵的方法，特别是在应用EM（期望最大化）算法进行最大似然估计（MLE）时。由于EM算法的复杂性和涉及隐变量的结构，直接估计参数的协方差矩阵可能会很困难。Louis’s Method 提供了一种系统化的方法来从EM算法的输出中计算参数估计的协方差矩阵。

2、核心思想

• Louis’s Method的核心思想是通过计算Hessian矩阵的逆矩阵来估计参数的方差；
• Louis’s Method将完整数据的Hessian矩阵拆分为两部分：一部分由完整数据的似然函数得到，另一部分由观察数据的似然函数得到，用完全信息减去缺失信息来得到观测信息，观测信息的逆就是想要的协方差；
• Louis’s Method使用了EM算法中的一些特性，比如E步中的期望和M步中的最大化，来简化Hessian矩阵的计算。

已知：
$$log{f}_X(x|\theta )=log{f}_Y(y|\theta )-log{f}_{Z|X}(z|x;\theta )$$
等式两边对隐变量，在观测数据和参数$w=\theta$给定下，求条件期望即可得：
$$log{f}_X(x|\theta )=Q(\theta |w){|}_{w=\theta }-H(\theta |w){|}_{w=\theta }$$
再对第一个$\theta$求二阶偏导数：
$$-l^{\prime \prime }(\theta |x)=-Q^{\prime \prime }(\theta |w){|}_{w=\theta }+H^{\prime \prime }(\theta |w){|}_{w=\theta }$$

进而可以得到缺失信息法则：观测信息等于完全信息减去缺失信息
$${\hat{i}}_X(\theta )={\hat{i}}_Y(\theta )-{\hat{i}}_{Z|X}(\theta )$$

• i ^ Y ( θ ) = − Q ′ ′ ( θ ∣ w ) ∣ w = θ = − E { l ′ ′ ( θ ∣ Y ) ∣ x , θ } \hat{i}_Y(\theta)=-Q''(\theta|w)|_{w=\theta}=-E\{l''(\theta|Y)|x,\theta\} i^Y​(θ)=−Q′′(θ∣w)∣w=θ​=−E{l′′(θ∣Y)∣x,θ}
• i ^ Z ∣ X ( θ ) = v a r { d l o g f Z ∣ X ( Z ∣ x , θ ) d θ } \hat{i}_{Z|X}(\theta)=var\{\frac{d\ log\ f_{Z|X}(Z|x,\theta)}{d\theta}\} i^Z∣X​(θ)=var{dθd log fZ∣X​(Z∣x,θ)​}

缺失信息法则使得我们能够用完全数据似然和给定观测数据下缺失数据的条件密度来表达${\hat{i}}_X(\theta )$， 而且可以避免包括观测数据的可能复杂的边际似然的计算.。在某些情况下该方法可较容易得到并编码, 但它并不总比直接计算 $-\ell \prime \prime (\theta |x)$容易：如果 ${\hat{i}}_Y(\theta )$或者 ${\hat{i}}_{Z|X}(\theta )$难于解析计算。

此时可以通过 Monte Carlo 方法来估计：

• ${\hat{i}}_Y(\theta )$的最简单的 Monte Carlo 估计为
  $$\frac{1}{m}\sum _{i=1}^m-\frac{d^{2}log{f}_Y(y_i|\theta )}{d\theta \cdot d\theta }$$

3、计算步骤

• 计算EM算法的期望步骤（E步）和最大化步骤（M步），得到每次迭代的参数估计
• 计算每次迭代中的对数似然函数关于参数的一阶和二阶导数
• 使用这些导数计算Hessian矩阵
• 计算Hessian矩阵的逆矩阵
• 使用逆Hessian矩阵来估计参数的方差

（二）删失的指数数据

对于处理删失的指数数据，Louis’s Method可以被应用于估计参数的方差。删失数据是指在统计分析中，一部分观测数据因为某些原因没有被完全观察到。在这种情况下，EM算法非常适用，因为它可以处理含有未观测数据（隐变量）的参数估计问题。

使用Louis’s Method来处理删失的指数数据并估计参数方差：

步骤1：定义模型

假设我们有一个指数分布的数据集，某些数据点是删失的。指数分布的概率密度函数为：$f(t;\lambda )=\lambda e^{-\lambda t}$，$\lambda$是待估计的参数。

步骤2：EM算法的E步

在E步中，我们计算期望的对数似然函数，其中包含隐变量的条件期望。对于删失数据，我们需要计算每个观测值$y_i$的完全数据对数似然函数的期望。

完整数据的对数似然函数可以写成：
$$l(\lambda |y_1,...,y_n)=nlog\lambda -\lambda \sum _{i=1}^n{y}_i$$

步骤3：EM算法的M步

Q ( λ ∣ λ ( t ) ) = E { l ( λ ∣ Y 1 , . . . , Y n ) ∣ x , λ ( t ) } = n l o g λ − λ ∑ i = 1 n E { Y i ∣ x i , λ ( t ) } = n l o g λ − λ ∑ i = 1 n [ x i δ i + ( c i + 1 / λ ( t ) ) ( 1 − δ i ) ] = n l o g λ − λ ∑ i = 1 n [ x i δ i + c i ( 1 − δ i ) ] − C λ / λ ( t ) \begin{equation*} \begin{aligned} Q(\lambda|\lambda^{(t)}) &= E\{l(\lambda|Y_1,...,Y_n)|x,\lambda^{(t)}\} \\ &= nlog\lambda-\lambda\sum_{i=1}^nE\{Y_i|x_i,\lambda^{(t)}\} \\ &= nlog\lambda-\lambda\sum_{i=1}^n[x_i\delta_i+(c_i+1/\lambda^{(t)})(1-\delta_i)] \\ &= nlog\lambda-\lambda\sum_{i=1}^n[x_i\delta_i+c_i(1-\delta_i)]-C\lambda/\lambda^{(t)} \\ \end{aligned} \end{equation*} Q(λ∣λ(t))​=E{l(λ∣Y1​,...,Yn​)∣x,λ(t)}=nlogλ−λi=1∑n​E{Yi​∣xi​,λ(t)}=nlogλ−λi=1∑n​[xi​δi​+(ci​+1/λ(t))(1−δi​)]=nlogλ−λi=1∑n​[xi​δi​+ci​(1−δi​)]−Cλ/λ(t)​​

其中$C={\sum }_{i=1}^n(1-{\delta }_i)$表示删失事件的个数，则我们有：
$$-Q^{\prime \prime }(\lambda |{\lambda }^{(t)})=n/{\lambda }^2$$

对于未观测到的变量$Z_i$，即$y_i$超过$c_i$的部分，有密度 f Z i ∣ x ( z i ∣ x , λ ) = λ e x p { − λ ( z i − c i ) } 1 { z i > c i } f_{Z_i|x}(z_i|x,\lambda)=\lambda exp\{-\lambda(z_i-c_i)\} 1_{\{z_i>c_i\} } fZi​∣x​(zi​∣x,λ)=λexp{−λ(zi​−ci​)}1{zi​>ci​}​，则：
d l o g f Z ∣ X ( Z ∣ x , λ ) d λ = C / λ − ∑ { i : δ i = 0 } ( Z i − c i ) \frac{d\ log f_{Z|X}(Z|x,\lambda)}{d\lambda}=C/\lambda-\sum_{\{i:\delta_i=0\}}(Z_i-c_i) dλd logfZ∣X​(Z∣x,λ)​=C/λ−{i:δi​=0}∑​(Zi​−ci​)

进而得到缺失信息为：
$${\hat{i}}_{Z|X}(\lambda )=-E\frac{d^{2}log{f}_{Z|X}(Z|x,\lambda )}{d{\lambda }^2}=var\frac{dlog{f}_{Z|X}(Z|x,\lambda )}{d\lambda }=C/{\lambda }^2$$

应用Louis方法：
$${\hat{i}}_X(\lambda )=n/{\lambda }^2-C/{\lambda }^2=U/{\lambda }^2$$
其中$U={\sum }_{i=1}^n{\delta }_i$表示未删失事件的个数.。对这个例子, 通过直接分析容易验证 $-l^{\prime \prime }(\lambda |x)=U/{\lambda }^2$。

（三）自助法 (Bootstrap)

Bootstrap 是非常经典的求估计量方差的方法，不仅仅适用于 EM 算法的估计量。 对独立同分布的观测数据 x1, …, xn 来说：

• 有放回地从$x_1,...,x_n$抽取$x_1^{\ast },...,x_n^{\ast }$
• 通过EM算法在这组新的数据上得到新的估计${\hat{\theta }}_b$
• 重复上述过程，直到采用出B组参数 { θ ^ 1 , θ ^ 2 , . . . . , θ ^ B } \{\hat{\theta}_1,\hat{\theta}_2,....,\hat{\theta}_B\} {θ^1​,θ^2​,....,θ^B​}

假定 { θ ^ 1 , θ ^ 2 , . . . . , θ ^ B } \{\hat{\theta}_1,\hat{\theta}_2,....,\hat{\theta}_B\} {θ^1​,θ^2​,....,θ^B​}是从原本EM估计量$\hat{\theta }$的分布中采样出来的样本，则可以方便地画出估计量的分布图以及计算方差或分位数

如果使用参数自助法，则第一步改成用第二步估计出来的参数生成一组数据，而不是 resampling。

五、EM变型

EM 有时候表达式太过复杂，难以精确地计算期望或者求出最值，在这种情况下，可以选择牺牲一些计算速度或精确性做一些近似计算。

（一）改进 E 步 (MCEM)

E 步需要找到在观测数据条件下完全数据的期望对数似然。 我们已经用$Q(\theta |{\theta }^{(t)})$ 表示该期望，当该期望难以解析计算时, 可以用 Monte Carlo方法来近似。

1、Monte Carlo EM

Wei and Tanner 提出第 t 个 E 步可以用下面的两步替代：

• 从 $f_{Z|X}(z|x,{\theta }^{(t)})$ 中抽取独立同分布的缺失数据集 $Z_1^{(t)},...,Z_{m^{(t)}}^{(t)}$，每个$Z_j^{(t)}$是用来补齐观测数据集的所有缺失值的一个向量，这样$Y_j=(x,Z_j)$表示一个补齐的数据集，缺失值由$z_j$代替
• 计算${\hat{Q}}^{(t+1)}(\theta |{\theta }^{(t)})=\frac{1}{m^{(t)}}log{f}_Y(Y_j^{(t)}|\theta )$

此时${\hat{Q}}^{(t+1)}(\theta |{\theta }^{(t)})$就是$Q(\theta |{\theta }^{(t)})$的Monte Carlo估计，M步改为最大化${\hat{Q}}^{(t+1)}(\theta |{\theta }^{(t)})$。

2、删失的指数数

假定试图在模型$Y_1,...,Y_n\sim i.i.d.Exp(\lambda )$下观测到完全数据，但有些情形是右删失的。这样，观测数据是$x=(x_1,...,x_n)$，其中$x_i=(min(y_i,c_i),{\delta }_i)$，$c_i$是删失水平，如果$y_i\le c_i$，${\delta }_i=1$，否则${\delta }_i=0$。

完整数据的对数似然函数可以写成：
$$l(\lambda |y_1,...,y_n)=nlog\lambda -\lambda \sum _{i=1}^n{y}_i$$
Q函数为：
$$Q(\lambda |{\lambda }^{(t)})=nlog\lambda -\lambda \sum _{i=1}^n[x_i{\delta }_i+c_i(1-{\delta }_i)]-C\lambda /{\lambda }^{(t)}$$
所以我们的标准EM更新为：
$${\lambda }^{(t+1)}=\frac{n}{{\sum }_{i=1}^n[x_i{\delta }_i+c_i(1-{\delta }_i)]+C/{\lambda }^{(t)}}$$

而MCEM，将用MC的方法估计Q函数：
$${\hat{Q}}^{(t+1)}(\lambda |{\lambda }^{(t)})=nlog\lambda -\frac{\lambda }{m^{(t)}}\sum _{j=2}^{m^{(t)}}{Y}_j^{T}1$$

令$\widehat{Q^{\prime }}(\lambda |{\lambda }^{(t)})=0$且对$\lambda$求解得到：
$${\lambda }^{(t+1)}=\frac{n}{{\sum }_{i=1}^{(t)}{Y}_j^{T}1/m^{(t)}}$$
作为 MCEM 的更新。

3、MCEM 和 EM 对比

• 计算复杂度：EM算法在E步中计算复杂度较高，但如果可以解析求解，则效率较高。MCEM算法则使用采样方法，计算复杂度通常更高。
• 适用范围：EM算法适用于期望值可以解析计算的情况，而MCEM算法则适用于无法解析计算期望值的复杂模型。
• 收敛性：两者都可以收敛，但MCEM的收敛速度和精度受采样质量影响更大。EM算法的收敛性更稳定，但可能会陷入局部极大值。

（二）改进 M 步

EM 算法的吸引力之一在于$Q(\theta |{\theta }^{(t)})$ 的求导和最大化通常比不完全数据极大似然的计算简单，这是因为$Q(\theta |{\theta }^{(t)})$与完全数据似然有关。然而，在某些情况下，即使导出$Q(\theta |{\theta }^{(t)})$的 E 步是直接了当的，M 步也不容易实施. 为此人们提出了多种策略以便于 M 步的实施。

1、ECM算法

Meng 和 Rubin 的 ECM 算法是用一系列计算较简单的条件极大化步骤（CM步）代替 M 步，每次条件极大化均被设计为一个简单的优化问题， 该优化问题把 $\theta$ 限制在某特殊子空间, 使得可以得到解析解或非常初等的数值解。

• 第t个E步：与EM算法相同，计算当前参数下隐藏变量的条件期望
  Q ( θ ∣ θ ( t ) ) = E { l o g L ( θ ; X , Z ) ∣ x , θ ( t ) } Q(\theta|\theta^{(t)}) = E\{log L(\theta;X,Z)|x,\theta^{(t)}\} Q(θ∣θ(t))=E{logL(θ;X,Z)∣x,θ(t)}
• 第t个CM循环：将M步分解为多个条件最大化步骤。在每个CM步中，只对一部分参数进行最大化，而保持其他参数固定。重复这个过程直到所有参数都更新完毕。令$S$表示每个CM循环中的CM步的数目，第t次循环里第s个CM步需要在约束：
  $$g_s(\theta )=g_s({\theta }^{(t+(s-1)/S)})$$
  下最大化$Q(\theta |{\theta }^{(t)})$，其中${\theta }^{(t+(s-1)/S}$是在当前循环的第$(s-1)$个CM步中求得的极大值点。当S个CM步的整个循环完成时，令${\theta }^{(t+1)}={\theta }^{(t+S/S)}$并进行第$(t+1)$次迭代的E步。

构造有效 ECM 算法的技巧在于巧妙地选择约束条件。通常，可自然地把 $\theta$ 分成 $S$ 个子向量 $\theta =({\theta }_1,...,{\theta }_S)$， 然后在第 $s$ 个 CM 步中，固定 $\theta$ 其余的元素而关于 ${\theta }_s$ 寻求最大化 $Q$：
$$g_s(\theta )=({\theta }_1,...,{\theta }_{s-1},{\theta }_{s+1},...,{\theta }_S)$$

另外，第 $s$ 个 CM 步也可以在固定 ${\theta }_s$ 下关于 $\theta$ 的其他元素最大化 $Q$：
$$g_s(\theta )={\theta }_s$$
也可根据特定的问题背景想象其他的约束体系。ECM 的一种变型是在每两个 CM 步之间插入一个 E 步，由此在 CM循环的每一个阶段均更新了 $Q$。

2、EM梯度算法

如果最大化不能用解析的方法来实现，那么可以考虑采用迭代优化方法来实施每个 M 步。 这将会产生一个有嵌套迭代循环的算法。ECM 算法在 EM 算法的每次迭代中搬入 S 个条件最大化步骤， 这也会产生嵌套迭代。

为避免嵌套循环的计算负担，Lange 提出用单步 Newton 法替代 M 步，从而可近似取得最大值而不用真正地精确求解。M 步是用由

$$\begin{array}{c}\begin{array}{rl}{\theta }^{(t+1)}& ={\theta }^{(t)}-Q^{\prime \prime }(\theta |{\theta }^{(t)}{)}^{-1}{|}_{\theta ={\theta }^{(t)}}{Q}^{\prime }(\theta |{\theta }^{(t)}){|}_{\theta ={\theta }^{(t)}}\\ & ={\theta }^{(t)}-Q^{\prime \prime }(\theta |{\theta }^{(t)}{)}^{-1}{|}_{\theta ={\theta }^{(t)}}{l}^{\prime }({\theta }^{(t)}|x)\end{array}\end{array}$$
给出的更新替代的，其中$l^{\prime }({\theta }^{(t)}|x)$是当前迭代得分函数的估值，第二个等式是由${\theta }^{(t)}$最大化$Q(\theta |{\theta }^{(t)})-l(\theta |x)$的结论得到。

也可以用最速梯度法代替Newton法，从而避免计算$Q$函数矩阵的逆矩阵：
$$\begin{array}{c}\begin{array}{rl}{\theta }^{(t+1)}& ={\theta }^{(t)}-{\alpha }^{(t)}{Q}^{\prime }(\theta |{\theta }^{(t)}){|}_{\theta ={\theta }^{(t)}}\\ & ={\theta }^{(t)}-{\alpha }^{(t)}{l}^{\prime }({\theta }^{(t)}|x)\end{array}\end{array}$$

这篇关于统计计算二|EM算法（Expectation-Maximization Algorithm，期望最大化算法）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---