贝叶斯网专题7：基于蒙特卡洛法的近似推理之重要性抽样

---

第一部分：贝叶斯网基础

1.1 信息论基础

1.2 贝叶斯网基本概念

1.3 变量独立性的图论分析

第二部分：贝叶斯网推理

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

2.2 团树传播算法

2.3 近似推理

前面两节介绍的都是贝叶斯网的精确推理方法，其核心是通过独立性和可复用性化简推理复杂度。然而在实际问题中，贝叶斯网的节点往往非常多，且连接稠密，即使通过变量消元、团树传播等等手段，依然会有很高的计算复杂度，从而精确推理不再适用。
其实，贝叶斯网解决的是与概率相关的推理问题，而概率是针对不确定性的统计特征，通过统计方法来近似解决不确定性推理，这才是贝叶斯网的价值所在。
贝叶斯网近似推理主要有两大类方法：一类是随机抽样方法，又叫蒙特卡洛方法，这是一类广泛应用于数值积分和统计物理的数值计算方法；另一类是变分法，这也是常用于统计物理的近似方法。
本讲先介绍基于蒙特卡洛法的贝叶斯网近似推理方法，下一讲将介绍基于变分的近似推理方法。

2.3.1 蒙特卡洛方法

蒙特卡洛法的基本思想是从某个概率分布随机抽样，生成一组符合概率分布的样本集，通过样本集来估计所要计算的统计量。蒙特卡洛方法又可分为两大类：重要性抽样法和马尔可夫蒙特卡洛法(MCMC)。它们的主要区别在于前者产生的样本之间相互独立，后者产生的样本相互关联。

2.3.1.1 重要性抽样法

1. 重要性抽样原理
   设$f(x)$是一组变量X在其定义域${\Omega }_X\subset R^n$上的可积函数。考虑计算积分：
   $$I={\int }_{{\Omega }_X}f(X)\mathrm{d}X$$
   通过重要性抽样来近似计算该积分：
   $$\begin{array}{rl}I& ={\int }_{{\Omega }_X}\frac{f(X)}{p(X)}p(X)\mathrm{d}X\\ & =E_{p(X)}\left[\frac{f(X)}{p(X)}\right]\\ & \approx \frac{1}{m}\sum _{i=1}^m\frac{f(D_i)}{p(D_i)},D_i\sim p(X)\end{array}\text{\hspace{1em}(1)}$$
   这里，X被看成一组随机变量，p(X)是X的一个联合概率分布，称为重要性分布，它满足以下条件：p的支撑包含f的支撑，即对X的任一取值x，若$f(X=x)\ne 0$，那么$p(X=x)\ne 0$。
   上式中约等号后面是抽样求和的过程，即从p(X)中独立抽取m个样本$D_1,\cdots ,D_m$，用$\frac{f(X)}{p(X)}$在该样本集上的平均值来近似其在分布p(X)上的期望。其符合大数定律，当样本容量足够大时，样本均值逼近总体期望。
   上式中的约等号的近似质量可以通过对约等号后面部分计算期望和方差来衡量。通过计算期望可以判断是否有偏，通过计算方差可以判断收敛性，以及推导收敛界与样本容量的关系。下面推导其无偏性，至于收敛界的推导与大数定律的几个界有关（Hoeffding界、Chernoff界），这里仅给出特殊情况下的定性判断，完整推导可参考Daphne所著《概率图模型》中第12章相关内容。
   首先分析期望：
   $$\begin{array}{rl}{E}_{p(X)}\left[\frac{1}{m}\sum _{i=1}^m\frac{f(D_i)}{p(D_i)}\right]& =\frac{1}{m}\sum _{i=1}^m{E}_{p(X)}\left[\frac{f(X)}{p(X)}\right]\\ & =\frac{1}{m}\sum _{i=1}^m{\int }_{{\Omega }_X}p(X)\frac{f(X)}{p(X)}dX\\ & =\frac{1}{m}m{\int }_{{\Omega }_X}f(X)dX\\ & =I\end{array}$$
   从而式(1)中的近似是无偏估计。
   接着分析方差：
   $$\begin{array}{rl}{D}_{p(X)}\left[\frac{1}{m}\sum _{i=1}^m\frac{f(D_i)}{p(D_i)}\right]& =\frac{1}{m^2}\sum _{i=1}^m{D}_{p(X)}\left[\frac{f(X)}{p(X)}\right]\\ & =\frac{1}{m}{D}_{p(X)}\left[\frac{f(X)}{p(X)}\right]\end{array}\text{\hspace{1em}(2)}$$
   由此可见，该估计的收敛性与采样样本容量m以及重要性分布p(X)与被积函数f(X)的相似度有关。当m足够大时，方差趋近于0，估计值收敛于真实值；当p(X)的变化趋势趋近于f(X)时，方差也趋近于0。我们不妨观察该特殊情况：$f(X)>0$，$p(X)=\frac{f(X)}{I}$。在该条件下，对式(2)可进一步计算：
   $$\begin{array}{rl}\frac{1}{m}{D}_{p(X)}\left[\frac{f(X)}{p(X)}\right]& =\frac{1}{m}{D}_{p(X)}[I]\\ & =0\end{array}$$
   此时的p(X)称为最优重要性分布。当然这是理论分析的情况，实际问题本就是求取积分I，构造$p(X)=\frac{f(X)}{I}$犹如坐在板凳上想将板凳搬起来一样不可行。但这给了我们一个方向性的指导，我们选择的重要性分布应该尽量与被积函数的趋势相近，从而可以用较少的抽样次数来达到足够的近似质量，加快抽样的收敛过程。
   以上讨论是针对连续变量的，它们可自然地推广到离散变量的情况，只需将积分运算相应地替换为求和运算即可。

2. 利用重要性抽样进行概率推理
   考虑一个贝叶斯网N，用X表示其中所有变量的集合，P(X)表示N的联合概率分布。设观测到证据变量E=e。下面将讨论如何近似计算一组查询变量Q取某值q的后验概率P(Q=q|E=e)。
   按照条件概率的定义，有：
   $$P(Q=q|E=e)=\frac{P(Q=q,E=e)}{P(E=e)}$$
   我们定义一个指示函数，其中X为该函数的自变量集合，Y为其子集，$Y\subset X$。当给定一个自变量的取值X=x，若其中的自变量子集Y=y，则该指示函数值为1，否则为0：
   $$I_{Y=y}(X)=\{\begin{array}{r}1,ifY=y\\ 0,ifY\ne y\end{array}$$
   利用该指示函数，P(Q=q,E=e)和P(E=e)可以表示成如下期望的形式，并利用抽样样本均值近似期望：
   $$\begin{array}{rl}P(Q=q,E=e)& =\sum _X{I}_{Q=q}(X)I_{E=e}(X)p(X)\\ & =E_{p(X)}\left[I_{Q=q}(X)I_{E=e}(X)\right]\\ & \approx \frac{1}{m}\sum _{i=1}^m{I}_{Q=q}(D_i)I_{E=e}(D_i),D_i\sim p(X)\end{array}\text{\hspace{1em}(3)}$$
   $$\begin{array}{rl}P(E=e)& =\sum _X{I}_{E=e}(X)p(X)\\ & =E_{p(X)}\left[I_{E=e}(X)\right]\\ & \approx \frac{1}{m}\sum _{i=1}^m{I}_{E=e}(D_i),D_i\sim p(X)\end{array}\text{\hspace{1em}(4)}$$
   从而可得后验概率的抽样估计：
   $$P(Q=q|E=e)\approx \frac{{\sum }_{i=1}^m{I}_{Q=q}(D_i)I_{E=e}(D_i)}{{\sum }_{i=1}^m{I}_{E=e}(D_i)},D_i\sim p(X)\text{\hspace{1em}(5)}$$
   前面已经分析过，对P(Q=q,E=e)和P(E=e)的抽样估计是无偏的，但对P(Q=q|E=e)的抽样估计是有偏的。由于P(Q=q|E=e)是两个累加相除的形式，计算其在p(X)上的期望比较困难，我们不妨通过一个特殊情况来观察它的有偏性，即只采样一个样本的情况，此时m=1：
   $$\begin{array}{rl}{E}_{p(X)}\left[\frac{{\sum }_{i=1}^m{I}_{Q=q}(D_i)I_{E=e}(D_i)}{{\sum }_{i=1}^m{I}_{E=e}(D_i)}\right]& =E_{p(X)}\left[\frac{I_{Q=q}(D_1)I_{E=e}(D_1)}{I_{E=e}(D_1)}\right]\\ & =E_{p(X)}\left[I_{Q=q}(D_1)\right]\\ & =\sum _X{I}_{Q=q}(X)p(X)\end{array}$$
   另一方面：
   $$P(Q=q|E=e)=\sum _X\frac{I_{Q=q}(X)I_{E=e}(X)p(X)}{p(E=e)}\text{\hspace{1em}(6)}$$
   两者并不相等，因此式(5)对P(Q=q|E=e)的抽样估计是有偏的。所幸的是，随着样本容量的增加，该抽样估计逐渐收敛于真实值，因为当m趋于无穷时，分子和分母都趋近于各自的真实值。

3. 接受-拒绝抽样
   在上式(3)和(4)中，我们很自然地选择了p(X)为重要性分布进行抽样估计，这种抽样方法称为『接受-拒绝抽样』。当然，我们也可以选择其它的重要性分布进行抽样估计，甚至会产生更好的效果，这将在下一小节的似然加权抽样中介绍。这里先介绍『接受-拒绝抽样』这种最直观的抽样方式。
   贝叶斯网是一个有向无圈图，其中的变量存在一个拓扑序，其中每一个节点都在它的子节点之前出现。
   接受-拒绝抽样需要从p(X)中抽取样本。注意到p(X)可以分解为：
   $$p(X)=\prod _{x\in X}p(X|\pi (X))$$
   因此可以按照贝叶斯网N的拓扑序对其中的变量逐个进行抽样：若抽样变量Y为根节点，则按分布p(Y)进行抽样；若抽样变量Y不是根节点，则按分布$p(Y|\pi (Y)=r)$进行抽样，其中$\pi (Y)=r$是Y的父节点的抽样结果，在对Y抽样时是已知的。该抽样过程称为顺序抽样。
   所谓抽样，可通过一个随机数发生器来产生一个[0,1]之间均匀分布的随机数，若某变量Y的概率分布如下：
   $$\begin{gathered} p(Y=y_1)=0.2 \\ p(Y=y_2)=0.4 \\ p(Y=y_3)=0.4 \end{gathered}$$
   若所产生的随机数在0~0.2之间，则对Y的本次采样结果为$y_1$；若随机数在0.2~0.6之间，则本次采样结果为$y_2$，若随机数在0.6~1之间，则采样结果为$y_3$。
   下面我们用一个具体案例来演示该抽样推理的过程。
   例1 对下图所示的贝叶斯网顺序抽样。图中标出了按照<CRSW>拓扑序的一次抽样过程。
   
   假设这样的抽样共进行了m次，获得了m个独立的样本$D_1,\cdots ,D_m$，其中满足E=e的有$m_e$个，在这$m_e$个样本中，进一步满足Q=q的有$m_{q,e}$个，根据式(5)有：
   $$p(Q=q|E=e)\approx \frac{m_{q,e}}{m_e}$$
   假设所查询的变量Q为R=T，证据变量E为S=T，总样本数m=100，其中有75个满足S=F，与证据不一致，被舍弃。在剩余的25个样本中，有18个满足R=T，因此得到的后验概率估计为：$p(R=T|S=T)\approx 18/25=0.72$，接近其精确值为0.7。
   在接受-拒绝抽样中，所产生的与证据E=e不一致的那些样本相当于被舍弃，这也是『接受-拒绝』名字的由来。
   该抽样方法简单易行，缺点是当p(E=e)很小时，大量样本会被丢弃，收敛速度慢。实际上，从式(6)可以看出，后验概率问题的最优重要性分布是$\frac{I_{E=e}(X)p(X)}{p(E=e)}$，而抽样使用的重要性分布为p(X)，两者的差别显著，前者概率集中在E=e的取值处，而后者在在这些区域的概率值却很小。随着p(E=e)的减小，所抽得的与证据一致的样本数也会减少，因此大量样本被舍弃，造成计算资源浪费。
   为了解决该问题，可以使用下面介绍的似然加权法。

4. 似然加权抽样
   为了避免接受-拒绝抽样中因舍弃样本而影响算法效率，提出了似然加权抽样，其是重要性抽样的一个特例。其与接受-拒绝抽样的区别在于，不再把p(X)当做重要性分布进行抽样，而是精心构造了另一个重要性分布，具体的构造过程如下。
   首先，我们定义一个证据变量E上的概率分布：
   $$p^{\prime }(E)=I_{E=e}(E)$$
   当E=e时，$p^{\prime }(E=e)=1$；当$E\ne e$时，$p^{\prime }(E\ne e)=0$。从而$p^{\prime }(E)$满足概率归一性：${\sum }_E{p}^{\prime }(E)=1$.
   贝叶斯网上的所有变量为X，Z是除去证据变量E的其它变量，Z=X\E。我们再定义一个条件概率分布：
   $$p^{\prime }(Z|E)=\prod _{Y\in Z}p(Y|\pi (Y))$$
   于是有：
   $$\begin{array}{rl}{p}^{\prime }(X)& =p^{\prime }(Z,E)\\ & =p^{\prime }(E)p^{\prime }(Z|E)\\ & =I_{E=e}(E)\prod _{Y\in Z}p(Y|\pi (Y))\end{array}$$
   由于$I_{E=e}(E)=I_{E=e}(Z,E)=I_{E=e}(X)$，上式可进一步变为：
   $$p^{\prime }(X)=I_{E=e}(X)\prod _{Y\in Z}p(Y|\pi (Y))$$
   以$p^{\prime }(X)$为重要性分布，我们可以从$p^{\prime }(X)$抽样m个样本$D_1,\cdots ,D_m$。从上式可知，抽样方法为：按照拓扑序对贝叶斯网中每个变量Y进行抽样，当$Y\in Z$，不是证据变量时，抽样方式与接受-拒绝抽样的抽样方法一致；当$Y\in E$，为证据变量时，则直接给Y赋E中的值。按照该方法抽样生成的样本符合$p^{\prime }(X)$分布，且能保证抽样生成的每一个样本都与证据一致，无需舍弃。
   根据该重要性分布，我们重新对式(3)和式(4)进行抽样估计：
   $$\begin{array}{rl}P(Q=q,E=e)& =\sum _X{I}_{Q=q}(X)I_{E=e}(X)p(X)\\ & =\sum _X\frac{I_{Q=q}(X)I_{E=e}(X)p(X)}{p^{\prime }(X)}{p}^{\prime }(X)\\ & =E_{p^{\prime }(X)}\left[\frac{I_{Q=q}(X)I_{E=e}(X)p(X)}{p^{\prime }(X)}\right]\\ & \approx \frac{1}{m}\sum _{i=1}^m\frac{I_{Q=q}(D_i)I_{E=e}(D_i)p(D_i)}{p^{\prime }(D_i)}\end{array}\text{\hspace{1em}(7)}$$
   $$\begin{array}{rl}P(E=e)& =\sum _X{I}_{E=e}(X)p(X)\\ & =\sum _X\frac{I_{E=e}(X)p(X)}{p^{\prime }(X)}{p}^{\prime }(X)\\ & =E_{p(X^{\prime })}\left[\frac{I_{E=e}(X)p(X)}{p^{\prime }(X)}\right]\\ & \approx \frac{1}{m}\sum _{i=1}^m\frac{I_{E=e}(D_i)p(D_i)}{p^{\prime }(D_i)}\end{array}\text{\hspace{1em}(8)}$$
   其中：
   $$\begin{array}{rl}{I}_{E=e}(D_i)& =1\\ p(D_i)& =\prod _{Y\in X}p(Y|\pi (Y)){|}_{X=D_i}\\ p^{\prime }(D_i)& =I_{E=e}(D_i)\prod _{Y\in Z}p(Y|\pi (Y)){|}_{X=D_i}\\ & =\prod _{Y\in Z}p(Y|\pi (Y)){|}_{X=D_i}\end{array}$$
   代入式(7)和式(8)可得：
   $$\begin{array}{rl}P(Q=q,E=e)& \approx \frac{1}{m}\sum _{i=1}^m\frac{I_{Q=q}(D_i)I_{E=e}(D_i)p(D_i)}{p^{\prime }(D_i)}\\ & =\frac{1}{m}\sum _{i=1}^m\left(I_{Q=q}(D_i)\frac{{\prod }_{Y\in X}p(Y|\pi (Y)){|}_{X=D_i}}{{\prod }_{Y\in Z}p(Y|\pi (Y)){|}_{X=D_i}}\right)\\ & =\frac{1}{m}\sum _{i=1}^m\left(I_{Q=q}(D_i)\prod _{Y\in E}p(Y|\pi (Y)){|}_{X=D_i}\right)\\ & =\frac{1}{m}\sum _{i=1}^m{I}_{Q=q}(D_i)\omega (D_i)\end{array}\text{\hspace{1em}(9)}$$
   $$\begin{array}{rl}P(E=e)& \approx \frac{1}{m}\sum _{i=1}^m\frac{I_{E=e}(D_i)p(D_i)}{p^{\prime }(D_i)}\\ & =\frac{1}{m}\sum _{i=1}^m\frac{{\prod }_{Y\in X}p(Y|\pi (Y)){|}_{X=D_i}}{{\prod }_{Y\in Z}p(Y|\pi (Y)){|}_{X=D_i}}\\ & =\frac{1}{m}\sum _{i=1}^m\prod _{Y\in E}p(Y|\pi (Y)){|}_{X=D_i}\\ & =\frac{1}{m}\sum _{i=1}^m\omega (D_i)\end{array}\text{\hspace{1em}(10)}$$
   其中$\omega (D_i)={\prod }_{Y\in E}p(Y|\pi (Y)){|}_{X=D_i}$，是观察到E=e的似然率。所谓似然率，即给定概率分布下对观察到的现象的预测。比如抛硬币，给定正面向上的概率为0.5，如果抛5次硬币，预测能观察到3次朝上这一现象的概率则称为似然率。
   从而可的后验概率分布的估计：
   $$p(Q=q|E=e)\approx \frac{{\sum }_{i=1}^m{I}_{Q=q}(D_i)\omega (D_i)}{{\sum }_{i=1}^m\omega (D_i)},D_i\sim p^{\prime }(X)\text{\hspace{1em}(11)}$$
   与式(5)接受-拒绝抽样相比，似然加权算法相当于为每个样本$D_i$赋予了一个似然率权重$\omega (D_i)$，这也正是似然加权算法名称的由来。
   例2 仍以例1中的贝叶斯网为例，如下图为一次似然加权抽样过程。其中证据变量E={R=T,W=T}，查询变量Q={S=T}。若进行了m次似然加权抽样，得到m个样本$D_i$和对应的似然率权重$\omega (D_i)$以及加权指示$I_{S=T}(D_i)\omega (D_i)$，根据式(11)则可计算出后验概率P(S=T|R=T,W=T)的估计值。
   

本讲介绍了蒙特卡洛抽样中的重要性抽样方法，包括接受-拒绝抽样和似然加权抽样；下一讲将接着介绍蒙特卡洛抽样中的另一种方法：马尔可夫蒙特卡洛抽样。