PRML学习总结（11）——Sampling Methods

对于⼤多数实际应⽤中的概率模型来说，精确推断是不可⾏的，因此我们不得不借助与某种形式的近似。在第10章中，我们讨论了基于确定性近似的推断⽅法，它包括诸如变分贝叶斯⽅法以及期望传播。这⾥，我们考虑基于数值采样的近似推断⽅法，也被称为蒙特卡罗（Monte Carlo）⽅法。
在预测时，我们往往需要计算$$\mathbb{E}[f]=\int f(\mathbf{z})p(\mathbf{z})\mathrm{d}\mathbf{z}$$这个积分往往太复杂，无法计算。

计算上面的期望，可以从分布$p(\mathbf{z})$中独立抽样${\mathbf{z}}^{(l)}$，$l=1,\dots ,L$，则$$\hat{f}=\frac{1}{L}\sum _{l=1}^{L}f\left({\mathbf{z}}^{(l)}\right)$$使得$\mathbb{E}[\hat{f}]=\mathbb{E}[f]$，这样估计的方差为$$\mathrm{var}[\hat{f}]=\frac{1}{L}\mathbb{E}\left[(f-\mathbb{E}[f]{)}^2\right]$$值得强调的⼀点是，估计的精度不依赖于$\mathbf{z}$的维度，并且原则上，对于数量相对较少的样本${\mathbf{z}}^{(l)}$，可能会达到较⾼的精度。在实际应⽤中，10个或者20个独⽴的样本就能够以⾜够⾼的精度对期望做出估计。
有时，${\mathbf{z}}^{(l)}$不是独立的，从而导致有效样本⼤⼩可能远远⼩于表⾯上的样本⼤⼩。在上图中，我们还发现如果$f(\mathbf{z})$在$p(\mathbf{z})$很大的地方值很小，反之亦然，那么期望可能由⼩概率的区域控制，表明为了达到⾜够的精度，需要相对较⼤的样本⼤⼩。

11.1 Basic Sampling Algorithms

11.1.1 Standard distributions

首先考虑如何从均匀分布产生一些简单的分布，假设$z$为（0，1）均匀分布，通过$y=f(z)$得到$y$的分布$$p(y)=p(z)\left|\frac{dz}{dy}\right|$$由于$p(z)=1$，则$$z=h(y)\equiv {\int }_{-\infty }^{y}p(\hat{y})\mathrm{d}\hat{y}$$则$y=h^{-1}(z)$，如下图所示

下面举两个具体的例子，指数分布$$p(y)=\lambda \exp (-\lambda y)$$和柯西分布$$p(y)=\frac{1}{\pi }\frac{1}{1+y^2}$$对于多元变量，则有$$p\left(y_1,\dots ,y_M\right)=p\left(z_1,\dots ,z_M\right)\left|\frac{\partial \left(z_1,\dots ,z_M\right)}{\partial \left(y_1,\dots ,y_M\right)}\right|$$作为变换⽅法的最后⼀个例⼦，我们考虑Box-Muller⽅法，⽤于⽣成⾼斯概率分布的样本。⾸先，假设我们⽣成⼀对均匀分布的随机变量$z_1,z_2\in (-1,1)$，我们可以这样⽣成：对$(0,1)$上的均匀分布的变量使⽤$z\to 2z-1$的⽅式进⾏变换。接下来，我们丢弃那些不满⾜$z_1^2+{\dot{z}}_2^2⩽1$的点对。这产⽣出单位圆内部的⼀个均匀分布，且$p\left(z_1,z_2\right)=1/\pi$ ，然后，对于每对$z_1,z_2$，我们计算$$\begin{array}{rl}{y}_1& =z_1{\left(\frac{-2\ln z_1}{r^2}\right)}^{1/2}\\ y_2& =z_2{\left(\frac{-2\ln z_2}{r^2}\right)}^{1/2}\end{array}$$其中$r^2=z_1^2+z_2^2$则$$\begin{array}{rl}p\left(y_1,y_2\right)& =p\left(z_1,z_2\right)\left|\frac{\partial \left(z_1,z_2\right)}{\partial \left(y_1,y_2\right)}\right|\\ & =\left[\frac{1}{\sqrt{2\pi }}\exp \left(-y_1^2/2\right)\right]\left[\frac{1}{\sqrt{2\pi }}\exp \left(-y_2^2/2\right)\right]\end{array}$$得到标准高斯后，如重参数那样直接得到一般的高斯分布，对于多元高斯来说，先进行分解$$\mathbf{\Sigma }=\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$$则$\mathbf{y}=\mathbf{\mu }+\mathbf{L}\mathbf{z}$就会得到均值为$\mathbf{\mu }$，方差为$\mathbf{\Sigma }$的多元高斯，其中$\mathbf{z}$为标准多元高斯分布。
显然，变换⽅法依赖于它能够进⾏计算所需的概率分布，并且能够求所需的概率分布的不定积分的反函数。这样的计算只对于⼀些⾮常有限的简单的概率分布可⾏，因此我们必须寻找⼀些更加⼀般的⽅法。

11.1.2 Rejection sampling

拒绝采样框架使得我们能够在满⾜某些限制条件的情况下，从相对复杂的概率分布中采样。⾸先，我们考虑单变量分布，然后接下来讨论对于多维情形的推⼴。
考虑分布$p(\mathbf{z})$，直接对它进行采样无从下手。进一步考虑$$p(z)=\frac{1}{Z_p}\tilde{p}(z)$$我们可以很轻易计算$\tilde{p}(z)$，但$Z_p$不知道。为了使用拒绝采样，我们引入一个简单的可采样的分布$q(z)$，被叫做proposal distribution。然后使得常数$k$，有$kq(z)⩾\tilde{p}(z)$

计算可得接受概率为 p (  accept  ) = ∫ { p ~ ( z ) / k q ( z ) } q ( z ) d z = 1 k ∫ p ~ ( z ) d z \begin{aligned} p(\text { accept }) &=\int\{\widetilde{p}(z) / k q(z)\} q(z) \mathrm{d} z \\ &=\frac{1}{k} \int \widetilde{p}(z) \mathrm{d} z \end{aligned} p( accept )​=∫{p ​(z)/kq(z)}q(z)dz=k1​∫p ​(z)dz​为了让接受率大，就必须保证$k$尽可能小，但又必须满足上面的约束。

11.1.3 Adaptive rejection sampling

在许多我们希望应⽤拒绝采样的情形中，确定概率分布$q(z)$的⼀个合适的解析形式是很困难的。另⼀种确定其函数形式的⽅法是基于概率分布$p(z)$的值直接构建函数形式.对于$p(z)$是对数凹函数的情形，即$lnp(z)$的导数是$z$的单调⾮增函数时，界限函数的构建是相当简单的。下图给出了⼀个合适的界限函数的构建的例⼦。

$\ln p(z)$和它的梯度在预先设定的一些点上计算得到，从而可以确定该直线的方程式子，则取exp后得到分段分布刚好为指数分布 q ( z ) = k i λ i exp ⁡ { − λ i ( z − z i − 1 ) } z i − 1 &lt; z ⩽ z i q(z)=k_{i} \lambda_{i} \exp \left\{-\lambda_{i}\left(z-z_{i-1}\right)\right\} \quad z_{i-1}&lt;z \leqslant z_{i} q(z)=ki​λi​exp{−λi​(z−zi−1​)}zi−1​<z⩽zi​采样可以采取之前的inverse CDF方法即可。

11.1.4 Importance sampling

之所以进行采样，是为了计算期望，重要采样（importance sampling）的⽅法提供了直接近似期望的框架，但是它本⾝并没有提供从概率分布$p(\mathbf{z})$中采样的⽅法。
一般来说$p(\mathbf{z})$的PDF是知道的，就是很难进行采样，一个很简单的策略是将$\mathbf{z}$空间网格化，则$$\mathbb{E}[f]\simeq \sum _{l=1}^{L}p\left({\mathbf{z}}^{(l)}\right)f\left({\mathbf{z}}^{(l)}\right)$$但是这个方法对于维度高的数据来说就很费时，而且很大部分分布其实是仅仅集中在$\mathbf{z}$空间的一部分，那么所有空间均匀网格化十分没必要，因此如果能采样到这些集中的区域，或者更理想的是$p(\mathbf{z})f(\mathbf{z})$大的区域，那么效率将会提升！

因此可以考虑一个易于采样的proposal 分布$q(\mathbf{z})$，则$$\begin{array}{rl}\mathbb{E}[f]& =\int f(\mathbf{z})p(\mathbf{z})\mathrm{d}\mathbf{z}\\ & =\int f(\mathbf{z})\frac{p(\mathbf{z})}{q(\mathbf{z})}q(\mathbf{z})\mathrm{d}\mathbf{z}\\ & \simeq \frac{1}{L}\sum _{l=1}^L\frac{p\left({\mathbf{z}}^{(l)}\right)}{q\left({\mathbf{z}}^{(l)}\right)}f\left({\mathbf{z}}^{(l)}\right)\end{array}$$其中$r_l=p\left({\mathbf{z}}^{(l)}\right)/q\left({\mathbf{z}}^{(l)}\right)$为importance weights。修正了由于从错误的概率分布中采样引⼊的偏差。注意，与拒绝采样不同，所有⽣成的样本都被保留。
同样地对于非归一化的分布$p(\mathbf{z})=\tilde{p}(\mathbf{z})/Z_p$，其中$Z_p$未知$$\begin{array}{rl}\mathbb{E}[f]& =\int f(\mathbf{z})p(\mathbf{z})\mathrm{d}\mathbf{z}\\ & =\frac{Z_q}{Z_p}\int f(\mathbf{z})\frac{\tilde{p}(\mathbf{z})}{\tilde{q}(\mathbf{z})}q(\mathbf{z})\mathrm{d}\mathbf{z}\\ & \simeq \frac{Z_q}{Z_p}\frac{1}{L}\sum _{l=1}^L{\tilde{r}}_{l}f\left({\mathbf{z}}^{(l)}\right)\end{array}$$其中${\tilde{r}}_l=\tilde{p}\left({\mathbf{z}}^{(l)}\right)/\tilde{q}\left({\mathbf{z}}^{(l)}\right)$，用同样的采样点，可以得到$$\begin{array}{rl}\frac{Z_p}{Z_q}& =\frac{1}{Z_q}\int \tilde{p}(\mathbf{z})\mathrm{d}\mathbf{z}=\int \frac{\tilde{p}(\mathbf{z})}{\tilde{q}(\mathbf{z})}q(\mathbf{z})\mathrm{d}\mathbf{z}\\ & \simeq \frac{1}{L}\sum _{l=1}^L{\tilde{r}}_l\end{array}$$因此$$\mathbb{E}[f]\simeq \sum _{l=1}^L{w}_{l}f\left({\mathbf{z}}^{(l)}\right)$$其中$$w_l=\frac{{\tilde{r}}_l}{{\sum }_m{\tilde{r}}_m}=\frac{\tilde{p}\left({\mathbf{z}}^{(l)}\right)/q\left({\mathbf{z}}^{(l)}\right)}{{\sum }_m\tilde{p}\left({\mathbf{z}}^{(m)}\right)/q\left({\mathbf{z}}^{(m)}\right)}$$与拒绝采样一样，其成功严重依赖proposal分布选择的合适程度。因此，重要性采样⽅法的⼀个主要的缺点是它具有产⽣任意错误的结果的可能性，并且这种错误⽆法检测。这也强调了采样分布$q(\mathbf{z})$的⼀个关键的要求，即它不应该在$p(\mathbf{z})$可能较⼤的区域中取得较⼩的值或者为零的值。

11.1.5 Sampling-importance-resampling

拒绝采样⽅法部分依赖于它能够成功确定常数$k$的⼀个合适的值。对于许多对概率分布$p(\mathbf{z})$和$q(\mathbf{z})$来说，确定⼀个合适的$k$值是不现实的，因为任意的⾜够⼤的$k$值都能够保证产⽣所求的分布的上界，但是这会产⽣相当⼩的接受率。
sampling-importance-resampling (SIR)跟RS一样，利用$q(\mathbf{z})$，但是不需要求$k$值。首先从$q(\mathbf{z})$抽样${\mathbf{z}}^{(1)},\dots ,{\mathbf{z}}^{(L)}$，利用$$w_l=\frac{{\tilde{r}}_l}{{\sum }_m{\tilde{r}}_m}=\frac{\tilde{p}\left({\mathbf{z}}^{(l)}\right)/q\left({\mathbf{z}}^{(l)}\right)}{{\sum }_m\tilde{p}\left({\mathbf{z}}^{(m)}\right)/q\left({\mathbf{z}}^{(m)}\right)}$$得到$\left(w_1,\dots ,w_L\right)$。最后，$L$个样本的第⼆个集合从离散概率分布(${\mathbf{z}}^{(1)},\dots ,{\mathbf{z}}^{(L)}$)中抽取，概率由权值$\left(w_1,\dots ,w_L\right)$给定。⽣成的$L$个样本只是近似地服从$p(\mathbf{z})$，但是在极限$L\to \infty$的情况下，分布变为了正确的分布。为了说明这⼀点，考虑⼀元变量的情形，并且注意重新采样的值的累积分布为$$\begin{array}{rl}p(z⩽a)& =\sum _{l:z^{(l)}⩽a}{w}_l\\ & =\frac{{\sum }_{l}I\left(z^{(l)}⩽a\right)\tilde{p}\left(z^{(l)}\right)/q\left(z^{(l)}\right)}{{\sum }_l\tilde{p}\left(z^{(l)}\right)/q\left(z^{(l)}\right)}\end{array}$$当$L\to \infty$，并且假设概率分布进⾏了适当的正则化，我们可以将求和替换为积分，权值为原始的采样分布$q(z)$，即 p ( z ⩽ a ) = ∫ I ( z ⩽ a ) { p ~ ( z ) / q ( z ) } q ( z ) d z ∫ { p ~ ( z ) / q ( z ) } q ( z ) d z = ∫ I ( z ⩽ a ) p ~ ( z ) d z ∫ p ~ ( z ) d z = ∫ I ( z ⩽ a ) p ( z ) d z \begin{aligned} p(z \leqslant a) &amp;=\frac{\int I(z \leqslant a)\{\widetilde{p}(z) / q(z)\} q(z) \mathrm{d} z}{\int\{\widetilde{p}(z) / q(z)\} q(z) \mathrm{d} z} \\ &amp;=\frac{\int I(z \leqslant a) \widetilde{p}(z) \mathrm{d} z}{\int \widetilde{p}(z) \mathrm{d} z} \\ &amp;=\int I(z \leqslant a) p(z) \mathrm{d} z \end{aligned} p(z⩽a)​=∫{p ​(z)/q(z)}q(z)dz∫I(z⩽a){p ​(z)/q(z)}q(z)dz​=∫p ​(z)dz∫I(z⩽a)p ​(z)dz​=∫I(z⩽a)p(z)dz​这个跟拒绝采样一样，proposal分布越接近$p(z)$，则越精确。利用此求各阶矩时，$$\begin{array}{rl}\mathbb{E}[f(\mathbf{z})]& =\int f(\mathbf{z})p(\mathbf{z})\mathrm{d}\mathbf{z}\\ & =\frac{\int f(\mathbf{z})[\tilde{p}(\mathbf{z})/q(\mathbf{z})]q(\mathbf{z})\mathrm{d}\mathbf{z}}{\int [\tilde{p}(\mathbf{z})/q(\mathbf{z})]q(\mathbf{z})\mathrm{d}\mathbf{z}}\\ & \simeq \sum _{l=1}^L{w}_{l}f\left({\mathbf{z}}_l\right)\end{array}$$

11.1.6 Sampling and the EM algorithm

蒙特卡罗⽅法除了为贝叶斯框架的直接实现提供了原理，还在频率学家的框架内起着重要的作⽤，例如寻找最⼤似然解。特别地，对于EM算法中的E步骤⽆法解析地计算的模型，采样⽅法也可以⽤来近似E步骤。考虑⼀个模型，它的隐含变量为$\mathbf{Z}$，可见（观测）变量为$\mathbf{X}$，参数为$\mathbf{\theta }$。在M步骤中关于$\mathbf{\theta }$最⼤化的步骤为完整数据对数似然的期望，形式为$$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)=\int p\left(\mathbf{Z}|\mathbf{X},{\mathbf{\theta }}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)\ln p(\mathbf{Z},\mathbf{X}|\mathbf{\theta })\mathrm{d}\mathbf{Z}$$可以利用采样的方式去近似这个结果$$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{\text{ old }}\right)\simeq \frac{1}{L}\sum _{l=1}^L\ln p\left({\mathbf{Z}}^{(l)},\mathbf{X}|\mathbf{\theta }\right)$$然后在最大化上式，这个方法称为Monte Carlo EM algorithm。跟之前一样，这个还可以引入$\mathbf{\theta }$的先验，$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{\text{ old }}\right)+\ln p(\mathbf{\theta })$。
蒙特卡罗EM算法的⼀个特定的情形，被称为随机EM（stochastic EM）。如果我们考虑有限数量的概率分布组成的混合模型，并且在每个E步骤中只抽取⼀个样本时，我们就会⽤到这种算法。这⾥，潜在变量$\mathbf{Z}$描述了$K$个混合分量中的哪个分量被⽤于⽣成每个数据点。在E步骤中，$\mathbf{Z}$的样本从后验概率分布$p\left(\mathbf{Z}|\mathbf{X},{\mathbf{\theta }}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)$中抽取，其中$\mathbf{X}$是数据集。这⾼效地将每个数据点硬性地分配到混合分布中的⼀个分量中。在M步骤中，对于后验概率分布的这个采样的近似被⽤于按照平常的⽅式更新模型的参数。
现在假设我们从最⼤似然的⽅法转移到纯粹的贝叶斯⽅法，其中我们希望从参数向量$\mathbf{\theta }$上的后验概率分布中进⾏采样。原则上，我们希望从联合后验分布$p(\mathbf{\theta },\mathbf{Z}|\mathbf{X})$中抽取样本，但是我们假设这个计算⼗分困难。进⼀步地，我们假设从完整数据参数的后验概率分布$p(\mathbf{\theta }|\mathbf{Z},\mathbf{X})$中进⾏采样相对简单。这就产⽣了数据增⼴算法（data augmentation algorithm），它在两个步骤之间交替进⾏，这两个步骤被称为I步骤（归咎（imputation）步骤，类似于E步骤）和P步骤（后验（posterior）步骤，类似于M步骤）。

11.2 Markov Chain Monte Carlo

前⼀节中，我们讨论了计算函数期望的拒绝采样⽅法和重要采样⽅法，我们看到在⾼维空间中，这两种⽅法具有很⼤的局限性。因此，我们在本节中讨论⼀个⾮常⼀般的并且强⼤的框架，被称为马尔科夫链蒙特卡罗（Markov chain Monte Carlo, MCMC），它使得我们可以从⼀⼤类概率分布中进⾏采样，并且可以很好地应对样本空间维度的增长。
与拒绝采样和重要采样相同，我们再⼀次从提议分布中采样。但是这次我们记录下当前状态${\mathbf{z}}^{(\tau )}$，以及依赖于这个当前状态的提议分布$q\left(\mathbf{z}|{\mathbf{z}}^{(\tau )}\right)$，从⽽样本序列${\mathbf{z}}^{(1)},{\mathbf{z}}^{(2)},\dots$组成了⼀个马尔科夫链。与之前⼀样，如果我们有$p(\mathbf{z})=\tilde{p}(\mathbf{z})/Z_p$ ，那么我们会假定对于任意的$\mathbf{z}$值都可以计算$\tilde{p}(\mathbf{z})$，虽然$Z_p$的值可能未知。提议分布本⾝被选择为⾜够简单，从⽽直接采样很容易。在算法的每次迭代中，我们从提议分布中⽣成⼀个候选样本${\mathbf{z}}^{\star }$，然后根据⼀个恰当的准则接受这个样本。
在基本的Metropolis algorithm中，proposal分布为对称的$q\left({\mathbf{z}}_A|{\mathbf{z}}_B\right)=q\left({\mathbf{z}}_B|{\mathbf{z}}_A\right)$，则接受概率为$$A\left({\mathbf{z}}^{\star },{\mathbf{z}}^{(\tau )}\right)=\min \left(1,\frac{\tilde{p}\left({\mathbf{z}}^{\star }\right)}{\tilde{p}\left({\mathbf{z}}^{(\tau )}\right)}\right)$$当接受时，${\mathbf{z}}^{(\tau +1)}={\mathbf{z}}^{\star }$，否则${\mathbf{z}}^{(\tau +1)}={\mathbf{z}}^{\tau }$，然后再抽样$q\left(\mathbf{z}|{\mathbf{z}}^{(\tau +1)}\right)$，可以看出，当⼀个候选点被拒绝时，前⼀个样本点会被包含到最终的样本的列表中，从⽽产⽣了样本点的多个副本。正如我们将看到的那样，只要对于任意的${\mathbf{z}}_A$和${\mathbf{z}}_B$都有$q({\mathbf{z}}_A|{\mathbf{z}}_B)$为正（这是⼀个充分条件但不是必要条件），那么当$\tau \to \infty$时，${\mathbf{z}}^{(\tau )}$趋近于$p(\mathbf{z})$。然⽽，应该强调的是，序列${\mathbf{z}}^{(1)},{\mathbf{z}}^{(2)},\dots$不是来⾃$p(\mathbf{z})$的⼀组独⽴的样本，因为连续的样本是⾼度相关的。如果我们希望得到独⽴的样本，那么我们可以丢弃序列中的⼤部分样本，每$M$个样本中保留⼀个样本。对于充分⼤的$M$，保留的样本点对于所有的实际⽤途来说都是独⽴的。下图给出了⼀个简单的例⼦，这个例⼦使⽤Metropolis算法从⼀个⼆维⾼斯分布中采样，其中提议分布是⼀个各向同性的⾼斯分布。

通过考察⼀个具体的例⼦，即简单的随机游⾛的例⼦，我们可以对马尔科夫链蒙特卡罗算法的本质得到更深刻的认识。考虑⼀个由整数组成的状态空间$z$，概率为$$\begin{array}{rl}p\left(z^{(\tau +1)}=z^{(\tau )}\right)& =0.5\\ p\left(z^{(\tau +1)}=z^{(\tau )}+1\right)& =0.25\\ p\left(z^{(\tau +1)}=z^{(\tau )}-1\right)& =0.25\end{array}$$初始值为$z^{(1)}=0$，则$\mathbb{E}\left[z^{(\tau )}\right]=0$，$\mathbb{E}\left[{\left(z^{(\tau )}\right)}^2\right]=\tau /2$。因此，在$\tau$步骤之后，随机游⾛所经过的平均距离正⽐于$\tau$的平⽅根。这个平⽅根依赖关系是随机游⾛⾏为的⼀个典型性质，表明了随机游⾛在探索状态空间时是很低效的。正如我们会看到的那样，设计马尔科夫链蒙特卡罗⽅法的⼀个中⼼⽬标就是避免随机游⾛⾏为。

11.2.1 Markov chains

在详细讨论马尔科夫链蒙特卡罗⽅法之前，仔细研究马尔科夫链的⼀些⼀般的性质是很有⽤的。特别地，我们考察在什么情况下马尔科夫链会收敛到所求的概率分布上。一阶马尔科夫链$$p\left({\mathbf{z}}^{(m+1)}|{\mathbf{z}}^{(1)},\dots ,{\mathbf{z}}^{(m)}\right)=p\left({\mathbf{z}}^{(m+1)}|{\mathbf{z}}^{(m)}\right)$$我们可以按照下⾯的⽅式具体化⼀个马尔科夫链：给定初始变量的概率分布$p\left({\mathbf{z}}^{(0)}\right)$，以及后续变量的条件概率，⽤转移概率（transition probability）$T_m\left({\mathbf{z}}^{(m)},{\mathbf{z}}^{(m+1)}\right)\equiv p\left({\mathbf{z}}^{(m+1)}|{\mathbf{z}}^{(m)}\right)$的形式表⽰。如果对于所有的$m$，转移概率都相同，那么这个马尔科夫链被称为同质（homogeneous）。
边缘分布$$p\left({\mathbf{z}}^{(m+1)}\right)=\sum _{{\mathbf{z}}^{(m)}}p\left({\mathbf{z}}^{(m+1)}|{\mathbf{z}}^{(m)}\right)p\left({\mathbf{z}}^{(m)}\right)$$对于⼀个概率分布来说，如果马尔科夫链中的每⼀步都让这个概率分布保持不变，那么我们说这个概率分布关于这个马尔科夫链是不变的，或者静⽌的。因此，对于⼀个转移概率为$T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)$的同质的马尔科夫链来说，如果$$p^{\star }(\mathbf{z})=\sum _{{\mathbf{z}}^{\prime }}T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)p^{\star }\left({\mathbf{z}}^{\prime }\right)$$分布$p^{\star }(\mathbf{z})$就是不变的。注意，⼀个给定的马尔科夫链可能有多个不变的概率分布。例如，如果转移概率由恒等变换给出，那么任意的概率分布都是不变的。
确保所求的概率分布$p(\mathbf{z})$不变的⼀个充分（⾮必要）条件是令转移概率满⾜细节平衡（detailed balance）性质，定义为$$p^{\star }(\mathbf{z})T\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)=p^{\star }\left({\mathbf{z}}^{\prime }\right)T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)$$很容易看到，满⾜关于特定概率分布的细节平衡性质的转移概率会使得那个概率分布具有不变性$$\sum _{{\mathbf{z}}^{\prime }}{p}^{\star }\left({\mathbf{z}}^{\prime }\right)T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)=\sum _{{\mathbf{z}}^{\prime }}{p}^{\star }(\mathbf{z})T\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)=p^{\star }(\mathbf{z})\sum _{{\mathbf{z}}^{\prime }}p\left({\mathbf{z}}^{\prime }|\mathbf{z}\right)=p^{\star }(\mathbf{z})$$满⾜细节平衡性质的马尔科夫链被称为可翻转的（reversible）。我们的⽬标是使⽤马尔科夫链从⼀个给定的概率分布中采样。如果我们构造⼀个马尔科夫链使得所求的概率分布是不变的，那么我们就可以达到这个⽬标。然⽽，我们还要要求对于$m\to \infty$，概率分布$p\left({\mathbf{z}}^{(m)}\right)$收敛到所求的不变的概率分布$p^{\star }(\mathbf{z})$，与初始概率分布$p\left({\mathbf{z}}^{(0)}\right)$⽆关。这种性质被称为各态历经性（ergodicity），这个不变的概率分布被称为均衡（equilibrium）分布。很明显，⼀个具有各态历经性的马尔科夫链只能有唯⼀的⼀个均衡分布。可以证明，同质的马尔科夫链具有各态历经性，只需对不变的概率分布和转移概率做出较弱的限制即可（Neal, 1993）。
在实际中，我们经常可以从⼀组“基”转移$B_1,\dots ,B_K$中构建转移概率，⽅法为：将各个“基”转移表⽰为混合概率分布，形式为$$T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)=\sum _{k=1}^K{\alpha }_k{B}_k\left({\mathbf{z}}^{\prime },\mathbf{z}\right)$$mixing coefficients ${\alpha }_1,\dots ,{\alpha }_K$ satisfying ${\alpha }_k⩾0$ and ${\sum }_k{\alpha }_k=1$此外，基转移可以通过连续的应⽤组合到⼀起，$$T\left({\mathbf{z}}^{\prime },\mathbf{z}\right)=\sum _{{\mathbf{z}}_1}\dots \sum _{{\mathbf{z}}_{n-1}}{B}_1\left({\mathbf{z}}^{\prime },{\mathbf{z}}_1\right)\dots B_{K-1}\left({\mathbf{z}}_{K-2},{\mathbf{Z}}_{K-1}\right)B_K\left({\mathbf{z}}_{K-1},\mathbf{z}\right)$$如果⼀个概率分布关于每个基转移都是不变的，那么显然它关于公式（1）和（2）也是不变的。对于公式（1）的混合分布，如果每个基转移满⾜细节平衡，那么混合转移$T$也满⾜细节平衡。这对于使⽤公式（2）构造的转移概率不成⽴，虽然通过将基转移的顺序对称化，即采⽤$B_1,B_2,\dots ,B_K,B_K,\dots ,B_2,B_1$的形式，细节平衡的性质可以被恢复。使⽤组合转移概率的⼀个常见的例⼦是每个基转移只改变变量的⼀个⼦集的情形。

11.2.2 The Metropolis-Hastings algorithm

之前介绍了基本的Metropolis算法，现在说一个更加一般的方法，在当前得到状态${\mathbf{z}}^{(\tau )}$，从$q_k\left(\mathbf{z}|{\mathbf{z}}^{(\tau )}\right)$抽样出${\mathbf{z}}^{\star }$，则接受概率为$$A_k\left({\mathbf{z}}^{\star },{\mathbf{z}}^{(\tau )}\right)=\min \left(1,\frac{\tilde{p}\left({\mathbf{z}}^{\star }\right)q_k\left({\mathbf{z}}^{(\tau )}|{\mathbf{z}}^{\star }\right)}{\tilde{p}\left({\mathbf{z}}^{(\tau )}\right)q_k\left({\mathbf{z}}^{\star }|{\mathbf{z}}^{(\tau )}\right)}\right)$$与之前一样，也可不用归一化的$p(\mathbf{z})$。而Metropolis算法，取的是对称的proposal分布。下面将证明$p(\mathbf{z})$是构造出来的马尔科夫链的平稳分布。$$\begin{array}{rl}p(\mathbf{z})q_k\left({\mathbf{z}}^{\prime }|\mathbf{z}\right)A_k\left({\mathbf{z}}^{\prime },\mathbf{z}\right)& =\min \left(p(\mathbf{z})q_k\left({\mathbf{z}}^{\prime }|\mathbf{z}\right),p\left({\mathbf{z}}^{\prime }\right)q_k\left(\mathbf{z}|{\mathbf{z}}^{\prime }\right)\right)\\ & =\min \left(p({\mathbf{z}}^{\prime })q_k\left(\mathbf{z}|{\mathbf{z}}^{\prime }\right),p\left(\mathbf{z}\right)q_k\left({\mathbf{z}}^{\prime }|\mathbf{z}\right)\right)\\ & =p\left({\mathbf{z}}^{\prime }\right)q_k\left(\mathbf{z}|{\mathbf{z}}^{\prime }\right)A_k\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)\end{array}$$
提议分布的具体的选择会对算法的表现产⽣重要的影响。对于连续状态空间来说，⼀个常见的选择是⼀个以当前状态为中⼼的⾼斯分布，这会在确定分布的⽅差参数时需要进⾏⼀个重要的折中。如果⽅差过⼩，那么接受的转移的⽐例会很⾼，但是遍历状态空间的形式是⼀个缓慢的随机游⾛过程，导致较长的时间开销。然⽽，如果⽅差过⼤，那么拒绝率会很⾼，因为在我们考虑的这种复杂问题中，许多的步骤会到达$p(\mathbf{z})$很低的状态。考虑⼀个多元概率分布$p(\mathbf{z})$，它在$\mathbf{z}$的元素之间具有很强的相关性，如图所⽰。提议分布的标度$\rho$应该尽可能⼤，同时要避免达到较⾼的拒绝率。这表明，$\rho$应该与最⼩的长度标度${\sigma }_{\min }$是同⼀个量级的。然后，系统通过随机游⾛的⽅式探索伸长的⽅向，因此到达⼀个与原始状态或多或少独⽴的状态所需的步骤数量是${\left({\sigma }_{\max }/{\sigma }_{\min }\right)}^2$量级的。事实上，在⼆维的情形下，随着$\rho$的增加，拒绝率的增加会被接收的转移步骤数的增加所抵消。更⼀般地，对于多元⾼斯分布，得到独⽴样本所需的步骤的数量的增长量级是${\left({\sigma }_{\max }/{\sigma }_2\right)}^2$的，其中${\sigma }_2$是第⼆⼩的标准差（Neal, 1993）。抛开这些细节不谈，如果概率分布在不同的⽅向上的差异⾮常⼤，那么Metropolis-Hastings算法的收敛速度会⾮常慢。

11.3 Gibbs Sampling

我们需要从分布$p(\mathbf{z})=p\left(z_1,\dots ,z_M\right)$中采样，假定我们已经初始化了马尔科夫链，Gibbs采样是每次从条件分布$p\left(z_i|{\mathbf{z}}_{\backslash i}\right)$中采样代替$z_i$。
具体地，采样$p\left(z_1,z_2,z_3\right)$，在某步有$z_1^{(\tau )},z_2^{(\tau )}$ and $z_3^{(\tau )}$，从$p\left(z_1|z_2^{(\tau )},z_3^{(\tau )}\right)$中采出$z_1^{(\tau +1)}$替代$z_1^{(\tau )}$；然后从$p\left(z_2|z_1^{(\tau +1)},z_3^{(\tau )}\right)$采出$z_2^{(\tau +1)}$替代$z_2^{(\tau )}$；然后从$p\left(z_3|z_1^{(\tau +1)},z_2^{(\tau +1)}\right)$
采出$z_3^{(\tau +1)}$。

只需要证明符合细致平稳过程即可，很容易。还需要说明这个过程为各态历经性。各态历经性的⼀个充分条件是没有条件概率分布处处为零。如果这个要求满⾜，那么$\mathbf{z}$空间中的任意⼀点都可以从其他的任意⼀点经过有限步骤达到，这些步骤中每次对⼀个变量进⾏更新。如果这个要求没有满⾜，即某些条件概率分布为零，那么在这种情况下应⽤吉布斯采样时，必须显式地证明各态历经性。
为了完成算法，初始状态的概率分布也应该被指定，虽然在多轮迭代之后，样本与初始状态的分布⽆关。当然，马尔科夫链中的连续的样本是⾼度相关的，因此为了得到近似独⽴的样本，需要对序列进⾏下采样。
我们可以将吉布斯采样步骤看成Metropolis-Hastings算法的⼀个特定的情况，如下所述。取proposal分布为$q_k\left({\mathbf{z}}^{\star }|\mathbf{z}\right)=p\left(z_k^{\star }|{\mathbf{z}}_{\backslash k}\right)$，其中${\mathbf{z}}_{\backslash k}^{\star }={\mathbf{z}}_{\backslash k}$。此时接受率为$$A\left({\mathbf{z}}^{\star },\mathbf{z}\right)=\frac{p\left({\mathbf{z}}^{\star }\right)q_k\left(\mathbf{z}|{\mathbf{z}}^{\star }\right)}{p(\mathbf{z})q_k\left({\mathbf{z}}^{\star }|\mathbf{z}\right)}=\frac{p\left(z_k^{\star }|{\mathbf{z}}_k^{\star }\right)p\left({\mathbf{z}}_{\backslash k}^{\star }\right)p\left(z_k|{\mathbf{z}}_{\backslash k}^{\star }\right)}{p\left(z_k|{\mathbf{z}}_{\backslash k}\right)p\left({\mathbf{z}}_{\backslash k}\right)p\left(z_k^{\star }|{\mathbf{z}}_{\backslash k}\right)}=1$$此时接受率为1。与Metropolis算法⼀样，我们可以通过研究吉布斯采样算法在⾼斯分布上的应⽤，更深刻地认识算法的原理。考虑两个相关变量上的⼀个⾼斯分布，如图所⽰。这个⾼斯分布的条件概率分布的宽度为$l$，边缘概率分布的宽度为$L$。典型的步长由条件概率分布确定，从⽽量级为$l$。由于状态按照随机游⾛的⽅式进⾏转移，因此得到这个分布中的独⽴样本所需的步骤数量的量级为$(L/l{)}^2$。当然，如果⾼斯分布不是相关的，那么吉布斯采样的效率是最⾼的。对于这个简单的问题，我们可以将坐标系旋转，从⽽解除变量之间的相关关系。然⽽，在实际应⽤中，通常找到这种变换是不可⾏的。

⼀种减⼩吉布斯采样过程中的随机游⾛⾏为的⽅法被称为过松弛（over-relaxation）（Adler,1981）。在这种⽅法的最初的形式中，它被⽤于处理条件概率分布是⾼斯分布的情形，这种情形要⽐多元⾼斯分布更⼀般，因为诸如⾮⾼斯分布$p(z,y)\propto \exp \left(-z^2{y}^2\right)$具有⾼斯条件分布的形式。在吉布斯采样算法的每个步骤中，对于⼀个特定的分量$z_i$，条件概率分布具有均值${\mu }_i$和⽅差${\sigma }_i^2$。在过松弛框架中，$z_i$被替换为$$z_i^{\prime }={\mu }_i+\alpha \left(z_i-{\mu }_i\right)+{\sigma }_i{\left(1-{\alpha }^2\right)}^{1/2}\nu$$其中$\nu$为标准高斯分布，$-1<\alpha <1$，当$\alpha =0$时，此为标准的Gibbs采样，当$\alpha <0$时，步骤会偏向于与均值相反的⼀侧。这个步骤使得所求的概率分布具有不变性，因为如果$z_i$的均值为${\mu }_i$，⽅差为${\sigma }_i^2$ ，那么$z_i^{\prime }$也是。过松弛的效果是当变量⾼度相关时，⿎励在状态空间中的直接移动。有序过松（ordered over-relaxation）框架（Neal, 1999）将这种⽅法推⼴到了⾮⾼斯分布的情形。
具体运用Gibbs采样在图模型中，都是从$p\left(z_k|{\mathbf{z}}_{\backslash k}\right)$采样，对于剩下的节点主要就是马尔科夫毯，如图所示

对于有向图来说，以某个结点的⽗结点为条件，这个结点的⼀⼤类条件概率分布都会使得⽤于吉布斯采样的概率分布是对数凹函数。于是，11.1.3节讨论的可调节拒绝采样⽅法提供了有向图的蒙特卡罗采样⽅法的⼀个框架，这种⽅法具有⼴泛的适⽤性。
如果图是使⽤指数族分布构建的，并且⽗结点-⼦结点关系保持共轭，那么吉布斯采样中的完整的条件概率分布会与定义在每个结点的原始的条件概率分布（以⽗结点为条件）具有相同的函数形式，因此可以使⽤标准的采样⽅法。通常，完整的条件概率分布的形式会很复杂，从⽽⽆法使⽤标准的采样⽅法。然⽽，如果这些条件概率分布是对数凹函数，那么使⽤可调整的拒绝采样⽅法，采样可以⾼效地完成（假设对应的变量是标量）。
如果在吉布斯采样算法的每个阶段，我们不从对应的条件概率分布中抽取样本，⽽是对变量进⾏⼀个点估计，这个点估计由条件概率分布的最⼤值给出，那么我们就得到了8.3.3节讨论的迭代条件峰值（ICM）算法。因此，ICM可以看成是吉布斯采样的⼀种贪⼼近似。
由于基本的吉布斯采样⽅法每次只考虑⼀个变量，因此它在连续样本之间具有很强的依赖性。在另⼀个极端情况下，如果我们直接从联合概率分布中采样（我们⼀直假定这种操作⽆法完成），那么连续的样本点之间就是独⽴的。我们可以采⽤⼀种折中的⽅法来提升简单的吉布斯采样的效果，即我们连续地对⼀组变量进⾏采样，⽽不是对⼀个变量进⾏采样。这就是分块吉布斯（blocking Gibbs）采样算法。这种算法中，将变量集合分块（未必互斥），然后在每个块内部联合地采样，采样时以剩余的变量为条件（Jensen et al., 1995）。

11.4 Slice Sampling

我们已经看到，Metropolis算法的⼀个困难之处是它对于步长的敏感性。如果步长过⼩，那么由于随机游⾛⾏为，算法会很慢。⽽如果步长过⼤，那么由于较⾼的拒绝率，算法会很低效。切⽚采样（slice sampling）⽅法（Neal, 2003）提供了⼀个可以⾃动调节步长来匹配分布特征的⽅法。与之前⼀样，它需要我们能够计算未归⼀化的概率分布$\tilde{p}(\mathbf{z})$。
首先考虑单变量的情况，切片采样对$z$进行增广，即从联合分布$(z,u)$中采样。之后hybrid Monte Carlo也是采用这种形式。从该分布空间中均匀采样$$\hat{p}(z,u)=\{\begin{array}{ll}1/Z_p& \text{ if }0⩽u⩽\tilde{p}(z)\\ 0& \text{ otherwise }\end{array}$$其中$Z_p=\int \tilde{p}(z)\mathrm{d}z$，则边缘分布为$$\int \hat{p}(z,u)\mathrm{d}u={\int }_0^{\tilde{p}(z)}\frac{1}{Z_p}\mathrm{d}u=\frac{\tilde{p}(z)}{Z_p}=p(z)$$则我们可以从$\hat{p}(z,u)$采出$z,u$并且忽略$u$，则得到了$p(z)$采出的样本。具体来说，可以迭代地采样$z,u$，首先给定$z$，得到$\tilde{p}(z)$，然后从$[0,\tilde{p}(z)]$中均匀采出$u$；然后，我们固定$u$，在由 { z : p ~ ( z ) &gt; u } \{z : \widetilde{p}(z)&gt;u\} {z:p ​(z)>u}定义的分布的“切⽚”上，对$z$进⾏均匀地采样。如下图(a)所示。

在实际应⽤中，直接从穿过概率分布的切⽚中采样很困难，因此我们定义了⼀个采样⽅法，它保持$\hat{p}(z,u)$下的均匀分布具有不变性，这可以通过确保满⾜细节平衡的条件来实现。假设$z$的当前值记作$z^{(\tau )}$，并且我们已经得到了⼀个对应的样本$u$。$z$的下⼀个值可以通过考察包含$z^{(\tau )}$的区域$z_{\min }⩽z⩽z_{\max }$来获得。根据概率分布的特征长度标度来对步长进⾏的调节就发⽣在这⾥。我们希望区域包含尽可能多的切⽚，从⽽使得$z$空间中能进⾏较⼤的移动，同时希望切⽚外的区域尽可能⼩，因为切⽚外的区域会使得采样变得低效。
⼀种选择区域的⽅法是，从⼀个包含$z^{(\tau )}$的具有某个宽度$w$的区域开始，然后测试每个端点，看它们是否位于切⽚内部。如果有端点没在切⽚内部，那么区域在增加$w$值的⽅向上进⾏扩展，直到端点位于区域外。然后，$z^{\prime }$的⼀个样本被从这个区域中均匀抽取。如果它位于切⽚内，那么它就构成了$z^{(\tau +1)}$。如果它位于切⽚外，那么区域收缩，使得$z^{\prime }$组成⼀个端点，并且区域仍然包含$z^{(\tau )}$。然后，另⼀个样本点从这个缩⼩的区域中均匀抽取，以此类推，直到找到位于切⽚内部的⼀个$z$值。
切⽚采样可以应⽤于多元分布中，⽅法是按照吉布斯采样的⽅式重复地对每个变量进⾏采样。这要求对于每个元素$z_i$，我们能够计算⼀个正⽐于$p\left(z_i|{\mathbf{z}}_{\backslash i}\right)$的函数。

11.5 The Hybrid Monte Carlo Algorithm

正如我们已经注意到的那样，Metropolis算法的⼀个主要的局限是它具有随机游⾛的⾏为，⽽在状态空间中遍历的距离与步骤数量只是平⽅根的关系。仅仅通过增⼤步长的⽅式是⽆法解决这个问题的，因为这会使得拒绝率变⾼。
本节中，我们介绍⼀类更加复杂的转移⽅法。这些⽅法基于对物理系统的⼀个类⽐，能够让系统状态发⽣较⼤的改变，同时让拒绝的概率较低。它适⽤于连续变量上的概率分布，对于连续变量，我们已经能够计算对数概率关于状态变量的梯度。我们会在11.5.1节讨论动态系统框架，然后在11.5.2节，我们会解释这个框架如何与Metropolis算法结合，产⽣出⼀个强⼤的混合蒙特卡罗算法。这⾥不需要物理学的背景，因为本节是⾃洽的，并且关键的结果全部从基本的原理中推导出。

11.5.1 Dynamical systems

我们考虑的动⼒学对应于在连续时刻（记作$\tau$）下的状态变量 z = { z i } \mathbf{z}=\left\{z_{i}\right\} z={zi​}的演化。经典的动⼒学由⽜顿第⼆定律描述，即物体的加速度正⽐于施加的⼒，对应于关于时间的⼆阶微分⽅程。我们可以将⼀个⼆阶微分⽅程分解为两个相互偶合的⼀阶⽅程，⽅法是引⼊中间的动量（momentum）变量$\mathbf{r}$，对应于状态变量$\mathbf{z}$的变化率，元素为$$r_i=\frac{\mathrm{d}{z}_i}{\mathrm{d}\tau }$$其中$z_i$看作位置变量，对于每个位置有个动量变量，这两个变量组成了相空间(phase space)。不失一般性，我们可以将概率分布$p(\mathbf{z})$写成下⾯的形式$$p(\mathbf{z})=\frac{1}{Z_p}\exp (-E(\mathbf{z}))$$其中$E(\mathbf{z})$称为势能量，系统的加速度是动量的变化率，通过施
加⼒（force）的⽅式确定，它本⾝是势能的负梯度，即$$\frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }=-\frac{\partial E(\mathbf{z})}{\partial z_i}$$使⽤哈密顿框架重新写出这个动态系统的公式是⽐较⽅便的。为了完成这⼀点，我们⾸先将动能（kinetic energy）定义为$$K(\mathbf{r})=\frac{1}{2}\Vert \mathbf{r}{\Vert }^2=\frac{1}{2}\sum _i{r}_i^2$$系统的总能量是势能和动能之和，即$$H(\mathbf{z},\mathbf{r})=E(\mathbf{z})+K(\mathbf{r})$$其中$H$为汉密尔顿函数，我们现在可以将系统的动⼒学⽤哈密顿⽅程的形式表⽰出来，形式为$$\begin{array}{rl}\frac{\mathrm{d}{z}_i}{\mathrm{d}\tau }& =\frac{\partial H}{\partial r_i}\\ \frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }& =-\frac{\partial H}{\partial z_i}\end{array}$$在动态系统的变化过程中，哈密顿函数$H$的值是⼀个常数，这⼀点通过求微分的⽅式很容易看出来。$$\begin{array}{rl}\frac{\mathrm{d}H}{\mathrm{d}\tau }& =\sum _i\left\{\frac{\partial H}{\partial z_i}\frac{\mathrm{d}{z}_i}{\mathrm{d}\tau }+\frac{\partial H}{\partial r_i}\frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }\right\}\\ & =\sum _i\left\{\frac{\partial H}{\partial z_i}\frac{\partial H}{\partial r_i}-\frac{\partial H}{\partial r_i}\frac{\partial H}{\partial z_i}\right\}=0\end{array}$$哈密顿动态系统的第⼆个重要性质是动态系统在相空间中体积不变，这被称为Liouville定理（Liouville’s Theorem）。换句话说，如果我们考虑变量$(\mathbf{z},\mathbf{r})$空间中的⼀个区域，那么当这个区域在哈密顿动态⽅程下的变化时，它的形状可能会改变，但是它的体积不会改变。可以这样证明：我们注意到流场（位置在相空间的变化率）为$$\mathbf{V}=\left(\frac{\mathrm{d}\mathbf{z}}{\mathrm{d}\tau },\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}\tau }\right)$$这个场的散度为零，即$$\begin{array}{rl}\mathrm{div}\mathbf{V}& =\sum _i\left\{\frac{\partial }{\partial z_i}\frac{\mathrm{d}{z}_i}{\mathrm{d}\tau }+\frac{\partial }{\partial r_i}\frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }\right\}\\ & =\sum _i\left\{-\frac{\partial }{\partial z_i}\frac{\partial H}{\partial r_i}+\frac{\partial }{\partial r_i}\frac{\partial H}{\partial z_i}\right\}=0\end{array}$$现在考虑相空间上的联合概率分布，它的总能量是哈密顿函数，即概率分布的形式为$$p(\mathbf{z},\mathbf{r})=\frac{1}{Z_H}\exp (-H(\mathbf{z},\mathbf{r}))$$使⽤体系的不变性和$H$的守恒性，可以看到哈密顿动态系统会使得$p(\mathbf{z},\mathbf{r})$保持不变。可以这样证明：考虑相空间的⼀个⼩区域，区域中$H$近似为常数。如果我们跟踪⼀段有限时间内的哈密顿⽅程的变化，那么这个区域的体积不会发⽣改变，从⽽这个区域的$H$的值不会发⽣改变，因此概率密度（只是$H$的函数）也不会改变。
虽然$H$是不变的，但是$\mathbf{z}$和$\mathbf{r}$会发⽣变换，因此通过在⼀个有限的时间间隔上对哈密顿动态系统积分，我们就可以让$\mathbf{z}$以⼀种系统化的⽅式发⽣较⼤的变化，避免了随机游⾛的⾏为。
然⽽，哈密顿动态系统的变化对$p(\mathbf{z},\mathbf{r})$的采样不具有各态历经性，因为$H$的值是⼀个常数。为了得到⼀个具有各态历经性的采样⽅法，我们可以在相空间中引⼊额外的移动，这些移动会改变$H$的值，同时也保持了概率分布$p(\mathbf{z},\mathbf{r})$的不变性。达到这个⽬标的最简单的⽅式是将$\mathbf{r}$的值替换为⼀个从以$\mathbf{z}$为条件的概率分布中抽取的样本。这可以被看成吉布斯采样的步骤，因此根据11.3节，我们看到这也使得所求的概率分布保持了不变性。注意，$\mathbf{z}$和$\mathbf{r}$在概率分布$p(\mathbf{z},\mathbf{r})$中是独⽴（根据定义可以分解）的，我们看到条件概率分布$p(\mathbf{r}|\mathbf{z})$是⾼斯分布，从中我们可以很容易地进⾏采样。
在这种⽅法的⼀个实际应⽤中，我们必须解决计算哈密顿⽅程的数值积分的问题。这会引⼊⼀些数值的误差，因此我们要设计⼀种⽅法来最⼩化这些误差产⽣的影响。事实上，可以证明，能够在Liouville定理仍然精确成⽴的条件下，对积分⽅法进⾏修改。这个性质在11.5.2节讨论混合蒙特卡罗算法时很重要。完成这件事的⼀种⽅法是蛙跳（leapfrog）离散化。这种⽅法使⽤下⾯的公式对位置变量和动量变量的离散时间近似$\hat{\mathbf{z}}$ and $\hat{\mathbf{r}}$进⾏交替地更新。$$\begin{array}{rl}{\hat{r}}_i(\tau +ϵ/2)& ={\hat{r}}_i(\tau )-\frac{ϵ}{2}\frac{\partial E}{\partial z_i}(\hat{\mathbf{z}}(\tau ))\\ {\hat{z}}_i(\tau +ϵ)& ={\hat{z}}_i(\tau )+ϵ{\hat{r}}_i(\tau +ϵ/2)\\ {\hat{r}}_i(\tau +ϵ)& ={\hat{r}}_i(\tau +ϵ/2)-\frac{ϵ}{2}\frac{\partial E}{\partial z_i}(\hat{\mathbf{z}}(\tau +ϵ))\end{array}$$我们看到，这种⽅法对动量变量的更新形式是半步更新，步长为$ϵ/2$，接着是对位置变量的整步更新，步长为$ϵ$，然后是对动量变量的第⼆个半步更新。如果我们连续地使⽤⼏次蛙跳，那么可以看到，对动量变量的半步更新可以结合到步长为$ϵ$的整步更新中。于是，位置变量的更新和动量变量的更新互相之间以蛙跳的形式结合。为了将动态系统推进进⼀个时间间隔$\tau$，我们需要进⾏$\tau /ϵ$个步骤。对连续时间动态系统的离散化近似引⼊的误差会在极限$ϵ\to 0$的情况下趋于零，假设函数$E(\mathbf{z})$是光滑的。然⽽，对于实际应⽤中使⽤的⼀个⾮零的$ϵ$，⼀些保留的误差仍然会存在。我们会在11.5.2节看到在混合蒙特卡罗算法中，这些误差的影响如何被消除。
总结⼀下，哈密顿动⼒学⽅法涉及到交替地进⾏⼀系列蛙跳更新以及根据动量变量的边缘分布进⾏重新采样。
注意，与基本的Metropolis⽅法不同，哈密顿动⼒学⽅法能够利⽤对数概率分布的梯度信息以及概率分布本⾝的信息。在函数最优化领域有⼀个类似的情形。⼤多数可以得到梯度信息的情况下，使⽤哈密顿动⼒学⽅法是很有优势的。⾮形式化地说，这种现象是由于下⾯的事实造成的：在$D$维空间中，与计算函数本⾝的代价相⽐，计算梯度所带来的额外的计算代价通常是⼀个与$D$⽆关的固定因⼦。⽽与函数本⾝只能传递⼀条信息相⽐，$D$维梯度向量可以传递$D$条信息。

11.5.2 Hybrid Monte Carlo

正如我们在前⼀节讨论的那样，对于⼀个⾮零的步长$ϵ$，蛙跳算法的离散化会在哈密顿动⼒学⽅程的积分过程中引⼊误差。混合蒙特卡罗（hybrid Monte Carlo）（Duane et al., 1987; Neal,1996）将哈密顿动态系统与Metropolis算法结合在⼀起，因此消除了与离散化过程关联的任何偏差。
具体来说，算法使⽤了⼀个马尔科夫链，它由对动量变量$\mathbf{r}$的随机更新以及使⽤蛙跳算法对哈密顿动态系统的更新交替组成。在每次应⽤蛙跳算法之后，基于哈密顿函数$H$的值，确定Metropolis准则，确定⽣成的候选状态被接受或者拒绝。因此，如果$(\mathbf{z},\mathbf{r})$是初始状态，$({\mathbf{z}}^{\star },{\mathbf{r}}^{\star })$是蛙跳积分后的状态，那么候选状态被接受的概率为 min ⁡ ( 1 , exp ⁡ { H ( z , r ) − H ( z ⋆ , r ⋆ ) } ) \min(1,\exp\{ H(\mathbf z,\mathbf r)-H(\mathbf z^{\star},\mathbf r^{\star})\}) min(1,exp{H(z,r)−H(z⋆,r⋆)})如果蛙跳积分完美地模拟了哈密顿动态系统，那么每个这种候选状态都会⾃动地被接受，因为$H$的值会保持不变。由于数值误差，$H$的值有时可能会减⼩，因此我们希望Metropolis准则将这种效果引发的任何偏差都消除，并且确保得到的样本确实是从所需的概率分布中抽取的。为了完成这件事，我们需要确保对应于蛙跳积分的更新⽅程满⾜细节平衡。通过按照下⾯的⽅式修改蛙跳⽅法，这个⽬标很容易实现。
在开始蛙跳积分序列之前，我们等概率地随机选择是沿着时间向前的⽅向积分（步长为$ϵ$）还是沿着时间向后的⽅向积分（步长为$-ϵ$）。我们⾸先注意到，蛙跳积分⽅法（三个步骤，如上一节）是时间可翻转的，即$L$步使⽤步长为$-ϵ$的积分会抵消$L$步使⽤步长为$ϵ$的积分。接下来我们证明蛙跳积分精确地保持了相空间的体积不变性。这是因为，蛙跳⽅法中的每⼀步对$z_i$或者$r_i$的更新都只是另⼀个变量的函数。如图所⽰，这个现象产⽣的效果是将相空间的⼀个区域进⾏形变⽽不改变它的体积。

11.6 Estimating the Partition Function

$$p_E(\mathbf{z})=\frac{1}{Z_E}\exp (-E(\mathbf{z}))$$
之前我们都没考虑归一化常数，因为采样并不需要，但是这一项对于模型对比十分重要，也称为model evidence，因此十分有必要得到该值，一般来说直接求和积分是不现实的。
对于模型⽐较来说，我们所需的实际是两个模型的划分函数的⽐值。将这个⽐值与先验概率的⽐值相乘可以得到后验概率的⽐值。之后可以⽤这个⽐值来进⾏模型选择或者模型平均。
⼀种估计划分函数⽐值的⽅法是使⽤概率分布的重要采样，这个概率分布的能量函数为$G(\mathbf{z})$，即$$\begin{array}{rl}\frac{Z_E}{Z_G}& =\frac{{\sum }_{\mathbf{z}}\exp (-E(\mathbf{z}))}{{\sum }_{\mathbf{z}}\exp (-G(\mathbf{z}))}\\ & =\frac{{\sum }_{\mathbf{z}}^{}\exp (-E(\mathbf{z})+G(\mathbf{z}))\exp (G(\mathbf{z}))}{{\sum }_{\mathbf{z}}\exp (-G(\mathbf{z}))}\\ & ={\mathbb{E}}_{G(\mathbf{z})}[\exp (-E(\mathbf{z})+G(\mathbf{z}))]\\ & \simeq \sum _l\exp (-E({\mathbf{z}}^{(l)})+G({\mathbf{z}}^{(l)}))\end{array}$$其中 { z ( l ) } \{\mathbf z^{(l)}\} {z(l)}采样于$p_G(\mathbf{z})$，如果$p_G(\mathbf{z})$为一个简单的分布，那么可以将$Z_G$得到，从而得到$Z_E$。
这个方法的前提也是$p_G$能够很接近$p_E$，即⽐值$p_E/p_G$变化不⼤，那么这种⽅法会产⽣准确的结果。在实际应⽤中，对于本书中考察的复杂的模型，我们⽆法找到⼀个可以很容易地解析计算的重要采样分布。
于是，另⼀种⽅法是使⽤从马尔科夫链中得到的样本来定义重要采样分布。如果马尔科夫链的转移概率为$T(\mathbf{z},{\mathbf{z}}^{\prime })$，样本集合为${\mathbf{z}}^{(1)},\dots ,{\mathbf{z}}^{(L)}$，那么采样分布可以写成$$\frac{1}{Z_G}\exp (-G(\mathbf{z}))=\sum _{l=1}^{L}T({\mathbf{z}}^{(l)},\mathbf{z})$$然后将上式运用到刚刚提到的地方。