Chapter 11: Sampling Methods

一：Introduction

1：我们在做推理时，主要想要解决的是在后验概率$p(\vec z)$的情况下，计算$f(\vec z)$的期待值，即

$$E[f]=\int f(\vec z)p(\vec z)d\vec z$$ 当 $\vec z$是离散变量时，积分符号换成求和符号

2：为了解决上述问题，抽样方法(sampling method)的观点是从分布$p(\vec z)$中独立地获得一系列抽样值$\vec z^{l}$(l=1,…,L)，因此上述积分可以被近似为

$$\widehat{f}=\frac{1}{L}\sum_{l=1}^{L}f(\vec z^{l})$$ 不难获得 $E[\widehat{f}]=E[f]，var[\widehat{f}]=E[(f-E[f])^{2}]$

3：上述抽样方法的问题在于以下两点：第一点是提取的样品${\vec z^{l}}$或许不是相互独立地，因此我们实际上需要提取比有效样品数目更大的样品集合；第二点是如果在$p(\vec z)$很大的区域中$f(\vec z)$的值很小，反之亦然，则$f(\vec z)$的期待值被小概率的区域主导着，这就意味着为了实现足够大的准确度，我们需要比较大的样品数目$\vec z^{l}$

二：Basic Sampling Algorithms

1：Standard distributions

1）：假设我们有一个单变量z，其在(0,1)区间内均匀分布，同时我们想要从单变量分布p(y)中y的提取样品值，首先

$$p(y)=p(z)|\frac{dz}{dy}|\ \ \ \ \ (p(z)=1)$$
然后 $$z=h(y)=\int_{-\infty}^{y}p(\widehat{y})d\widehat{y}$$ 最后求其逆函数 $y=h^{-1}(z)$，有了这个关系式，当我们提取z的杨品值时，我们自然地也就提取了y的样品值

2）：对于多变量分布，我们则有以下关系式：

$$p(y_{1},...,y_{M})=p(z_{1},..,z_{M})|\frac{\partial(z_{1},...,z_{M})}{\partial(y_{1},...,y_{M})}|$$

3）：很明显该提取样品的方法需要目标概率分布的不确定积分很容易计算，以及对应的逆函数也很容易计算，因此该提取样品的方法只对有限数目的概率分布适用。

2：Rejection sampling

1）：我们想要从概率分布$p(z)$中提取样品，并且$p(z)=\frac{1}{Z_{p}}\widetilde{p}(z)$，$\widetilde{p}(z)$的值很容易被计算，但$Z_{p}$是未知的。Rejection sampling的步骤如下：

$\quad$(1)：需要有一个比较简单的易于抽样的proposal distribution q(z)， 同时我们需要有一个常数值k，使得$kq(z)>=\widetilde{p}(z)$，kq(z)被称之为comparsion function；

$\quad$(2)：首先从q(z)中提取一个样品$z_{0}$，然后再在区间$[0,kq(z_{0})]$均匀提取$u_{0}$；

$\quad$(3）：如果$u_{0}>\widetilde{p}(z_{0})$，则该样品值$z_{0}$被拒绝，否则被接收；

$\quad$(4）：可以证明，获得的样品值$z_{0}$满足概率分布p(z)。

2）：从q(z)中提取的样品被接受的概率：

$$p(accept)=\int \{\widetilde{p}(z)/kq(z)\}q(z)dz\\=\frac{1}{k}\int \widetilde{p}(z)dz$$
因此为了提高样品接受率，我们需要使得k在满足 $kq(z)>=\widetilde{p}(z)$的条件尽可能地小。

3）：Rejection sampling方法有如下几个缺点：第一个是在实际情况下一个好的proposal distribution和comparsion function是比较难找到的；第二个是随着输入变量维数的增加，样品接受率呈指数下降。因此rejection sampling方法虽然在一维或两维情况下是一个比较好的方法，但是不适合求解高维问题，然而该方法能够在求解高维问题中扮演一个子程序的角色。

3：Importance sampling

1）：Importance sampling解决的是怎么去近似$f(\vec z)$相对于概率分布$p(\vec z)$的期待值，但是该方法并不能够被应用于怎么从概率分布$p(\vec z)$中提取样品；

2）：类似于Rejection sampling，Importance sampling同样要求要有一个易于抽样的proposal distribution $q(\vec z)$，现在假设我们要从分布$p(\vec z)$中提取样品，并且$p(\vec z)=\widetilde{p}(\vec z)/Z_{p}$，$Z_{p}$的值未知，则

$$E[f]=\int f(\vec z)p(\vec z)d\vec z \\ \quad \ \approx \sum_{l=1}^{L}w_{l}f(\vec z^{(l)}) \\ w_{l}=\frac{\widetilde{p}(\vec z^{(l)})/q(\vec z^{(l)})}{\sum_{m}\widetilde{p}(\vec z^{(m)})/q(\vec z^{(m)})}$$

3)：Importance sampling方法成功的关键在于proposal distribution $q(\vec z)$ 与 $p(\vec z)f(\vec z)$的变化趋势匹配的好坏程度。该方法有如下两个缺点：第一点是$p(\vec z)f(\vec z)$经常只是在一个非常小的$\vec z$空间中值比较大，因此$f(\vec z)$的期待值或许只是被该一小片$\vec z$空间区域主导，也就是说有效样品的数目要远小于实际的样品数目，在实际中我们提取的大部分样品可能对于期待值的计算只起到非常小的的作用；第二点更糟的是如果在$p(\vec z)f(\vec z)$值很大的那一部分很小的$\vec z$空间区域内，$q(\vec z)$值很小，则没有或者很少的样品落入到那一部分区域，这样计算出来的期待值误差比较大。因此我们希望在$p(\vec z)f(\vec z)$值很大的区域内，抽样分布$q(\vec z)$的值不要太小。

4：Sampling-importance-resampling(SIR)

1)：类似于Rejection samplinging，SIR也要求要有一个易于抽样的分布$q(\vec z)$，SIR方法要求以下三个步骤：

$\quad$（1）：从$q(\vec z)$中提取L个样品${\vec z}^{(1)},...,{\vec z}^{(L)}$；

$\quad$（2）：使用下式计算权重$w_{1},...,w_{L}$

$$w_{l}=\frac{\widetilde{p}(\vec z^{(l)})/q(\vec z^{(l)})}{\sum_{m}\widetilde{p}(\vec z^{(m)})/q(\vec z^{(m)})}$$

$\quad$（3）：现在我们有一个离散分布(${\vec z}^{(1)},...,{\vec z}^{(L)}$)，每个成分对应的概率是($w_{1},...,w_{L}$)，我们就从该离散分布中提取L个样品作为从分布$p(\vec z)$中提取的样品

可以证明：在L趋于无穷大情况下，最终提取的L个样品服从的分布是$p(\vec z)$

2）：如果我们只是需要计算$f(\vec z)$相对于$p(\vec z)$的期待值，不需要从$p(\vec z)$中提取样品，则我们只需要计算下式

$$E[f(\vec z)] \approx \sum_{l=1}^{L}w_{l}f(\vec z_{l})$$

5：Sampling and the EM algorithm

1）：EM算法的E步骤

$$Q(\vec \theta,{\vec \theta}^{old})=\int p(Z|X,{\vec \theta}^{old})\ln{p(Z,X|\vec \theta)}dZ$$ 利用样品方法能够利用下式近似 $$Q(\vec \theta,{\vec \theta}^{old}) \approx \frac{1}{L}\sum_{l=1}^{L}\ln{p(Z^{l},X|\vec \theta)}$$
M步骤按照通常的方法最大化Q函数，这称之为Monte Carlo EM algorithm；

2）：在上面获得的Q函数中增加一项$\ln{p(\vec \theta)}$我们即可以找到使得$\vec \theta$后验分布最大值的$\vec \theta$。

三：Markov Chain Monte Carlo(MCMC)

1: Introdcution to MCMC

1)：像之前一样，假设要抽样的分布$p(\vec z)=\widetilde{p}(\vec z)/Z_{p}$，$Z_{p}$值未知，但$\widetilde{p}(\vec z)$的值很容易被计算。同时也和之前的抽样方法一样，MCMC也要求要有一个易于抽样的proposal distribution $q(\vec z)$。非常广泛的$p(\vec z)$适合于用MCMC方法来抽样

2)：MCMC方法大致思路如下：假设现在有了一个目前的抽样态${\vec{z}}^{(\tau)}$，则接下来从$q(\vec z|{\vec{z}}^{(\tau)})$中提取新的样品态${\vec z}^{*}$，如果该样品态满足一定的标准则接受该样品态，令${\vec z}^{(\tau+1)}={\vec z}^{*}$，如果不满足则抛弃该样品态${\vec z}^{*}$，令${\vec z}^{(\tau+1)}={\vec z}^{(\tau)}$，接下来再从$q(\vec z|{\vec{z}}^{(\tau+1)})$中提取新的样品态，如此重复下去，可以证明当$\tau->\infty$时，样品${\vec z}^{\tau}$服从的概率分布趋于$p(\vec z)$。从这里我们可以看出获得的一些列样品${\vec z}^{(1)},{\vec z}^{(2)},...$形成了一个Markov Chain。

3)：在Metropolis算法中，假定proposal distribution $q(\vec z)$是对称分布，即$q({\vec z}_{A}|{\vec z}_{B})=q({\vec z}_{B}|{\vec z}_{A})$，然后获得的样品${\vec z}^{*}$被接受的概率为

$$A({\vec z}^{*},{\vec z}^{(\tau)})=min(1,\frac{\widetilde{p}({\vec z}^{*})}{\widetilde{p}({\vec z}^{(\tau)})})$$ 因此我们可以在(0,1)区间均匀选择一个随机数u,如果 $u<A({\vec z}^{*},{\vec z}^{(\tau)})$则接受该样品 ${\vec z}^{*}$，否则拒绝该样品。

4)：我们应注意到在MCMC方法中，未来的样品态是条件在目前样品态中才得到的，目前的样品态是条件在过去的样品态中才得到的，所以我们获得的样品集合${\vec z}^{(1)},{\vec z}^{(2)},...$它们之间不是相互独立的，因为连续的样品态是高度相关联的。如果我们希望获得相互独立的样品，则我们可以每隔M个样品取一个样品态，如果M足够大，在实际用途中获得的样品集合可以认为是相互独立的。

2: Markov chains

1)：对于一系列随机变量${\vec z}^{(1)},...,{\vec z}^{(M)}$来说，其一阶Markov chain满足如下性质：

$$p({\vec z}^{(m+1)}|{\vec z}^{(1)},...,{\vec z}^{(m)})=p({\vec z}^{(m+1)}|{\vec z}^{(m)}) \ (m={1,...,M-1})$$
该公式意味着如果条件着系统目前的态，则系统未来的态和过去的态是独立不相关的；

2)：Markov chain结点变量的边际概率能够被表达成如下形式：

$$p({\vec z}^{(m+1)})=\sum_{{\vec z}^{(m)}}p({\vec z}^{(m+1)}|{\vec z}^{(m)})p({\vec z}^{(m)})$$ 如果在Markov chain每一步中都使得变量 ${\vec z}^{(m+1)}，{\vec z}^{(m)}$都遵循相同形式的概率分布，我们就说该概率分布相对于该Markov chain来说是invariant或者是stationary；

3)：使得概率分布$p(\vec z)$是Markov chain invariant分布的一个充分但不是必要的条件是下列等式的成立

$$p(\vec z)T(\vec z,{\vec z}^{'})=p({\vec z}^{'})T({\vec z}^{'},\vec z)$$ 其中 $T(\vec z,{\vec z}^{'})=p({\vec z}^{'}|\vec z)$称之为Markov chain的transition probability；

4)：当我们使用Markov chain从某一个概率分布$p^{*}(\vec z)$中抽取样品时，我们需要建立一个满足如下两个条件的Markov chain：第一个是使得$p^{*}(\vec z)$是该Markov chain的invariant distribution；第二个是当m->$\infty$时，$p({\vec z}^{(m)})$收敛于要求的invariant分布$p^{*}(\vec z)$。并且当这两个条件满足时，我们称此invariant分布为该Markov chain的平衡分布(equilibrium distribution)；

5)：实际中，我们经常利用一系列’base’ transitions $B_{1},...,B_{K}$通过如下的线性叠加方式来构造transition probability：

$$T({\vec z}^{'},\vec z)=\sum_{k=1}^{K}\alpha_{k}B_{k}({\vec z}^{'},\vec z)$$ 其中 $\alpha_{k}>=0$以及 $\sum_{k=1}^{K}\alpha_{k}=1$

四：Gibbs Sampling

1：Gibbs sampling是一个简单但被广泛使用的Markov chain Monte Carlo算法，其过程如下：

(1)：初始化$\{z_{i},i=1,...,M\}$

(2)：对于$\tau=1,...,T$：

• 从$p(z_{1}|{z_{2}}^{(\tau)},{z_{3}}^{(\tau)},...,{z_{M}}^{\tau})$中提取样品${z_{1}}^{(\tau+1)}$
• 从从$p(z_{2}|{z_{1}}^{(\tau+1)},{z_{3}}^{(\tau)},...,{z_{M}}^{\tau})$中提取样品${z_{2}}^{(\tau+1)}$
• …
• 从$p(z_{j}|{z_{1}}^{(\tau+1)},...,{z_{j-1}}^{(\tau)+1},{z_{j+1}}^{(\tau)}...,{z_{M}}^{\tau})$中提取样品${z_{j}}^{(\tau+1)}$
• …
• 从$p(z_{M}|{z_{1}}^{(\tau+1)},{z_{2}}^{(\tau+1)},...,{z_{M-1}}^{\tau+1})$中提取样品${z_{M}}^{(\tau+1)}$

2：Gibbs sampling的应用性依赖于从条件分布$p(z_{k}|{\vec z}_{\setminus k})$中提取样品的难易程度