第一部分：贝叶斯网基础

1.1 信息论基础

1.2 贝叶斯网基本概念

1.3 变量独立性的图论分析

第二部分：贝叶斯网推理

2.1 概率推理中的变量消元方法

2.2 团树传播算法

2.3 近似推理

2.3.1 蒙特卡洛方法

2.3.1.1 重要性抽样法

2.3.1.2 马尔可夫蒙特卡洛抽样法(MCMC)

2.3.2 变分推理

第三部分：贝叶斯网参数学习

3.1 理论基础-分布的分布

3.2 极大似然估计

贝叶斯网参数估计有两种方法，一种是极大似然估计，另一种是贝叶斯估计。两种参数估计方法使用不同的思想，前者来自于频率派，认为参数是固定的，我们要做的事情就是根据已经掌握的数据来估计这个参数；而后者属于贝叶斯派，认为参数服从某种概率分布，已有的数据只是在这种参数的分布下产生的。所以，极大似然估计就是假设一个参数 θ \theta θ，然后根据观测到的样本数据来求出这个 θ \theta θ；而贝叶斯估计则需要先设定参数的先验概率 p ( θ ) p(\theta) p(θ)，然后再结合观测到的样本数据的似然度和贝叶斯公式计算后验概率，并以 θ \theta θ在该后验概率分布下的期望来估计 θ \theta θ。
所以极大似然估计与贝叶斯估计最大的不同就在于是否考虑了先验，而两者适用范围也变成了：极大似然估计适用于样本量足够大，可以用样本集估计整体的情况；而贝叶斯估计则适用于对待估参数已有先验知识，只需通过较少的样本量来修正先验知识的情况。
本节将先介绍极大似然估计（Maximum Likelihood Estimation, MLE）。贝叶斯估计将在下一节介绍。
MLE是经典算法，大部分人在初高中时就已经接触过。那是通过MLE估计单个参数的最优取值。我们先通过这种最简单情况下的MLE来回顾其概念。然后再逐步深入讨论MLE应用于贝叶斯网参数的估计。主要按如下思路展开：

• MLE回顾：单变量贝叶斯网参数估计
• MLE推广：多变量贝叶斯网参数估计
• 缺值数据MLE：EM算法

3.2.1 MLE回顾：单变量贝叶斯网参数估计

3.2.1.1 二项分布案例引入

我们从一个最简单的经典例子开始。假设有一个只含单个变量 X X X的简单贝叶斯网， X X X表示投掷图钉的结果，它有两个可能的取值：帽朝上(head, h)和针朝上(tail, t)。此贝叶斯网只有一个独立参数，即帽朝上的概率 θ = P ( X = h ) \theta=P(X=h) θ=P(X=h)，针朝上的概率与之互补，为 1 − θ 1-\theta 1−θ。
这里之所以选择投掷图钉的实验，而不是选择更常见的抛硬币实验，是因为对于抛硬币事件，我们有更充分的先验知识，即硬币正反的概率相等，而投掷图钉的结果对我们而言则没有充分的先验知识。因此估算投掷图钉事件的参数我们只能完全依赖于所观测到的数据，即采用MLE顺理成章。而对于抛硬币事件的参数估计，则更适合采用贝叶斯估计。
假设我们投掷了6次图钉，结果依次是 t , h , t , t , h , t t,h,t,t,h,t t,h,t,t,h,t。基于这个结果，直观上可以用图钉帽朝上出现的频率 2 / 6 2/6 2/6作为对参数 θ \theta θ的估计。
该直观思想背后的原则是这样的：假如我们考虑 θ = 0 \theta=0 θ=0的情况，这个取值与观测到的数据 D = ( t , h , t , t , h , t ) D=(t,h,t,t,h,t) D=(t,h,t,t,h,t)明显是矛盾的，相比之下， θ = 2 / 6 \theta=2/6 θ=2/6与数据更相容。所以，如果要在0和 2 / 6 2/6 2/6之间选择一个对 θ \theta θ的估计，应该选择 2 / 6 2/6 2/6。这意味着，在判断某个值作为 θ \theta θ的估计是否合适时，应该考虑它与数据的拟合程度。
一般而言，参数对数据的拟合程度可以通过条件概率 P ( D ∣ θ = θ 0 ) P(D|\theta=\theta_0) P(D∣θ=θ0​)来判断，这个概率越大，说明 θ 0 \theta_0 θ0​对数据 D D D的拟合程度越高。在上面的例子中， P ( D ∣ θ = 0 ) = 0 , P ( D ∣ θ = 2 / 6 ) > 0 P(D|\theta=0)=0,P(D|\theta=2/6)>0 P(D∣θ=0)=0,P(D∣θ=2/6)>0，所以后者对数据的拟合程度更高，更适合作为参数 θ \theta θ的估计。
下面给出几个概念。

似然度（likelihood）: 给定参数 θ \theta θ，数据 D D D的条件概率 P ( D ∣ θ ) P(D|\theta) P(D∣θ)称为 θ \theta θ的似然度，记为：
L ( θ ∣ D ) = P ( D ∣ θ ) L(\theta|D)=P(D|\theta) L(θ∣D)=P(D∣θ)
似然函数（likelihood function）: 当数据给定后（ D D D固定），让 θ \theta θ在其定义域上变动， L ( θ ∣ D ) L(\theta|D) L(θ∣D)就是 θ \theta θ的函数，称为 θ \theta θ的似然函数。
最大似然估计（Maximum Likelihood Estimation, MLE）： 用使得 L ( θ ∣ D ) L(\theta|D) L(θ∣D)达到最大的那个取值 θ ∗ \theta^* θ∗估计 θ \theta θ，称为 θ \theta θ的最大似然估计，即：
θ ∗ = arg max ⁡ θ L ( θ ∣ D ) \theta^*=\argmax_\theta L(\theta|D) θ∗=θargmax​L(θ∣D)

显然，MLE的概念可以推广到多变量和多参数的情况。

3.2.1.2 MLE的计算

数据 D D D由样本 ( D 1 , D 2 , ⋯   , D m ) (D_1,D_2,\cdots,D_m) (D1​,D2​,⋯,Dm​)组成，这里有关于数据的两个基本假设：
a）给定参数 θ \theta θ时， D D D中样本相互独立，即：
L ( θ ∣ D ) = P ( D ∣ θ ) = ∏ i = 1 m P ( D i ∣ θ ) L(\theta|D)=P(D|\theta)=\prod_{i=1}^mP(D_i|\theta) L(θ∣D)=P(D∣θ)=i=1∏m​P(Di​∣θ)
b）每个样本 D i D_i Di​的条件概率分布 P ( D i ∣ θ ) P(D_i|\theta) P(Di​∣θ)相同。在投掷图钉的例子中，即设对每一个 i i i都有：
P ( D i = h ∣ θ ) = θ , P ( D i = t ∣ θ ) = 1 − θ P(D_i=h|\theta)=\theta,P(D_i=t|\theta)=1-\theta P(Di​=h∣θ)=θ,P(Di​=t∣θ)=1−θ
这两个假设是统计学中的基本假设，称为独立同分布假设（Independent and Identically Distributed），简称i.i.d.
在投掷图钉的例子中，根据i.i.d.假设，似然函数为：
L ( θ ∣ D ) = θ m h ( 1 − θ ) m t (1) L(\theta|D)=\theta^{m_h}(1-\theta)^{m_t} \tag{1} L(θ∣D)=θmh​(1−θ)mt​(1)
其中 m h m_h mh​和 m t m_t mt​分别是数据集 D D D中图钉帽朝上和针朝上出现的次数。具有式(1)形式的似然函数称为二项似然函数(binomial likelihood function).
在投掷图钉的例子中， D = ( t , h , t , t , h , t ) D=(t,h,t,t,h,t) D=(t,h,t,t,h,t)，从而有：
L ( θ ∣ D ) = P ( D ∣ θ ) = P ( D 1 = t ∣ θ ) P ( D 2 = h ∣ θ ) ⋯ P ( D 6 = t ∣ θ ) = θ 2 ( 1 − θ ) 4 \begin{aligned} L(\theta|D)&=P(D|\theta)\\ &=P(D_1=t|\theta)P(D_2=h|\theta)\cdots P(D_6=t|\theta)\\ &=\theta^2(1-\theta)^4 \end{aligned} L(θ∣D)​=P(D∣θ)=P(D1​=t∣θ)P(D2​=h∣θ)⋯P(D6​=t∣θ)=θ2(1−θ)4​
对似然函数 L ( θ ∣ D ) L(\theta|D) L(θ∣D)取对数，就得到对数似然函数函数(loglikelihood function)，即：
l ( θ ∣ D ) = log ⁡ L ( θ ∣ D ) (2) l(\theta|D)=\log L(\theta|D) \tag{2} l(θ∣D)=logL(θ∣D)(2)
由于 log ⁡ \log log函数的单调性， l ( θ ∣ D ) l(\theta|D) l(θ∣D)的最大值点也是 L ( θ ∣ D ) L(\theta|D) L(θ∣D)的最大值点，但前者通过取对数将连乘变为累加，计算导数时更方便。
在投掷图钉的例子中， l ( θ ∣ D ) = m h log ⁡ θ + m t log ⁡ ( 1 − θ ) l(\theta|D)=m_h\log\theta+m_t\log(1-\theta) l(θ∣D)=mh​logθ+mt​log(1−θ)，其导数为0的点即为最大值点，从而可求出：
θ ∗ = m h m h + m t = m h m (3) \theta^*=\frac{m_h}{m_h+m_t}=\frac{m_h}{m} \tag{3} θ∗=mh​+mt​mh​​=mmh​​(3)
由此可见，在上一小节中我们给出的直观估计是合理的，参数 θ \theta θ的最大似然估计就是数据中图钉帽朝上出现的频率。

3.2.1.3 多项分布推广

前面介绍了仅取二值的单变量贝叶斯网的单参数极大似然估计，本节将推广到可取多值的单变量贝叶斯网的多参数极大似然估计。
考虑由一个多值变量 X X X组成的贝叶斯网。设 X X X有 r r r个取值， Ω X = { x 1 , x 2 , ⋯   , x r } \Omega_X=\{x_1,x_2,\cdots,x_r\} ΩX​={x1​,x2​,⋯,xr​}。网络的参数包括 θ i = P ( X = x i ) , i = 1 , 2 , ⋯   , r \theta_i=P(X=x_i),i=1,2,\cdots,r θi​=P(X=xi​),i=1,2,⋯,r。用 θ \mathbb\theta θ记向量 ( θ 1 , ⋯   , θ r ) (\theta_1,\cdots,\theta_r) (θ1​,⋯,θr​)。由于概率归一性，实际上网络只有 r − 1 r-1 r−1个独立参数。
设有一组i.i.d.数据 D = ( D 1 , ⋯   , D m ) D=(D_1,\cdots,D_m) D=(D1​,⋯,Dm​)，其中满足 X = x i X=x_i X=xi​的样本个数是 m i m_i mi​，则有：
L ( θ ∣ D ) = ∏ i = 1 r θ i m i (4) L(\theta|D)=\prod_{i=1}^r\theta_i^{m_i} \tag{4} L(θ∣D)=i=1∏r​θimi​​(4)
具有这个形式的似然函数称为多项似然函数(multinomial likelihood function)。相应的对数似然函数为：
l ( θ ∣ D ) = ∑ i = 1 r m i log ⁡ θ i (5) l(\theta|D)=\sum_{i=1}^rm_i\log\theta_i \tag{5} l(θ∣D)=i=1∑r​mi​logθi​(5)
通过拉格朗日乘数法引入概率归一化条件，并求偏导数为0的点即为最大值点，从而可求出 θ \theta θ的最大似然估计：
θ ∗ = m i m (6) \theta^*=\frac{m_i}{m} \tag{6} θ∗=mmi​​(6)

3.2.2 MLE推广：多变量贝叶斯网参数估计

考虑一个由 n n n个变量 X = { X 1 , X 2 , ⋯   , X n } X=\{X_1,X_2,\cdots,X_n\} X={X1​,X2​,⋯,Xn​}组成的贝叶斯网 N N N。不失一般性，设其中的节点 X i X_i Xi​共有 r i r_i ri​个取值 1 , 2 , ⋯   , r i 1,2,\cdots,r_i 1,2,⋯,ri​，其父节点 π ( X i ) \pi(X_i) π(Xi​)的取值共有 q i q_i qi​种组合 1 , 2 , ⋯   , q i 1,2,\cdots,q_i 1,2,⋯,qi​。若 X i X_i Xi​无父节点，则 q i = 1 q_i=1 qi​=1。那么，网络的参数为：
θ i j k = P ( X i = k ∣ π ( X i ) = j ) \theta_{ijk}=P(X_i=k|\pi(X_i)=j) θijk​=P(Xi​=k∣π(Xi​)=j)
其中 i i i的取值范围是 1 , 2 , ⋯   , n 1,2,\cdots,n 1,2,⋯,n，表示 n n n个变量中的第 i i i个。对于第 i i i个变量 X i X_i Xi​， j j j的取值范围是 1 , 2 , ⋯   , q i 1,2,\cdots,q_i 1,2,⋯,qi​，表示其父节点 π ( X i ) \pi(X_i) π(Xi​)的第 j j j种取值组合； k k k的取值范围是 1 , 2 , ⋯   , r i 1,2,\cdots,r_i 1,2,⋯,ri​，表示 X i X_i Xi​的第 k k k种取值。用 θ \bm\theta θ记所有 θ i j k \theta_{ijk} θijk​组成的向量，由于概率归一性，该贝叶斯网 N N N的独立参数个数为 ∑ i = 1 n q i ( r i − 1 ) \sum_{i=1}^nq_i(r_i-1) ∑i=1n​qi​(ri​−1)。

下面给出一个具体的示例。如上图所示的贝叶斯网，所有变量都取两个值：1或2。将变量父节点的取值组合按字典序排列，并从1开始编号。例如 X 3 X_3 X3​有两个父节点 X 1 , X 2 X_1,X_2 X1​,X2​，它们的取值共有4种组合，排序为 ( X 1 = 1 , X 2 = 1 ) , ( X 1 = 1 , X 2 = 2 ) , ( X 1 = 2 , X 2 = 1 ) , ( X 1 = 2 , X 2 = 2 ) (X_1=1,X_2=1),(X_1=1,X_2=2),(X_1=2,X_2=1),(X_1=2,X_2=2) (X1​=1,X2​=1),(X1​=1,X2​=2),(X1​=2,X2​=1),(X1​=2,X2​=2)，编号分别为 1 , 2 , 3 , 4 1,2,3,4 1,2,3,4，则参数 θ \bm\theta θ如下：
θ 111 = P ( X 1 = 1 ) , θ 112 = P ( X 1 = 2 ) , θ 211 = P ( X 2 = 1 ) , θ 212 = P ( X 2 = 2 ) , θ 311 = P ( X 3 = 1 ∣ X 1 = 1 , X 2 = 1 ) , θ 312 = P ( X 3 = 2 ∣ X 1 = 1 , X 2 = 1 ) , θ 321 = P ( X 3 = 1 ∣ X 1 = 1 , X 2 = 2 ) , θ 322 = P ( X 3 = 2 ∣ X 1 = 1 , X 2 = 2 ) , θ 331 = P ( X 3 = 1 ∣ X 1 = 2 , X 2 = 1 ) , θ 332 = P ( X 3 = 2 ∣ X 1 = 2 , X 2 = 1 ) , θ 341 = P ( X 3 = 1 ∣ X 1 = 2 , X 2 = 2 ) , θ 342 = P ( X 3 = 2 ∣ X 1 = 2 , X 2 = 2 ) , ⋯ \begin{aligned} &\theta_{111}=P(X_1=1),&\theta_{112}&=P(X_1=2),\\ &\theta_{211}=P(X_2=1),&\theta_{212}&=P(X_2=2),\\ &\theta_{311}=P(X_3=1|X_1=1,X_2=1),&\theta_{312}&=P(X_3=2|X_1=1,X_2=1),\\ &\theta_{321}=P(X_3=1|X_1=1,X_2=2),&\theta_{322}&=P(X_3=2|X_1=1,X_2=2),\\ &\theta_{331}=P(X_3=1|X_1=2,X_2=1),&\theta_{332}&=P(X_3=2|X_1=2,X_2=1),\\ &\theta_{341}=P(X_3=1|X_1=2,X_2=2),&\theta_{342}&=P(X_3=2|X_1=2,X_2=2),\\ &\cdots \end{aligned} ​θ111​=P(X1​=1),θ211​=P(X2​=1),θ311​=P(X3​=1∣X1​=1,X2​=1),θ321​=P(X3​=1∣X1​=1,X2​=2),θ331​=P(X3​=1∣X1​=2,X2​=1),θ341​=P(X3​=1∣X1​=2,X2​=2),⋯​θ112​θ212​θ312​θ322​θ332​θ342​​=P(X1​=2),=P(X2​=2),=P(X3​=2∣X1​=1,X2​=1),=P(X3​=2∣X1​=1,X2​=2),=P(X3​=2∣X1​=2,X2​=1),=P(X3​=2∣X1​=2,X2​=2),​
设 D = ( D 1 , ⋯   , D m ) D=(D_1,\cdots,D_m) D=(D1​,⋯,Dm​)是一组关于 N N N的i.i.d.数据，则 θ \bm\theta θ的对数似然函数为：
l ( θ ∣ D ) = log ⁡ ∏ l = 1 m P ( D l ∣ θ ) = ∑ l = 1 m log ⁡ P ( D l ∣ θ ) l(\bm\theta|D)=\log\prod_{l=1}^mP(D_l|\bm\theta)=\sum_{l=1}^m\log P(D_l|\bm\theta) l(θ∣D)=logl=1∏m​P(Dl​∣θ)=l=1∑m​logP(Dl​∣θ)
为了得到关于 log ⁡ P ( D l ∣ θ ) \log P(D_l|\bm\theta) logP(Dl​∣θ)的表达式，定义样本 D l D_l Dl​的特征函数 χ ( i , j , k : D l ) \chi(i,j,k:D_l) χ(i,j,k:Dl​)如下：
χ ( i , j , k : D l ) = { 1 , 若 在 D l 中 X i = k 且 π ( X i ) = j 0 , 若 否 \chi(i,j,k:D_l)=\left\{ \begin{aligned} &1,\quad 若在D_l中X_i=k且\pi(X_i)=j\\ &0,\quad 若否 \end{aligned} \right. χ(i,j,k:Dl​)={​1,若在Dl​中Xi​=k且π(Xi​)=j0,若否​
从而有：
log ⁡ P ( D l ∣ θ ) = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i χ ( i , j , k : D l ) log ⁡ θ i j k (7) \log P(D_l|\theta)=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\chi(i,j,k:D_l)\log\theta_{ijk} \tag{7} logP(Dl​∣θ)=i=1∑n​j=1∑qi​​k=1∑ri​​χ(i,j,k:Dl​)logθijk​(7)
比如，在上例的贝叶斯网中，样本 D 1 = ( 1 , 1 , 2 , 2 , 1 ) D_1=(1,1,2,2,1) D1​=(1,1,2,2,1)，其对应的特征函数为：
χ ( 1 , 1 , 1 : D 1 ) = 1 χ ( 2 , 1 , 1 : D 1 ) = 1 χ ( 3 , 1 , 2 : D 1 ) = 1 χ ( 4 , 1 , 2 : D 1 ) = 1 χ ( 5 , 4 , 1 : D 1 ) = 1 其 它 = 0 \chi(1,1,1:D_1)=1\\ \chi(2,1,1:D_1)=1\\ \chi(3,1,2:D_1)=1\\ \chi(4,1,2:D_1)=1\\ \chi(5,4,1:D_1)=1\\ 其它=0 χ(1,1,1:D1​)=1χ(2,1,1:D1​)=1χ(3,1,2:D1​)=1χ(4,1,2:D1​)=1χ(5,4,1:D1​)=1其它=0
我们按照贝叶斯网的链式法则对 P ( D l ∣ θ ) P(D_l|\bm\theta) P(Dl​∣θ)进行分解：
log ⁡ P ( D l ∣ θ ) = log ⁡ P ( X 1 = 1 , X 2 = 1 , X 3 = 2 , X 4 = 2 , X 5 = 1 ∣ θ ) = log ⁡ ( P ( X 1 = 1 ∣ θ ) P ( X 2 = 1 ∣ θ ) P ( X 3 = 2 ∣ X 1 = 1 , X 2 = 1 , θ ) P ( X 4 = 2 ∣ X 1 = 1 , X 2 = 1 , θ ) P ( X 5 = 1 ∣ X 3 = 2 , X 4 = 2 , θ ) ) = log ⁡ ( θ 111 θ 211 θ 312 θ 412 θ 541 ) = log ⁡ θ 111 + log ⁡ θ 211 + log ⁡ θ 312 + log ⁡ θ 412 + log ⁡ θ 541 = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i χ ( i , j , k : D l ) log ⁡ θ i j k \begin{aligned} \log P(D_l|\theta)&=\log P(X_1=1,X_2=1,X_3=2,X_4=2,X_5=1|\bm\theta)\\ &=\log(P(X_1=1|\bm\theta)P(X_2=1|\bm\theta)P(X_3=2|X_1=1,X_2=1,\bm\theta)P(X_4=2|X_1=1,X_2=1,\bm\theta)P(X_5=1|X_3=2,X_4=2,\bm\theta))\\ &=\log(\theta_{111}\theta_{211}\theta_{312}\theta_{412}\theta_{541})\\ &=\log\theta_{111}+\log\theta_{211}+\log\theta_{312}+\log\theta_{412}+\log\theta_{541}\\ &=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\chi(i,j,k:D_l)\log\theta_{ijk} \end{aligned} logP(Dl​∣θ)​=logP(X1​=1,X2​=1,X3​=2,X4​=2,X5​=1∣θ)=log(P(X1​=1∣θ)P(X2​=1∣θ)P(X3​=2∣X1​=1,X2​=1,θ)P(X4​=2∣X1​=1,X2​=1,θ)P(X5​=1∣X3​=2,X4​=2,θ))=log(θ111​θ211​θ312​θ412​θ541​)=logθ111​+logθ211​+logθ312​+logθ412​+logθ541​=i=1∑n​j=1∑qi​​k=1∑ri​​χ(i,j,k:Dl​)logθijk​​
从而验证了上式(7)。
我们定义：
m i j k = ∑ l = 1 m χ ( i , j , k : D l ) m_{ijk}=\sum_{l=1}^m\chi(i,j,k:D_l) mijk​=l=1∑m​χ(i,j,k:Dl​)
从而， l ( θ ∣ D ) l(\bm\theta|D) l(θ∣D)的表达式可进一步表示为：
l ( θ ∣ D ) = log ⁡ ∏ l = 1 m P ( D l ∣ θ ) = ∑ l = 1 m log ⁡ P ( D l ∣ θ ) = ∑ l = 1 m ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i χ ( i , j , k : D l ) log ⁡ θ i j k = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i ( ∑ l = 1 m χ ( i , j , k : D l ) ) log ⁡ θ i j k = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i m i j k log ⁡ θ i j k (8) \begin{aligned} l(\bm\theta|D)&=\log\prod_{l=1}^mP(D_l|\bm\theta)=\sum_{l=1}^m\log P(D_l|\bm\theta)\\ &=\sum_{l=1}^m\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\chi(i,j,k:D_l)\log\theta_{ijk}\\ &=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\left(\sum_{l=1}^m\chi(i,j,k:D_l)\right)\log\theta_{ijk}\\ &=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}m_{ijk}\log\theta_{ijk} \end{aligned}\tag{8} l(θ∣D)​=logl=1∏m​P(Dl​∣θ)=l=1∑m​logP(Dl​∣θ)=l=1∑m​i=1∑n​j=1∑qi​​k=1∑ri​​χ(i,j,k:Dl​)logθijk​=i=1∑n​j=1∑qi​​k=1∑ri​​(l=1∑m​χ(i,j,k:Dl​))logθijk​=i=1∑n​j=1∑qi​​k=1∑ri​​mijk​logθijk​​(8)
对任意 i , j i,j i,j，存在概率归一化条件 ∑ k = 1 r i θ i j k = 1 \sum_{k=1}^{r_i}\theta_{ijk}=1 ∑k=1ri​​θijk​=1，通过拉格朗日乘数法引入概率归一化条件，并求偏导为0的点即为最大值点，从而可求出 θ \bm\theta θ的极大似然估计：
θ i j k ∗ = { m i j k ∑ k = 1 r i m i j k , 若 ∑ k = 1 r i m i j k > 0 1 r i , 若 否 (9) \theta_{ijk}^*=\left\{ \begin{aligned} &\frac{m_{ijk}}{\sum_{k=1}^{r_i}m_{ijk}},&\quad &若\sum_{k=1}^{r_i}m_{ijk}>0\\ &\frac{1}{r_i},&\quad &若否 \end{aligned} \right.\tag{9} θijk∗​=⎩⎪⎪⎪⎨⎪⎪⎪⎧​​∑k=1ri​​mijk​mijk​​,ri​1​,​​若k=1∑ri​​mijk​>0若否​(9)
直观上理解上式，即：
θ i j k ∗ = D 中 满 足 X i = k 和 π ( X i ) = j 的 样 本 数 目 D 中 满 足 π ( X i ) = j 的 样 本 数 目 \theta_{ijk}^*=\frac{D中满足X_i=k和\pi(X_i)=j的样本数目}{D中满足\pi(X_i)=j的样本数目} θijk∗​=D中满足π(Xi​)=j的样本数目D中满足Xi​=k和π(Xi​)=j的样本数目​
当D中不存在满足 π ( X i ) = j \pi(X_i)=j π(Xi​)=j的样本时，则令 θ i j k ∗ \theta_{ijk}^* θijk∗​服从均匀分布。

下面我们通过一个简单的例子来演示贝叶斯网参数的极大似然估计过程。
考虑如下图(a)所示的贝叶斯网 N N N，其中所有变量均取二值1或2。右侧表格(b)为关于 N N N的一组i.i.d.数据。根据式(8)， N N N的参数的MLE估计分别如下表格©所示。

3.2.3 缺值数据MLE：EM算法

3.2.3.1 缺值的类别

有很多原因可导致我们所获得的样本集中的数据有缺项，数据缺项可分为以下三类：

• 完全随机缺失（Missing Completely At Random, MCAR）：某一变量缺失值不依赖于其他任何原因的完全随机缺失。比如在电子显微镜中，随机出现的灰尘会遮挡某些区域，使得该区域的观测值缺失；
• 随机缺失（Missing At Random, MAR)：某一缺值变量的实际取值与其他变量相关，且服从某种分布。比如在问卷调查中，如果年龄选项为未成年，则身高不填。若在所给出的样本集中，某个变量全部未观测到，则这样的缺值变量可称为隐变量，关于隐变量的分析在以后的专题中再单独展开，其中包含了很多迄今为止仍未解决的难题。
• 非随机缺失（Missing Not At Random, MNAR)：某一变量的缺失和该变量本身的数值相关。比如仪器的最低检测线，某被检测物质的含量低于该检测线则会产生非随机缺失。

对于不同类型的缺值，应该采取不同的处理方式。
针对MCAR和MAR，可通过某种精心设计的插值方法补齐缺值，或者在不影响样本集质量的情况下直接删去缺值样本。
针对MNAR，需要使用特殊值填补缺值，比如0，或者干脆重新进行测量（仪器误差太大）。
EM算法则是一种通过迭代插值进行优化逼近的算法，可用于处理具有MCAR/MAR缺值的样本集。

3.2.3.2 EM算法基本思想

在完整数据情况下，MLE可以用式(9)计算。但在现实应用中往往会遇到缺值的样本集。当数据有缺值时， log ⁡ P ( D l ∣ θ ) \log P(D_l|\theta) logP(Dl​∣θ)不能按式(7)分解，对数似然函数 l ( θ ∣ D ) l(\theta|D) l(θ∣D)也不能写成式(8)的形式。在实际中，人们通常使用期望最大（Exception Maximization, EM）算法进行逼近。
为了求得贝叶斯网 N = ( G , θ ) N=(G,\theta) N=(G,θ)的参数 θ \theta θ的MLE近似值，EM算法从 θ \theta θ的某个初值 θ 0 \theta^0 θ0出发开始迭代，初始值 θ 0 \theta^0 θ0可通过随机产生。设已经进行了 t t t次迭代，得到估计 θ t \theta^t θt。第 t + 1 t+1 t+1次迭代由以下E步和M步两个步骤组成：
（1）E步：基于 θ t \theta^t θt计算缺值数据的期望，并利用该期望值修补数据，使之完整；
（2）M步：基于修补后的完整数据计算 θ \theta θ的最大似然估计，得到 θ t + 1 \theta^{t+1} θt+1。
重复以上步骤，直至参数收敛。
EM算法最核心的两步分别是计算期望（E步）和计算最大似然估计（M步），这正是EM算法名字的由来。下面对这两步的计算方法展开介绍。
对于某一缺值样本 D l D_l Dl​，设 X l X_l Xl​为 D l D_l Dl​中缺值变量的集合，将 X l X_l Xl​的一个可能取值 x l x_l xl​加入 D l D_l Dl​，就得到一个完整样本 ( D l , X l = x l ) (D_l,X_l=x_l) (Dl​,Xl​=xl​)，这个过程称为数据修补。EM算法对 X l X_l Xl​的所有可能取值都进行考虑，并给每一种可能关联一个权重 ω x l \omega_{x_l} ωxl​​，得到加权样本 ( D l , X l = x l ) [ ω x l ] (D_l,X_l=x_l)[\omega_{x_l}] (Dl​,Xl​=xl​)[ωxl​​]。该权重由条件概率给定： ω x l = P ( X l = x l ∣ D l , θ t ) \omega_{x_l}=P(X_l=x_l|D_l,\theta^t) ωxl​​=P(Xl​=xl​∣Dl​,θt)。由于 ∑ x l ∈ Ω X l ω x l = 1 \sum_{x_l\in\Omega_{X_l}}\omega_{x_l}=1 ∑xl​∈ΩXl​​​ωxl​​=1，因此加权样本 ( D l , X l = x l ) [ ω x l ] (D_l,X_l=x_l)[\omega_{x_l}] (Dl​,Xl​=xl​)[ωxl​​]又称为碎权样本（fractional sample）。原本完整的样本无需修补，权重为1。
在数据修补过程中，每个缺值样本都被一系列碎权样本替代，因此修补后的样本集是完整的，且每个样本都有一个权重。对于完整样本集，利用式(9)即可计算出基于补后数据的极大似然估计 θ t + 1 \bm\theta^{t+1} θt+1，即：
θ i j k t + 1 = { m i j k t ∑ k = 1 r i m i j k t , 若 ∑ k = 1 r i m i j k t > 0 1 r i , 若 否 (10) \theta_{ijk}^{t+1}=\left\{ \begin{aligned} &\frac{m^t_{ijk}}{\sum_{k=1}^{r_i}m^t_{ijk}},&\quad &若\sum_{k=1}^{r_i}m^t_{ijk}>0\\ &\frac{1}{r_i},&\quad &若否 \end{aligned} \right.\tag{10} θijkt+1​=⎩⎪⎪⎪⎨⎪⎪⎪⎧​​∑k=1ri​​mijkt​mijkt​​,ri​1​,​​若k=1∑ri​​mijkt​>0若否​(10)
其中 m i j k t m^t_{ijk} mijkt​是补后数据中所有满足 X i = k , π ( X i ) = j X_i=k,\pi(X_i)=j Xi​=k,π(Xi​)=j的样本的权重之和。关于式(10)的具体推导过程，将在下一小节展开，此处暂先接受这个结论。
下面通过一个示例演示EM算法的迭代过程。
例 设有如下图(a)所示的一个贝叶斯网 N N N，其中所有变量均取二值1或2。下图(b)是从 N N N的联合分布中抽样得到的一组数据，其中有两个样本缺值。考虑用EM算法计算 N N N的参数的MLE估计。

首先随机选择一组初始参数 θ 0 \theta^0 θ0，设 θ 0 \theta^0 θ0如下：

• 第一轮迭代：
  执行E步：利用 θ 0 \theta^0 θ0修补数据：
  由于
  P ( X 2 = 1 ∣ D 3 , θ 0 ) = 4 / 5 P ( X 2 = 2 ∣ D 3 , θ 0 ) = 1 / 5 P(X_2=1|D_3,\bm\theta^0)=4/5\\ P(X_2=2|D_3,\bm\theta^0)=1/5 P(X2​=1∣D3​,θ0)=4/5P(X2​=2∣D3​,θ0)=1/5
  所以 D 3 D_3 D3​被如下两个碎权样本替代：
  D 3.1 = ( 1 , 1 , 1 ) [ 4 / 5 ] D 3.2 = ( 1 , 2 , 1 ) [ 1 / 5 ] D_{3.1}=(1,1,1)[4/5]\\ D_{3.2}=(1,2,1)[1/5] D3.1​=(1,1,1)[4/5]D3.2​=(1,2,1)[1/5]
  类似地， D 4 D_4 D4​也被两个碎权样本替换。修补后得到如下的碎权完整样本集：
  
  执行M步：根据式(10)，对该碎权完整样本集进行MLE估计，得到 θ 1 \theta^1 θ1如下：
  
• 第二轮迭代：
  执行E步：利用 θ 1 \theta^1 θ1修补数据：
  由于
  P ( X 2 = 1 ∣ D 3 , θ 1 ) = 81 / 82 P ( X 2 = 2 ∣ D 3 , θ 1 ) = 1 / 82 P ( X 2 = 1 ∣ D 4 , θ 1 ) = 1 / 82 P ( X 2 = 2 ∣ D 4 , θ 1 ) = 81 / 82 \begin{aligned} &P(X_2=1|D_3,\bm\theta^1)=81/82\\ &P(X_2=2|D_3,\bm\theta^1)=1/82\\ &P(X_2=1|D_4,\bm\theta^1)=1/82\\ &P(X_2=2|D_4,\bm\theta^1)=81/82 \end{aligned} ​P(X2​=1∣D3​,θ1)=81/82P(X2​=2∣D3​,θ1)=1/82P(X2​=1∣D4​,θ1)=1/82P(X2​=2∣D4​,θ1)=81/82​
  修补后得到如下的碎权完整样本集：
  
  执行M步：根据式(10)，对该碎权完整样本集进行MLE估计，得到 θ 2 \theta^2 θ2如下：
  

如此迭代下去，容易看到，这个过程快速向如下结果收敛：

3.2.3.3 EM算法推导

EM算法思想直观易懂，但关于『式(10)给出的是基于修补后碎权完整样本集的最大似然估计』这一结论并没有证明。此处将进行推导。
设 D = ( D 1 , ⋯   , D m ) D=(D_1,\cdots,D_m) D=(D1​,⋯,Dm​)是关于贝叶斯网 N = ( G , θ ) N=(G, \bm\theta) N=(G,θ)的一组i.i.d.样本集。对其中任一样本 D l D_l Dl​，设 X l X_l Xl​是 D l D_l Dl​中缺值变量的集合。设 θ t \bm\theta^t θt是关于参数 θ \bm\theta θ的当前估计， D t D^t Dt是基于 θ t \bm\theta^t θt将 D D D修补后得到的碎权完整样本集。定义 θ \bm\theta θ基于 D t D^t Dt的期望对数似然函数为：
l e ( θ ∣ D t ) = ∑ l = 1 m ∑ x l ∈ Ω X l P ( X l = x l ∣ D l t , θ t ) log ⁡ P ( D l t , X l = x l ∣ θ ) l_e(\bm\theta|D^t)=\sum_{l=1}^m\sum_{x_l\in\Omega_{X_l}}P(X_l=x_l|D^t_l,\bm\theta^t)\log P(D^t_l,X_l=x_l|\bm\theta) le​(θ∣Dt)=l=1∑m​xl​∈ΩXl​​∑​P(Xl​=xl​∣Dlt​,θt)logP(Dlt​,Xl​=xl​∣θ)
其中 P ( X l = x l ∣ D l t , θ t ) P(X_l=x_l|D^t_l,\bm\theta^t) P(Xl​=xl​∣Dlt​,θt)就是权重 ω x l \omega_{x_l} ωxl​​。当 X l = ∅ X_l=\empty Xl​=∅时，约定 P ( X l = x l ∣ D l t , θ t ) = 1 P(X_l=x_l|D^t_l,\bm\theta^t)=1 P(Xl​=xl​∣Dlt​,θt)=1，即完整样本的权重为1.
在EM算法第 t t t次迭代中：E步利用 θ t \bm\theta^t θt修补缺失数据，基于修补后的样本集可获得期望对数似然函数 l e ( θ ∣ D t ) l_e(\bm\theta|D^t) le​(θ∣Dt)；M步求使得 l e ( θ ∣ D t ) l_e(\bm\theta|D^t) le​(θ∣Dt)最大的 θ \bm\theta θ取值，即：
θ t + 1 = arg max ⁡ θ l ( θ ∣ D t ) \bm\theta^{t+1}=\argmax_{\bm\theta}l(\bm\theta|D^t) θt+1=θargmax​l(θ∣Dt)
下面讨论如何计算 θ t + 1 \bm\theta^{t+1} θt+1。首先，碎权样本 ( D l , X l = x l ) (D_l,X_l=x_l) (Dl​,Xl​=xl​)的特征函数为：
χ ( i , j , k : D l , X l = x l ) = { 1 , 若 在 ( D l , X l = x l ) 中 X i = k 且 π ( X i ) = j 0 , 若 否 \chi(i,j,k:D_l,X_l=x_l)=\left\{ \begin{aligned} 1,&\quad 若在(D_l,X_l=x_l)中X_i=k且\pi(X_i)=j\\ 0,&\quad 若否 \end{aligned} \right. χ(i,j,k:Dl​,Xl​=xl​)={1,0,​若在(Dl​,Xl​=xl​)中Xi​=k且π(Xi​)=j若否​
从而可得：
log ⁡ P ( D l , X l = x l ∣ θ ) = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i χ ( i , j , k : D l , X l = x l ) log ⁡ θ i j k \log P(D_l,X_l=x_l|\bm\theta)=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\chi(i,j,k:D_l,X_l=x_l)\log\theta_{ijk} logP(Dl​,Xl​=xl​∣θ)=i=1∑n​j=1∑qi​​k=1∑ri​​χ(i,j,k:Dl​,Xl​=xl​)logθijk​
于是，
l e ( θ ∣ D t ) = ∑ l = 1 m ∑ x l ∈ Ω X l P ( X l = x l ∣ D l t , θ t ) ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i χ ( i , j , k : D l , X l = x l ) log ⁡ θ i j k = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i ∑ l = 1 m ∑ x l ∈ Ω X l P ( X l = x l ∣ D l t , θ t ) χ ( i , j , k : D l , X l = x l ) log ⁡ θ i j k \begin{aligned} l_e(\bm\theta|D^t)&=\sum_{l=1}^m\sum_{x_l\in\Omega_{X_l}}P(X_l=x_l|D^t_l,\bm\theta^t)\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\chi(i,j,k:D_l,X_l=x_l)\log\theta_{ijk}\\ &=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}\sum_{l=1}^m\sum_{x_l\in\Omega_{X_l}}P(X_l=x_l|D^t_l,\bm\theta^t)\chi(i,j,k:D_l,X_l=x_l)\log\theta_{ijk} \end{aligned} le​(θ∣Dt)​=l=1∑m​xl​∈ΩXl​​∑​P(Xl​=xl​∣Dlt​,θt)i=1∑n​j=1∑qi​​k=1∑ri​​χ(i,j,k:Dl​,Xl​=xl​)logθijk​=i=1∑n​j=1∑qi​​k=1∑ri​​l=1∑m​xl​∈ΩXl​​∑​P(Xl​=xl​∣Dlt​,θt)χ(i,j,k:Dl​,Xl​=xl​)logθijk​​
定义
m i j k t = ∑ l = 1 m ∑ x l ∈ Ω X l P ( X l = x l ∣ D l t , θ t ) χ ( i , j , k : D l , X l = x l ) m^t_{ijk}=\sum_{l=1}^m\sum_{x_l\in\Omega_{X_l}}P(X_l=x_l|D^t_l,\bm\theta^t)\chi(i,j,k:D_l,X_l=x_l) mijkt​=l=1∑m​xl​∈ΩXl​​∑​P(Xl​=xl​∣Dlt​,θt)χ(i,j,k:Dl​,Xl​=xl​)
直观上， m i j k t m^t_{ijk} mijkt​是补后数据 D t D^t Dt中所有满足 X i = k , π ( X i ) = j X_i=k,\pi(X_i)=j Xi​=k,π(Xi​)=j的样本的权重之和。
从而有
l e ( θ ∣ D t ) = ∑ i = 1 n ∑ j = 1 q i ∑ k = 1 r i m i j k t log ⁡ θ i j k l_e(\bm\theta|D^t)=\sum_{i=1}^n\sum_{j=1}^{q_i}\sum_{k=1}^{r_i}m^t_{ijk}\log\theta_{ijk} le​(θ∣Dt)=i=1∑n​j=1∑qi​​k=1∑ri​​mijkt​logθijk​
由于有概率归一化条件 ∑ k = 1 r i θ i j k = 1 \sum_{k=1}^{r_i}\theta_{ijk}=1 ∑k=1ri​​θijk​=1，利用拉格朗日乘数法引入概率归一化条件，并求偏导为0的点即为最大值点，从而可求出：
θ i j k t + 1 = { m i j k t ∑ k = 1 r i m i j k t , 若 ∑ k = 1 r i m i j k t > 0 1 r i , 若 否 \theta_{ijk}^{t+1}=\left\{ \begin{aligned} &\frac{m^t_{ijk}}{\sum_{k=1}^{r_i}m^t_{ijk}},&\quad &若\sum_{k=1}^{r_i}m^t_{ijk}>0\\ &\frac{1}{r_i},&\quad &若否 \end{aligned} \right. θijkt+1​=⎩⎪⎪⎪⎨⎪⎪⎪⎧​​∑k=1ri​​mijkt​mijkt​​,ri​1​,​​若k=1∑ri​​mijkt​>0若否​
从而，式(10)推导完毕。

2.3.3.4 EM算法收敛性证明

EM算法在迭代过程中产生一个序列 { l ( θ 0 ∣ D ) , l ( θ 1 ∣ D ) , l ( θ 2 ∣ D ) , ⋯   } \{l(\bm\theta^0|D),l(\bm\theta^1|D),l(\bm\theta^2|D),\cdots\} {l(θ0∣D),l(θ1∣D),l(θ2∣D),⋯}，我们那来研究这个序列的性质。
首先给出一个定理。

定理1 设 θ t \bm\theta^t θt和 θ t + 1 \bm\theta^{t+1} θt+1分别是EM算法在第 t t t次和第 t + 1 t+1 t+1次迭代中所获得的结果，那么有：
l ( θ t + 1 ∣ D ) ≥ l ( θ t ∣ D ) l(\bm\theta^{t+1}|D)\geq l(\bm\theta^t|D) l(θt+1∣D)≥l(θt∣D)

证明： 由于
∑ X l P ( X l ∣ D l , θ ) = 1 P ( D l ∣ θ ) = P ( D l , X l ∣ θ ) P ( X l ∣ D l , θ ) \begin{aligned} &\sum_{X_l}P(X_l|D_l,\bm\theta)=1\\ &P(D_l|\bm\theta)=\frac{P(D_l,X_l|\bm\theta)}{P(X_l|D_l,\bm\theta)} \end{aligned} ​Xl​∑​P(Xl​∣Dl​,θ)=1P(Dl​∣θ)=P(Xl​∣Dl​,θ)P(Dl​,Xl​∣θ)​​
有
l ( θ ∣ D ) = ∑ l = 1 m log ⁡ P ( D l ∣ θ ) = ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ ) log ⁡ P ( D l ∣ θ ) = ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ ) log ⁡ P ( D l , X l ∣ θ ) P ( X l ∣ D l , θ ) = ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ ) log ⁡ P ( D l , X l ∣ θ ) − ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ ) log ⁡ P ( X l ∣ D l , θ ) = l e ( θ ∣ D ) − ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ ) log ⁡ P ( X l ∣ D l , θ ) \begin{aligned} l(\bm\theta|D)=&\sum_{l=1}^m\log P(D_l|\bm\theta)\\ =&\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta)\log P(D_l|\bm\theta)\\ =&\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta)\log\frac{P(D_l,X_l|\bm\theta)}{P(X_l|D_l,\bm\theta)}\\ =&\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta)\log P(D_l,X_l|\bm\theta)\\ &-\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta)\log P(X_l|D_l,\bm\theta)\\ =&l_e(\bm\theta|D)\\ &-\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta)\log P(X_l|D_l,\bm\theta) \end{aligned} l(θ∣D)=====​l=1∑m​logP(Dl​∣θ)l=1∑m​Xl​∑​P(Xl​∣Dl​,θ)logP(Dl​∣θ)l=1∑m​Xl​∑​P(Xl​∣Dl​,θ)logP(Xl​∣Dl​,θ)P(Dl​,Xl​∣θ)​l=1∑m​Xl​∑​P(Xl​∣Dl​,θ)logP(Dl​,Xl​∣θ)−l=1∑m​Xl​∑​P(Xl​∣Dl​,θ)logP(Xl​∣Dl​,θ)le​(θ∣D)−l=1∑m​Xl​∑​P(Xl​∣Dl​,θ)logP(Xl​∣Dl​,θ)​
由于在EM算法的M步中，
θ t + 1 = arg max ⁡ θ l e ( θ ∣ D t ) \bm\theta^{t+1}=\argmax_{\bm\theta}l_e(\bm\theta|D^t) θt+1=θargmax​le​(θ∣Dt)
因此有
l e ( θ t + 1 ∣ D t ) ≥   l e ( θ t ∣ D t ) l_e(\bm\theta^{t+1}|D^t)\geq\ l_e(\bm\theta^t|D^t) le​(θt+1∣Dt)≥ le​(θt∣Dt)
此外，根据信息不等式有：
∑ X l P ( X l ∣ D l , θ t ) log ⁡ P ( X l ∣ D l , θ t + 1 ) ≤ ∑ X l P ( X l ∣ D l , θ t ) log ⁡ P ( X l ∣ D l , θ t ) \sum_{X_l}P(X_l|D_l,\bm\theta^t)\log P(X_l|D_l,\bm\theta^{t+1})\leq\sum_{X_l}P(X_l|D_l,\bm\theta^t)\log P(X_l|D_l,\bm\theta^t) Xl​∑​P(Xl​∣Dl​,θt)logP(Xl​∣Dl​,θt+1)≤Xl​∑​P(Xl​∣Dl​,θt)logP(Xl​∣Dl​,θt)
从而有
l ( θ t + 1 ∣ D ) = l e ( θ t + 1 ∣ D t ) − − ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ t ) log ⁡ P ( X l ∣ D l , θ t + 1 ) ≥ l e ( θ t ∣ D t ) − ∑ l = 1 m ∑ X l P ( X l ∣ D l , θ t ) log ⁡ P ( X l ∣ D l , θ t ) = l ( θ t ∣ D ) \begin{aligned} l(\bm\theta^{t+1}|D)=&l_e(\bm\theta^{t+1}|D^t)-\\ &-\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta^t)\log P(X_l|D_l,\bm\theta^{t+1})\\ \geq&l_e(\bm\theta^t|D^t)\\ &-\sum_{l=1}^m\sum_{X_l}P(X_l|D_l,\bm\theta^t)\log P(X_l|D_l,\bm\theta^t)\\ =&l(\bm\theta^t|D) \end{aligned} l(θt+1∣D)=≥=​le​(θt+1∣Dt)−−l=1∑m​Xl​∑​P(Xl​∣Dl​,θt)logP(Xl​∣Dl​,θt+1)le​(θt∣Dt)−l=1∑m​Xl​∑​P(Xl​∣Dl​,θt)logP(Xl​∣Dl​,θt)l(θt∣D)​
定理得证。

推论1 设 { θ 0 , θ 1 , θ 2 , ⋯   } \{\bm\theta^0,\bm\theta^1,\bm\theta^2,\cdots\} {θ0,θ1,θ2,⋯}是EM算法所获得的参数估计序列，则序列 { l ( θ 0 ∣ D ) , l ( θ 1 ∣ D ) , l ( θ 2 ∣ D ) , ⋯   } \{l(\bm\theta^0|D),l(\bm\theta^1|D),l(\bm\theta^2|D),\cdots\} {l(θ0∣D),l(θ1∣D),l(θ2∣D),⋯}是单调递增的，并且收敛。

证明： 定理1已经证明了序列是单调递增的。
另一方面，
l ( θ ∣ D ) = ∑ l = 1 m log ⁡ P ( D l ∣ θ ) < 0 l(\bm\theta|D)=\sum_{l=1}^m\log P(D_l|\bm\theta)<0 l(θ∣D)=l=1∑m​logP(Dl​∣θ)<0
因此，0是该序列的一个上界，从而它必然收敛。推论得证。

EM算法虽然收敛，但不一定收敛于全局最优，很有可能收敛到局部最优。这与初始值的选择有很大关系。为了提高估计质量，可从不同的初始值出发，多次运行EM算法，并从所得结果中选择似然度最大的那个作为最终结果。

另外，EM算法收敛速度与数据缺失的多少有关。一般而言，数据缺失越多，收敛速度越慢。

3.2.4 小结

贝叶斯网参数估计有两种方法，极大似然估计和贝叶斯估计。极大似然估计适用于样本量足够大、可以用样本统计量估计整体统计量的场景下。而贝叶斯估计不仅适用于大样本量的场景，对于样本量不足但已有部分先验知识的场景，依然适用。本篇介绍了极大似然估计，下一篇将接着介绍贝叶斯估计。
在极大似然估计估计中，对于样本集存在缺值样本的情况，需要用EM算法进行迭代估计。对于EM算法的理解可分为几个层次：一是能套用算法解决基于缺值样本集的极大似然估计问题；二是能推导EM算法中的极大似然估计公式；三是能证明EM算法的收敛性；四是理解K-means、变分EM算法、Gibbs抽样、深度学习中的Wake-Sleep、对比散度算法等本质上都是EM算法的特例。本文仅达到第三个初级层次，对于EM算法的高级层次，这里给出一个进阶线索，可按图索翼去深入研究。