13 MCMC（Markov Chain Monte Carlo）

0 MCMC思想
1 采样方法
- 1.1概率分布采样
- 1.2 拒绝采样(Rejection Sampling)
- 1.3 重要性采样(Importance Sampling)
- 1.4 重要性重采样(Sampling Importance Resampling)
2 Markov Chain
- 2.1 基础概念介绍
- 2.2 平稳分布(Stationary Distribution)
3 MH采样（Metropolis Hastings Sampling）
- 3.1 Proposal Matrix
- 3.2 Metropolis-Hastings Sampling
- 3.3 小结
4 Gibbs Sampling（吉布斯采样）
- 4.1 A Example
- 4.2 接受率 α\alphaα 的计算
5 采样
- 5.1 采样的动机
- 5.2 什么样的是好样本
- 5.3 实际采样中的困难
- 5.4 采样方法
6 Method of MCMC
- 6.1 Markov Chain 收敛性介绍
- - 6.1.1 Markov Chain 状态转移计算
  - 6.1.2 Markov Chain 收剑性
- 6.2 Existing Problem
7 参考

0 MCMC思想

蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：

如果我们很难求解出f(x)的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢？假设我们函数图像如下图:

则一个简单的近似求解方法是在 [a,b][a,b][a,b] 之间随机的采样一个点。比如 x0x_0x0，然后用 f(x0)f(x_0)f(x0) 代表在 [a,b][a,b][a,b] 区间上所有的 f(x)f(x)f(x) 的值。那么上面的定积分的近似求解为:

当然，用一个值代表 [a,b][a,b][a,b] 区间上所有的 f(x)f(x)f(x) 的值，这个假设太粗糙。那么我们可以采样 [a,b][a,b][a,b] 区间的n个值： x0,x1,…,xn−1x_0,x_1,\dots,x_{n-1}x0,x1,…,xn−1 ，用它们的均值来代表 [a,b][a,b][a,b] 区间上所有的 f(x)f(x)f(x) 的值。这样我们上面的定积分的近似求解为:

虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即x在 [a,b][a,b][a,b] 之间是均匀分布的，而绝大部分情况，x在 [a,b][a,b][a,b] 之间不是均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远。

怎么解决这个问题呢？

如果我们可以得到 x 在 [a,b][a,b][a,b] 的概率分布函数 p(x)，那么我们的定积分求和可以这样进行：

上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。

可以看出，最上面我们假设 x 在 [a,b][a,b][a,b] 之间是均匀分布的时候，p(xi)=1/(b−a)p(x_i)=1/(b-a)p(xi)=1/(b−a)，带入我们有概率分布的蒙特卡罗积分的上式，可以得到：

也就是说，我们最上面的均匀分布也可以作为一般概率分布函数p(x)在均匀分布时候的特例。那么我们现在的问题转到了如何求出x的分布p(x)的若干和样本上来。

其实在之前的 Inference Variational 那一节中, 我们讲到过一些有关于 Markov Chain Monte Carlo (MCMC) 的知识。也就是我们有一些数据 X，看到这些数据 X，并且有一些隐变量 Z，我们给隐变量一些先验，根据观测数据来推后验知识，也就是 P(Z∣X)P(Z | X)P(Z∣X)

但是，很不幸的是 P(Z∣X)P(Z | X)P(Z∣X) 的计算非常的复杂，我们大致采用两种思路来解决这个问题，也就是 精确推断和近似推断。精确推断无法达到我们想要的结果时，就会采用近似推断的方法。而近似推断中我们又可以分成两大类，即为确定性近似 (VI) 和随机近似 (MCMC)。

Monte Carlo Method 是一种基于采样的随机近似算法。我们的目标是求解后验概率 P(Z∣X)P(Z | X)P(Z∣X) , 其中 Z 为 Latent data，X 为 Observed data。知道分布以后，我们通常的目标是求解：
EZ∣X[f(Z)]=∫ZP(Z∣X)f(Z)dZ≈1N∑i=1Nf(zi)\mathbb{E}_{Z | X}[f(Z)]=\int_{Z} P(Z | X) f(Z) d Z \approx \frac{1}{N} \sum_{i=1}^{N} f\left(z_{i}\right)EZ∣X[f(Z)]=∫ZP(Z∣X)f(Z)dZ≈N1i=1∑Nf(zi)然后，问题马上就来了，我们知道了后验分布 P(Z∣X)P(Z | X)P(Z∣X)，怎么去采样呢? 也就是如何通过采样得到 z(1),z(2),⋯,z(N)∼P(Z∣X)z^{(1)}, z^{(2)}, \cdots, z^{(N)} \sim P(Z | X)z(1),z(2),⋯,z(N)∼P(Z∣X)。那么，我们这一节将要主要介绍三种采样方法：

概率分布采样
拒绝采样
重要性采样

所以，MCMC要解决的问题是：

计算复杂积分 ∫p(x)q(x)q(x)dx\int \frac{p(x)}{q(x)} q(x)dx∫q(x)p(x)q(x)dx
求 f(x)f(x)f(x) 的期望，其中 x∼P(x)x \sim P(x)x∼P(x)
Ep(x)[f(x)]=∫P(x)f(x)dx\mathbb{E}_{p(x)}[f(x)]=\int P(x) f(x) d xEp(x)[f(x)]=∫P(x)f(x)dx

MCMC方法描述：
将所需求的式子描述成关于某个分布的期望的形式： Ep(x)[f(x)]\mathbb{E}_{p(x)}[f(x)]Ep(x)[f(x)]然后对分布 P(x)P(x)P(x) 进行采样，得到 x(1),x(2),⋯,x(N)∼P(x)x^{(1)}, x^{(2)}, \cdots, x^{(N)} \sim P(x)x(1),x(2),⋯,x(N)∼P(x) ,则有
Ep(x)[f(x)]≈1N∑i=1Nf(xi)\mathbb{E}_{p(x)}[f(x)] \approx \frac{1}{N} \sum_{i=1}^{N} f\left(x_{i}\right)Ep(x)[f(x)]≈N1i=1∑Nf(xi)

1 采样方法

1.1概率分布采样

常用写法

c.d.f：概率分布函数（从0到1的函数）
p.d.f ：概率密度函数
i.i.d ：独立同分布

为什么要有概率分布采样呢？因为我们直接根据概率分布来进行采样非常的复杂。如果我们知道概率分布的具体形式呢？我们可以直接求得概率累积的概率分布函数。

假设有一个概率密度函数（pdf） P(Z)P(Z)P(Z) ，在这个 pdf 上采样，通常需要先通过 pdf 求出其累积分布函数（cdf）。由于概率（累积）分布函数（cdf）的值一定是 [0,1] 之间的。所以，我们可以在均匀概率密度分布 U(0,1)U(0,1)U(0,1) 上采样, 得到 u(i)∼U(0,1)u^{(i)} \sim U(0,1)u(i)∼U(0,1)。然后通过 cdf 的反函数求 x(i)∼cdf−1(u(i))x(i)∼ c d f^{-1}\left(u^{(i)}\right)x(i)∼cdf−1(u(i)) 就可以计算得到我们想要的结果。这样就可以采样得到 {x(1),x(2),⋯,x(N)}N\left\{x^{(1)}, x^{(2)}, \cdots, x^{(N)}\right\} \mathrm{N}{x(1),x(2),⋯,x(N)}N 个样本点。当然，理论上这个方法好像很有效，但是实际上很多情况我们都根本不知道 p.d.f 的具体表现形式。就算知道，很多时候 c.d.f 也并不是那么的好求。所以很多情况下，概率分布采样并没有那么的好求。

为什么要按照 cdf 进行采样？

1.2 拒绝采样(Rejection Sampling)

由于对目标分布 p(Z)p(Z)p(Z) 的采样非常的困难，所以我们可以对一个比较筒单的分布 q(Z)q(Z)q(Z) 进行采样来辅助采样。那么我们具体做法怎么办呢？

我们可以设定一个 proposal distribution： q(Z)q(Z)q(Z) 。对于 ∀zi\forall z_{i}∀zi , 保证 M⋅q(zi)≥p(zi)M \cdot q\left(z^{i}\right) \geq p\left(z^{i}\right)M⋅q(zi)≥p(zi) , 那么我们为什么要引入 MMM 呢? 这是因为 ∫ZP(Z)dZ=∫Zq(Z)dZ=1\int_{Z} P(Z) d Z=\int_{Z} q(Z) d Z=1∫ZP(Z)dZ=∫Zq(Z)dZ=1。要使 q(zi)≥p(zi)q\left(z^{i}\right) \geq p\left(z^{i}\right)q(zi)≥p(zi) 是几乎不可能成立的。

为什么这样做呢？如下图，p(z)p(z)p(z) 是我们需要采样的原始分布，但是通常很难采样，或者说分布不清楚。所以我们可以⽤⼀个提议分布来进⾏采样，主要的思路是如果采样点落在 p(z)p(z)p(z) 下⽅，即蓝线下方的阴影部分则接收，反之拒绝。为了方便描述，我们画图来说明一下：

在这里我们需要定义一个接受率： α=P(z(l))M⋅q(z(t))\alpha=\frac{P\left(z^{(l)}\right)}{M \cdot q\left(z^{(t)}\right)}α=M⋅q(z(t))P(z(l)), 很显然 0≤α≤10 \leq \alpha \leq 10≤α≤1。这个实际就是上图中绿色的部分。我们来看看具体的步骤：

(1) 首先进行采样 z(i)∼q(z)z^{(i)} \sim q(z)z(i)∼q(z)
(2) 计算接受率： α=P(z(i))M⋅q(z(i))\alpha=\frac{P\left(z^{(i)}\right)}{M \cdot q\left(z^{(i)}\right)}α=M⋅q(z(i))P(z(i)),
(3) u∼U(0,1)u \sim U(0,1)u∼U(0,1); 如果 u≤αu \leq \alphau≤α, 我们就接收 z(i)z^{(i)}z(i) , 不然我们就拒绝。

所以，绿色的部分就被我们称为拒绝区域，就是这样来的，所以这个采样方法就是拒绝采样。

同样这样的采样方法也有缺点。如果 M⋅q(z)M \cdot q(z)M⋅q(z) 比 p(z)p(z)p(z) 大很多的话，那么我们的采样老是是失败的，这就涉及到一个采样效率低下的问题。而当 M⋅q(z)=p(z)M \cdot q(z)=p(z)M⋅q(z)=p(z) 的时候， α=1\alpha=1α=1, 我们每次采样的结果都是接受的。但是，实际上 p(z)p(z)p(z) 的分布形式非常的复杂，我们根本就没有办法来得到那么准确的结果，特别是采样 cost 非常高的话，经常性的采样失败带来的损失是很大的。

1.3 重要性采样(Importance Sampling)

拒绝采样中的接受率过小的时候，被丢弃的样本过多，导致我们得到的有效样本过少。重要性采样就是为了解决拒绝采样样本过少的问题。

重要性采样在我们的强化学习 (PPO) 中的应用非常的多。重要性采样并不是直接对概率进行采样，而是对概率分布的期望进行采样。也就是：
Ep(z)[f(z)]=∫p(z)f(z)dz=∫p(z)q(z)q(z)f(z)dz=∫f(z)p(z)q(z)q(z)dz≈1N∑i=1Nf(zi)p(zi)q(zi)zi∼q(z),i=1,2,⋯,N\begin{aligned} \mathbb{E}_{p(z)}[f(z)]=\int p(z) f(z) d z=& \int \frac{p(z)}{q(z)} q(z) f(z) d z \\ =& \int f(z) \frac{p(z)}{q(z)} q(z) d z \\ \approx & \frac{1}{N} \sum_{i=1}^{N} f\left(z_{i}\right) \frac{p\left(z_{i}\right)}{q\left(z_{i}\right)} \\ & z_{i} \sim q(z), i=1,2, \cdots, N \end{aligned}Ep(z)[f(z)]=∫p(z)f(z)dz==≈∫q(z)p(z)q(z)f(z)dz∫f(z)q(z)p(z)q(z)dzN1i=1∑Nf(zi)q(zi)p(zi)zi∼q(z),i=1,2,⋯,N 而这里的 p(zt)q(zt)\frac{p\left(z_{t}\right)}{q\left(z_{t}\right)}q(zt)p(zt) 也就是 Weight，用来平衡不同的概率密度值之间的差距，q(z)q(z)q(z)与p(z)p(z)p(z)越相似，则采样的效果越好。同样重要性采样也可能会出现一些问题，就是两个分布之间的差距太大了话，总是采样采不到重要的样本，采的可能都是实际分布概率值小的部分。重要性采样效率可能很低，可能出现其中一个样本的权重非常非常高，而其他的样本权重非常低，也就是采样效率不均匀的问题。

当使用重要性采样得到样本后，第二阶段是在 zi,…,zNz_i,\dots,z_Nzi,…,zN 上采样，以它们的权重为概率分布，即权重越大的概率越大，这样产生出的新的样本，权重大的数量越多。

在这个基础上，我们进一步提出了Sampling Importance Resampling。

没懂

1.4 重要性重采样(Sampling Importance Resampling)

经过重要性采样后，我们得到了 N 个样本点，以及对应的权重。那么我用权重来作为采样的概率，
重新测采样出 N 个样本。也就是如下图所示：

通过二次采样可以降低采样不平衡的问题。至于为什么呢？大家想一想，我在这里表达一下自己
的看法。 p(zt)q(zt)\frac{p\left(z_{t}\right)}{q\left(z_{t}\right)}q(zt)p(zt) 是Weight，如果Weight 比较大的话，说明 p(zi)p(z_i)p(zi) 比较大而 q(zi)q(z_i)q(zi) 比较的小，也就是我们通过 q(zi)q(z_i)q(zi) 采出来的数量比较少。那么我们按权重再来采一次，就可以增加采到重要性样本的概率，成功的弥补了重要性采样带来的缺陷，有效的弥补采样不均衡的问题。

2 Markov Chain

在上一小节中，我们描述了三种采样方法，

概率分布采样
拒绝采样法
重要性采样法。

这三种采样方法在高维情况下的采样效率很低，所以我们需要另找方法。

2.1 基础概念介绍

首先我们要明确什么是 Random Process，也就是它研究的变量是一个随机变量的序列 {xtx_txt}。通
俗的说就是，随机过程就是一个序列，而这个序列中的每一个元素都是一个随机变量。

而 Markov Chain 就是一个特殊的随机过程，它的时间和状态都是离散的。并且，Markov Chain 需要满足 Markov 性质，也就是未来和过去是无关的。我们用数学的语言表达就是：
P(xt+1=x∣x1,x2,⋯,xt)=P(xt+1∣x1,x2,⋯,xt−m)P\left(x_{t+1}=x | x_{1}, x_{2}, \cdots, x_{t}\right)=P\left(x_{t+1} | x_{1}, x_{2}, \cdots, x_{t-m}\right)P(xt+1=x∣x1,x2,⋯,xt)=P(xt+1∣x1,x2,⋯,xt−m)上述公式就是一个 m 阶马尔可夫性质。

齐次 (一阶) 马尔可夫假设：第 t+1t+1t+1 时刻 xxx 的状态仅与第 ttt 时刻 xxx 的状态有关，即
P(xt+1=x∣x1,x2,⋯,xt)=P(xt+1∣xt)P\left(x_{t+1}=x | x_{1}, x_{2}, \cdots, x_{t}\right)=P\left(x_{t+1} | x_{t}\right)P(xt+1=x∣x1,x2,⋯,xt)=P(xt+1∣xt)而 P(xt+1∣xt)P\left(x_{t+1} | x_{t}\right)P(xt+1∣xt) 这个概率我们用什么来表达呢？我们定义为一个转移矩阵 [Pij][P_{i j}][Pij] , 而 Pij=P(xt+1=j∣xt=i)P_{i j}=P(x_{t+1}=j | x_{t}=i)Pij=P(xt+1=j∣xt=i)。

马氏链定理

我们现在换一个初始概率分布试一试，现在我们用[0.7,0.1,0.2]作为初始概率分布，然后这个状态作为序列概率分布的初始状态t0，将其带入这个状态转移矩阵计算 t1,t2,t3…的状态。代码如下：

有上边的例⼦可以看到，⽆论初始状态是什么，只要转移状态矩阵确定了之后，在经过无数次的转移之后，每个点的概率都会保持不变，⽽最终的这个分布概率就是平稳分布。
|
如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

2.2 平稳分布(Stationary Distribution)

一个Markov Chain 可以用下图来表示：

此图就是一个时间序列， xix_ixi 就表示在第 iii 时刻的状态，而每一个状态都是一个随机变量。而 πi\pi_iπi 描述的就是第 iii 个随机变量的分布。对于一个马氏链来讲，它在第 t+1t+1t+1 时刻的概率分布，可以被我们表达为：
πt+1(x∗)=∫πt(x)⋅P(x↦x∗)dx\pi_{t+1}\left(x^{*}\right)=\int \pi_{t}(x) \cdot P\left(x \mapsto x^{*}\right) d xπt+1(x∗)=∫πt(x)⋅P(x↦x∗)dx熟悉强化学习的同学就会觉得这个公式非常的熟悉。通俗的讲，他实际上就是所有可能转移到状态 xt+1x_{t+1}xt+1 的概率分布的和。那么什么是随机分布呢?

假如这里存在一个 π\piπ , 这里的 π\piπ 和前面的 πt\pi_{t}πt 和 πt+1\pi_{t+1}πt+1 都没有一毛钱关系。假如 π\piπ 是一个概率分布，那么它可以被我们写成一个无限维向量的形式：
π=[π(1),π(2),⋯,π(t),⋯],∑i=1∞π(i)=1\pi=[\pi(1), \pi(2), \cdots, \pi(t), \cdots], \quad \sum_{i=1}^{\infty} \pi(i)=1π=[π(1),π(2),⋯,π(t),⋯],i=1∑∞π(i)=1如果存在上面的公式使得使得下式成立：
π(x∗)=∫π(x)⋅P(x↦x∗)dx\pi\left(x^{*}\right)=\int \pi(x) \cdot P\left(x \mapsto x^{*}\right) d xπ(x∗)=∫π(x)⋅P(x↦x∗)dx我们就称 π(k)k=1∞{\pi(k)}_{k=1}^{\infty}π(k)k=1∞ 是马氏链 {xk}\left\{x_{k}\right\}{xk} 的平稳分布。看了数学的描述我相信大部分同学还是不懂这个平稳分布时是个什么东西？用通俗的话讲就是，对于一个马氏链来说，每个时刻都符合一个概率分布，如果每一个时刻的概率的分布都是一样的都是 π(k)\pi(k)π(k) , 那么我们就可以称这个马氏链满足一个平稳分布。

那么下一个问题就是我们为什么要引入平稳分布呢？其实，我们想要去求的这个 P(Z)P(Z)P(Z)，可以被我们看成是一个平稳分布 π(k)\pi(k)π(k) , 那么我们就可以通过构建一系列的 {x1,x2,⋯,xt,⋯}\left\{x_{1}, x_{2}, \cdots, x_{t}, \cdots\right\}{x1,x2,⋯,xt,⋯} 的马氏链，让它来逼近这个平稳分布。那么我们构建这样的一个马氏链，包括随机变量和转移矩阵，如果它满足平稳分布的条件，确实是可以收剑到平稳分布的。那么，就可以让构建出来的这个马氏链收敘到平稳分布来求得 P(Z)P(Z)P(Z)

题然，已经知道了什么是平稳分布了，那么下一个问题就是，我们需要知道什么样的分布可以称为平稳分布，也就是我们怎样才能构建出一个马氏链让它收敘到一个平稳分布。这里我们需要引入一个条件，也就是 Detailed Balance（平衡分布 / 细致平稳条件）：
π(x)⋅P(x↦x∗)=π(x∗)⋅P(x∗↦x)\pi(x) \cdot P\left(x \mapsto x^{*}\right)=\pi\left(x^{*}\right) \cdot P\left(x^{*} \mapsto x\right)π(x)⋅P(x↦x∗)=π(x∗)⋅P(x∗↦x)大家直观的来想想这个公式，为什么满足它就满足了是一个平稳分布呢？其实并不难想到，对于任意两个状态之间，使用概率分布称为转移概率得到的结果都是可逆的，那么这两个状态之间的分布一定是一样的。而说如果一个分布，满足 Detailed Balance 那么它一定可以是一个平稳分布，但是反过来并不能成立。而证明过程也不难，如下所示：
∫π(x)⋅P(x↦x∗)dx=∫π(x∗)⋅P(x∗↦x)dx=π(x∗)∫P(x∗↦x)⏟∑j=1∞Pij=1dx=π(x∗)\begin{aligned} \int \pi(x) \cdot P\left(x \mapsto x^{*}\right) d x &=\int \pi\left(x^{*}\right) \cdot P\left(x^{*} \mapsto x\right) d x \\ &=\pi\left(x^{*}\right) \underbrace{\int P\left(x^{*} \mapsto x\right)}_{\sum_{j=1}^{\infty} P_{i j}=1} d x \\ &=\pi\left(x^{*}\right) \end{aligned}∫π(x)⋅P(x↦x∗)dx=∫π(x∗)⋅P(x∗↦x)dx=π(x∗)∑j=1∞Pij=1∫P(x∗↦x)dx=π(x∗)也就是 πP=π\pi P =\piπP=π，乘上转移矩阵没有变化。

这样的话，我们就可以不用从定义上来证明一个随机过程是马尔可夫链，直接看它满不满足 Detailed Balance 就可以了。而且这个公式中， π\piπ 是平稳分布， P(x↦x∗)P\left(x \mapsto x^{*}\right)P(x↦x∗) 是马尔可夫链的状态转移概率, 这样就成功的将平稳分布和马尔可夫链结合在了一起。

综上，MCMC采样的思想：对于我们希望采样的分布 P(x)P(x)P(x)，如果我们能找到一个转移矩阵 p∗p^{*}p∗，使得 p∗p^{*}p∗ 作为转移矩阵的马氏链的最终的平稳分布为 P(x)P(x)P(x) ，那么我们就可以通过对马氏链进行采样得到一组样本 X∼P(x)X \sim P(x)X∼P(x) 。由于直接构造转移矩阵 p∗p^{*}p∗ 是困难的，所以通过构造满足细致平稳条件的转移矩阵来进行采样。

但当P(X)是未知的怎么处理？

3 MH采样（Metropolis Hastings Sampling）

如何结合马尔科夫链与门特卡罗方法？

上一节中我们讲解了 Detailed Balance，这是平稳分布的充分必要条件。Detailed Balance 为：
π(x)P(x↦x∗)=π(x∗)P(x∗↦x)\pi(x) P\left(x \mapsto x^{*}\right)=\pi\left(x^{*}\right) P\left(x^{*} \mapsto x\right)π(x)P(x↦x∗)=π(x∗)P(x∗↦x)这里的 P(x↦x∗)P\left(x \mapsto x^{*}\right)P(x↦x∗) 实际上就是条件概率 P(z∗∣x∗)P\left(z^{*} | x^{*}\right)P(z∗∣x∗) , 这样写只是便于理解。首先，我们需要明确一点，我们想要求的是后验概率分布 P(Z)，也就是我们推断问题的核心目的。我们求 P(Z) 主要是为了求在 P(Z) 概率分布下的一个相关函数的期望，也就是：
EP(Z)[f(Z)]≈1N∑i=1Nf(z(i))\mathbb{E}_{P(Z)}[f(Z)] \approx \frac{1}{N} \sum_{i=1}^{N} f\left(z^{(i)}\right)EP(Z)[f(Z)]≈N1i=1∑Nf(z(i))而我们是通过采样来得到 P(Z)∼{z(1),z(2),⋯,z(N)}P(Z)\sim\left\{z^{(1)}, z^{(2)}, \cdots, z^{(N)}\right\}P(Z)∼{z(1),z(2),⋯,z(N)} 样本点。而直接在 P(Z) 上采样是很困难的，因此需要采用 MCMC 的方法去采样。

π(x)\pi(x)π(x) 是最终的平稳分布，可以看成我们这里的 P(Z)，下面的问题就是求出概率转移矩阵 PijP_{i j}Pij , 才能满足 Detailed Balance 条件。知道了上面的条件以后，我们每次这样进行采样， x1∼P(x∣x1),x2∼P(x∣x1),x3∼P(x∣x2),⋯,xNx_{1} \sim P\left(x | x_{1}\right), x_{2} \sim P\left(x | x_{1}\right), x_{3} \sim P\left(x | x_{2}\right), \cdots, x_{N}x1∼P(x∣x1),x2∼P(x∣x1),x3∼P(x∣x2),⋯,xN，最终就可以得到我们想要的 N 个样本。

实际上真正操作的时候，要周期性判断是否收敛（进入平稳分布），需要删除掉平稳分布前采的样本。

3.1 Proposal Matrix

接下来就需要解决问题：如何找到一个满足Detailed Balance的转移矩阵？

那我们怎么来找这个状态转移矩阵 PijP_{i j}Pij 呢？首先我们可以随机一个状态转移矩阵 QijQ_{i j}Qij ，也就是 Proposal Matrix 那么肯定是：
P(Z)Q(Z∗→Z)≠P(Z∗)Q(Z→Z∗)P(Z) Q\left(Z^{*}\to Z\right) \neq P\left(Z^{*}\right) Q\left(Z \to Z^{*}\right)P(Z)Q(Z∗→Z)=P(Z∗)Q(Z→Z∗)那么我们就要想办法找到 QijQ_{i j}Qij 使得：
P(Z)Q(Z∗→Z)=P(Z∗)Q(Z→Z∗)P(Z) Q\left(Z^{*} \to Z\right)=P\left(Z^{*}\right) Q\left(Z \to Z^{*}\right)P(Z)Q(Z∗→Z)=P(Z∗)Q(Z→Z∗)注意， Q(Z→Z∗)Q\left(Z \to Z^{*}\right)Q(Z→Z∗) 实际上就是条件概率 Q(Z∣Z∗)Q\left(Z | Z^{*}\right)Q(Z∣Z∗)，二者是一样的。

那么，我们怎么来解决这个问题呢？我们可以在左右两边乘上一个因子来解决这个问题。也就是,
P(Z)Q(Z∗→Z)α(Z∗,Z)⏟P(Z↦Z∗)=P(Z∗)Q(Z→Z∗)α(Z,Z∗)⏟P(Z∗↦Z)P(Z) \underbrace{Q\left(Z^{*} \to Z\right) \alpha\left(Z^{*}, Z\right)}_{P\left(Z \mapsto Z^{*}\right)}=P\left(Z^{*}\right) \underbrace{Q\left(Z \to Z^{*}\right) \alpha\left(Z, Z^{*}\right)}_{P\left(Z^{*} \mapsto Z\right)}P(Z)P(Z↦Z∗)Q(Z∗→Z)α(Z∗,Z)=P(Z∗)P(Z∗↦Z)Q(Z→Z∗)α(Z,Z∗)而 α(z,z∗)\alpha\left(z, z^{*}\right)α(z,z∗) 定义为接收率，大小为：
α(z,z∗)=min⁡(1,P(Z∗)Q(Z→Z∗)P(Z)Q(Z∗→Z))\alpha\left(z, z^{*}\right)=\min \left(1, \frac{P\left(Z^{*}\right) Q\left(Z\to Z^{*}\right)}{P(Z) Q\left(Z^{*}\to Z\right)}\right)α(z,z∗)=min(1,P(Z)Q(Z∗→Z)P(Z∗)Q(Z→Z∗))这样定义就行了？就可以满足 Detailed Balance 吗？我们可以证明一下,
P(Z)Q(Z∗→Z)α(Z,Z∗)=P(Z)Q(Z∗→Z)min⁡(1,P(Z∗)Q(Z→Z∗)P(Z)Q(Z∗→Z))=min⁡(P(Z)Q(Z∗→Z),P(Z∗)Q(Z→Z∗))=P(Z∗)Q(Z→Z∗)min⁡(P(Z)Q(Z∗→Z)P(Z∗)Q(Z→Z∗),1)=P(Z∗)Q(Z→Z∗)α(Z∗,Z)\begin{aligned} P(Z) Q\left(Z^{*} \to Z\right) \alpha\left(Z, Z^{*}\right) &=P(Z) Q\left(Z^{*} \to Z\right) \min \left(1, \frac{P\left(Z^{*}\right) Q\left(Z \to Z^{*}\right)}{P(Z) Q\left(Z^{*} \to Z\right)}\right) \\ &=\min \left(P(Z) Q\left(Z^{*} \to Z\right), P\left(Z^{*}\right) Q\left(Z\to Z^{*}\right)\right) \\ &=P\left(Z^{*}\right) Q\left(Z \to Z^{*}\right) \min \left(\frac{P(Z) Q\left(Z^{*}\to Z\right)}{P\left(Z^{*}\right) Q\left(Z \to Z^{*}\right)}, 1\right) \\ &=P\left(Z^{*}\right) Q\left(Z\to Z^{*}\right) \alpha\left(Z^{*}, Z\right) \end{aligned}P(Z)Q(Z∗→Z)α(Z,Z∗)=P(Z)Q(Z∗→Z)min(1,P(Z)Q(Z∗→Z)P(Z∗)Q(Z→Z∗))=min(P(Z)Q(Z∗→Z),P(Z∗)Q(Z→Z∗))=P(Z∗)Q(Z→Z∗)min(P(Z∗)Q(Z→Z∗)P(Z)Q(Z∗→Z),1)=P(Z∗)Q(Z→Z∗)α(Z∗,Z)那么我们就成功的证明了：
Q(Z∗→Z)α(Z∗,Z)⏟P(Z↦Z∗)=P(Z∗)Q(Z→Z∗)α(Z,Z∗)⏟P(Z∗↦Z)\underbrace{Q\left(Z^{*}\to Z\right) \alpha\left(Z^{*}, Z\right)}_{P\left(Z \mapsto Z^{*}\right)}=P\left(Z^{*}\right) \underbrace{Q\left(Z \to Z^{*}\right) \alpha\left(Z, Z^{*}\right)}_{P\left(Z^{*} \mapsto Z\right)}P(Z↦Z∗)Q(Z∗→Z)α(Z∗,Z)=P(Z∗)P(Z∗↦Z)Q(Z→Z∗)α(Z,Z∗)所以，P(Z) 在转移矩阵 Q(Z∣Z∗)α(Z∗,Z)Q\left(Z | Z^{*}\right) \alpha\left(Z^{*}, Z\right)Q(Z∣Z∗)α(Z∗,Z) 下是一个平稳分布，也就是一个马尔可夫链，通过在这个马尔可夫链中采样就可以得到我们的相应的数据样本点了。实际上这就是大名鼎鼎的 Metropolis-Hastings 采样法。

3.2 Metropolis-Hastings Sampling

第一步，我们从一个均匀分布中进行采样， u∼U(0,1)u \sim U(0,1)u∼U(0,1)
第二步，从 Q(Z∣Z(i−1))Q\left(Z | Z^{(i-1)}\right)Q(Z∣Z(i−1)) 中进行采样得到样本点 Z∗Z^{*}Z∗
第三步，计算接受率， α=min⁡(1,P(Z∗)Q(Z∣Z∗)P(Z)Q(Z∗∣Z))\alpha=\min \left(1, \frac{P\left(Z^{*}\right) Q\left(Z | Z^{*}\right)}{P(Z) Q\left(Z^{*} | Z\right)}\right)α=min(1,P(Z)Q(Z∗∣Z)P(Z∗)Q(Z∣Z∗))。注意, 这里的 P(Z)=P^(Z)ZpP(Z)=\frac{\hat{P}(Z)}{Z_{p}}P(Z)=ZpP^(Z) 。其中 ZpZ_{p}Zp 指的是归一化因子，几乎非常难以计算，所以一般是未知的。而 P^(Z)=likelihood×prior\hat{P}(Z)=likelihood \times priorP^(Z)=likelihood×prior。所以这里的 P(Z)P(Z)P(Z) 和 P(Z∗)P\left(Z^{*}\right)P(Z∗) 就是 P^(Z)\hat{P}(Z)P^(Z) 和 P^(Z∗)\hat{P}\left(Z^{*}\right)P^(Z∗) 。由于归一化因子被抵消了，干脆就直接写成了 P(Z)P(Z)P(Z) 和 P(Z∗)P\left(Z^{*}\right)P(Z∗)
第四步，如果 u≤αu \leq \alphau≤α 时 Zi=Z∗Z^{i}=Z^{*}Zi=Z∗, 不然 Zi=Z(i−1)Z^{i}=Z^{(i-1)}Zi=Z(i−1) 这样执行了 N 次，就可以得到 {Z(1),Z(2),⋯,Z(N)}\left\{Z^{(1)}, Z^{(2)}, \cdots, Z^{(N)}\right\}{Z(1),Z(2),⋯,Z(N)} 个样本点。

3.3 小结

4 Gibbs Sampling（吉布斯采样）

MH采样是使用一个状态转移矩阵构造出一个接受率，从转移矩阵对应的条件概率中采样，然后根据接受率决定要不要接受，若接受就让它作为一个新的样本，若不接受就将这一时刻的样本取自上一时刻样本（直接复制）

MH采样有两个缺点：

⼀是需要计算接收率。尤其是在特征维度是⾼维空间时，计算量较⼤，并且由于接受率需要计算时间导致收敛时间变⻓。
⼆是有些⾼维数据特征的条件概率分布可能会⽐较好求，但是特征联合分布不好求，

综上两点，因此需要⼀个好的⽅法的改进MH采样。Gibbs 采样⽤于解决高维情况下的采样问题。

如果我们要向一个高维的分布 P(Z)=P(Z1,Z2,⋯,ZN)P(Z)=P\left(Z_{1}, Z_{2}, \cdots, Z_{N}\right)P(Z)=P(Z1,Z2,⋯,ZN) 中进行采样。那么我们怎么来进行采样呢？我们的思想就是一维一维的来，在对每一维进行采样的时候固定住其他的维度，这就是 GibbsSampling 。

我们首先规定一个 z−iz_{-i}z−i 是去除 ziz_{i}zi 后的序列， {z1,z2,⋯,zi−1,zi+1,⋯,zN}\left\{z_{1}, z_{2}, \cdots, z_{i-1}, z_{i+1}, \cdots, z_{N}\right\}{z1,z2,⋯,zi−1,zi+1,⋯,zN}

4.1 A Example

假设 ttt 时刻，我们获得的样本为 z1(t),z2(t),z3(t)z_{1}^{(t)}, z_{2}^{(t)}, z_{3}^{(t)}z1(t),z2(t),z3(t) 那么 t+1t+1t+1 时刻，我们的采样顺序为：
z1(t+1)∼P(z1∣z2(t),z3(t))z2(t+1)∼P(z2∣z1(t+1),z3(t))z3(t+1)∼P(z3∣z1(t+1),z2(t+1))\begin{array}{l} z_{1}^{(t+1)} \sim P\left(z_{1} | z_{2}^{(t)}, z_{3}^{(t)}\right) \\ \\ z_{2}^{(t+1)} \sim P\left(z_{2} | z_{1}^{(t+1)}, z_{3}^{(t)}\right) \\ \\ z_{3}^{(t+1)} \sim P\left(z_{3} | z_{1}^{(t+1)}, z_{2}^{(t+1)}\right) \end{array}z1(t+1)∼P(z1∣z2(t),z3(t))z2(t+1)∼P(z2∣z1(t+1),z3(t))z3(t+1)∼P(z3∣z1(t+1),z2(t+1))其中每次用到了前面最新得到的数据，如下图中蓝色虚线：

从这个例子中，我们应该可以大致理解固定其他的维度然后进行一维一维采样的意思了。而实际上 Gibbs 是一种特殊的 MH 采样，为什么呢？我们来证明一下。

4.2 接受率 α\alphaα 的计算

假如我们需要采样的分布 P(Z)P(Z)P(Z) 中 ZZZ 是一个多维变量，即 Z=(Z1,Z2,⋯,Zp)Z=\left(Z_{1}, Z_{2}, \cdots, Z_{p}\right)Z=(Z1,Z2,⋯,Zp)
我们使 ZZZ 发生转移的时候仅转移其中一个坐标 ZiZ_iZi，令转移后的坐标为 Zi∗Z_i^*Zi∗，即
Z⇔Z∗Z \Leftrightarrow Z^*Z⇔Z∗(Z1,Z2,⋯,Zi,⋯,Zp)⇔(Z1,Z2,⋯,Zi∗,⋯,Zp)\left(Z_{1}, Z_{2}, \cdots, Z_i,\cdots,Z_{p}\right) \Leftrightarrow \left(Z_{1}, Z_{2}, \cdots, Z_i^*,\cdots,Z_{p}\right)(Z1,Z2,⋯,Zi,⋯,Zp)⇔(Z1,Z2,⋯,Zi∗,⋯,Zp)Z→Zi∗Z\to Z_i^*Z→Zi∗ 的转移概率是 P(Zi∗∣Z1,Z2,⋯,Zi−1,Zi+1,⋯,Zp)P\left(Z_i^*|Z_{1}, Z_{2}, \cdots, Z_{i-1},Z_{i+1},\cdots,Z_{p}\right)P(Zi∗∣Z1,Z2,⋯,Zi−1,Zi+1,⋯,Zp)
Zi∗→ZZ_i^*\to ZZi∗→Z 的转移概率是 P(Zi∣Z1,Z2,⋯,Zi−1,Zi+1,⋯,Zp)P\left(Z_i|Z_{1}, Z_{2}, \cdots, Z_{i-1},Z_{i+1},\cdots,Z_{p}\right)P(Zi∣Z1,Z2,⋯,Zi−1,Zi+1,⋯,Zp)

则关于 P(Z)P(Z)P(Z)的细致平稳条件方程是：
P(Zi∗∣Z1,Z2,…Zi−1Zi+1,…Zp)α(Z,Z∗)=P(Z∗)P(Zi∣Z1∗,Z2∗,…,Zi−1∗Zi+1∗,…Zp∗)P\left(Z_{i}^{*} | Z_{1}, Z_{2}, \ldots Z_{i-1} Z_{i+1}, \ldots Z_{p}\right) \alpha\left(Z, Z^{*}\right)=P\left(Z^{*}\right) P\left(Z_{i} | Z_{1}^{*}, Z_{2}^{*}, \ldots,Z_{i-1}^{*} Z_{i+1}^{*}, \ldots Z_{p}^{*}\right)P(Zi∗∣Z1,Z2,…Zi−1Zi+1,…Zp)α(Z,Z∗)=P(Z∗)P(Zi∣Z1∗,Z2∗,…,Zi−1∗Zi+1∗,…Zp∗)

令 P(Zi∗∣Z1,Z2,…Zi+1,Zi+1,…Zp)=Q(Z,Z∗)P\left(Z_{i}^{*} | Z_{1}, Z_{2}, \ldots Z_{i+1}, Z_{i+1,} \ldots Z_{p}\right)=Q\left(Z, Z^{*}\right)P(Zi∗∣Z1,Z2,…Zi+1,Zi+1,…Zp)=Q(Z,Z∗)

为精简记号，将 (Zi∗∣Z1,Z2,…Zi+1,Zi+1,…Zp)\left(Z_{i}^{*} | Z_{1}, Z_{2}, \ldots Z_{i+1}, Z_{i+1,} \ldots Z_{p}\right)(Zi∗∣Z1,Z2,…Zi+1,Zi+1,…Zp) 表示为 Z−iZ_{-i}Z−i

Q(Z,Z∗)=P(Zi∗∣Z−i)Q\left(Z, Z^{*}\right)=P\left(Z_{i}^{*} | Z_{-i}\right)Q(Z,Z∗)=P(Zi∗∣Z−i)

首先回顾一下，MH 采样的方法。我们的目的是从 Q(Z∣Z(t))Q\left(Z | Z^{(t)}\right)Q(Z∣Z(t)) 中采样获得 Z∗Z^{*}Z∗, 然后计算接受率
α=min⁡(1,P(Z∗)Q(Z∣Z∗)P(Z)Q(Z∗∣Z))\alpha=\min \left(1, \frac{P\left(Z^{*}\right) Q\left(Z | Z^{*}\right)}{P(Z) Q\left(Z^{*} | Z\right)}\right)α=min(1,P(Z)Q(Z∗∣Z)P(Z∗)Q(Z∣Z∗))首先我们来看 Q(Z∣Z(t))Q\left(Z | Z^{(t)}\right)Q(Z∣Z(t))
Q(Z∣Z(t))=Q(Zi,Z−i∣Zi(t),Z−i(t))Q\left(Z | Z^{(t)}\right)=Q\left(Z_{i}, Z_{-i} | Z_{i}^{(t)}, Z_{-i}^{(t)}\right)Q(Z∣Z(t))=Q(Zi,Z−i∣Zi(t),Z−i(t))假设我们现在是在对第 iii 维进行采样, 我们只要关注 P(Zi∗∣Z−i)P\left(Z_{i}^{*} | Z_{-i}\right)P(Zi∗∣Z−i) 所以, 我们可以得到 : Q(Z∣Z(t))=P(Zi∗∣Z−i(t))Q\left(Z | Z^{(t)}\right)= P\left(Z_{i}^{*} | Z_{-i}^{(t)}\right)Q(Z∣Z(t))=P(Zi∗∣Z−i(t))已经成功的将 Q(Z∣Z(t))Q\left(Z | Z^{(t)}\right)Q(Z∣Z(t)) 做了等价转换以后。那么我们想要求的 α\alphaα 可以被我们成功的转换成如下的形式：
α=min⁡(1,P(Zi∗∣Z−i∗)P(Z−i∗)P(Zi∣Z−i∗)P(Zi∣Z−i)P(Z−i)P(Zi∗∣Z−i))\alpha=\min \left(1, \frac{P\left(Z_{i}^{*} | Z_{-i}^{*}\right) P\left(Z_{-i}^{*}\right) P\left(Z_{i} | Z_{-i}^{*}\right)}{P\left(Z_{i} | Z_{-i}\right) P\left(Z_{-i}\right) P\left(Z_{i}^{*} | Z_{-i}\right)}\right)α=min(1,P(Zi∣Z−i)P(Z−i)P(Zi∗∣Z−i)P(Zi∗∣Z−i∗)P(Z−i∗)P(Zi∣Z−i∗))

计算到了这里，我们还是不好进行计算，上面和下面好像还是不好消除。如果我们可以得到 Z−iZ_{-i}Z−i 和 ZiZ_iZi 之间的关系就好了。下面我们会得出一个重要的结论来帮助我们计算 α\alphaα 的具体值。首先我们来举一个例子：

那么假设当 t=1 的时刻，有一个样本为： Z1(1),Z2(1),Z3(1)Z_{1}^{(1)}, Z_{2}^{(1)}, Z_{3}^{(1)}Z1(1),Z2(1),Z3(1)当 t=2t=2t=2 的时刻，我们假设先对第一维进行采样就可以得到： Z1(2),Z2(1),Z3(1)Z_{1}^{(2)}, Z_{2}^{(1)}, Z_{3}^{(1)}Z1(2),Z2(1),Z3(1) 很显然 Z2(1),Z3(1)Z_{2}^{(1)}, Z_{3}^{(1)}Z2(1),Z3(1) 根本没有发生变化。我们可以得到 Z−1=Z−1∘∗Z_{-1}=Z_{-1^{\circ}}^{*}Z−1=Z−1∘∗ 也就是在 Gibbs 采样时，采样前后只关注于一个维度，其他的维度我们都没有关注到。所以就可以得到结论：
Z−i=Z−i∗Z_{-i}=Z_{-i}^{*}Z−i=Z−i∗那么，我们根据这个结论就可以得到：
α=min⁡(1,P(Zi∗∣Z−i∗)P(Z−i∗)P(Zi∣Z−i∗)P(Zi∣Z−i∗)P(Z−i∗)P(Zi∗∣Z−i∗))=1\alpha=\min \left(1, \frac{P\left(Z_{i}^{*} | Z_{-i}^{*}\right) P\left(Z_{-i}^{*}\right) P\left(Z_{i} | Z_{-i}^{*}\right)}{P\left(Z_{i} | Z_{-i}^{*}\right) P\left(Z_{-i}^{*}\right) P\left(Z_{i}^{*} | Z_{-i}^{*}\right)}\right)=1α=min(1,P(Zi∣Z−i∗)P(Z−i∗)P(Zi∗∣Z−i∗)P(Zi∗∣Z−i∗)P(Z−i∗)P(Zi∣Z−i∗))=1那么计算出接受率为 1，也就是每次都必定被接受。所以，每次从 Q(Z∣Z(t))=P(Zi∗∣Z−i)Q\left(Z | Z^{(t)}\right)=P\left(Z_{i}^{*} | Z_{-i}\right)Q(Z∣Z(t))=P(Zi∗∣Z−i) 中进行采样得到 Zi∗Z_{i}^{*}Zi∗ 即可，一维一维的进行采样就可以采到整个高维的分布，各个维度上的样本。
所以，解释到了这里，大家基本就可以知道 Gibbs Samplings 是 α=1\alpha=1α=1 的 MH Sampling 的意义了。在 Gibbs Sampling 中 α=1\alpha=1α=1, 而且状态转移矩阵 Q(Z∣Z(t))=P(Zi∗∣Z−i(t))Q\left(Z | Z^{(t)}\right)=P\left(Z_{i}^{*} | Z_{-i}^{(t)}\right)Q(Z∣Z(t))=P(Zi∗∣Z−i(t)) , 所以 Gibbs Sampling 就是把目标分布 P 对应的条件概率当作状态转移分布 QQQ 。

这里我们需要额外提醒一下，使用 Gibbs Sampling 是有使用前提的，也就是固定其他维度后的一维分布时方便进行采样的，如果固定其他维度的时候得到的一维分布仍然是非常难进行采样的，那么使用 Gibbs Sampling 也是没有用的。

5 采样

在前面的章节中，我们已经基本介绍了Markov Chain Monte Carlo Sampling 的基本概念，基本思路和主要方法。那么这一小节中，我们将主要来介绍一下，什么是采样？我们为什么而采样？什么样的样本是好的样本？以及我们采样中主要会遇到哪些困难？

5.1 采样的动机

这一小节的目的就是我们要知道什么是采样的动机，我们为什么而采样？

首先第一点很简单，采样本身就是发出常见的任务，我们机器学习中经常需要进行采样来完成各种各样的任务。如果从一个 P(X)P(X)P(X) 中采出一堆样本。
求和求积分。包括大名鼎鼎的Monte Carlo 算法。我们求 P(X)P(X)P(X) 主要是为了求在 P(X)P(X)P(X) 概率分布下的一个相关函数的期望，也就是：
∫P(x)f(x)dx=EP(X)[f(X)]≈1N∑i=1Nf(x(i))\int P(x) f(x) d x=\mathbb{E}_{P(X)}[f(X)] \approx \frac{1}{N} \sum_{i=1}^{N} f\left(x^{(i)}\right)∫P(x)f(x)dx=EP(X)[f(X)]≈N1i=1∑Nf(x(i))而我们是通过采样来得到 P(X)∼{x(1),x(2),⋯,x(N)}P(X) \sim\left\{x^{(1)}, x^{(2)}, \cdots, x^{(N)}\right\}P(X)∼{x(1),x(2),⋯,x(N)} 样本点。

5.2 什么样的是好样本

既然，我们知道了采样的目的和动机，下一个问题就自然是，同样是采样，什么样的样本就是好样本呢？或者说是采样的效率更高一些。

首先样本的分布肯定是要趋向于原始的目标分布吧，也就是说样本要趋向于高概率选择区域。或者是说，采出来的样本出现的概率和实际的目标分布的概率保持一致。
样本和样本之间是相互独立的。这个就没有那么直观了。大家想一想就知道了，如果我采出来的一堆样本之间都差不多，那么就算采出来了趋向于高概率选择区域的样本，那采样效率太低了，样本中反映的信息量太有限了。

5.3 实际采样中的困难

实际采样中，采样是困难的，为什么呢？我们这里主要介绍两点：

Partation function is intractable. 我们的后验分布往往被写成 P(X)=1ZP^(X)P(X) = \frac 1Z \hat{P}(X)P(X)=Z1P^(X)，上面这个 P^(X)\hat{P}(X)P^(X) 都比较好求，就是等于Likelihood∗PriorLikelihood * PriorLikelihood∗Prior。而 ZZZ 就是我们要求的归一化常数，它非常的难以计算， Z=∫P^(X)dXZ=\int \hat{P}(X) d XZ=∫P^(X)dX , 这几乎就是不可计算的。所以，有很多采样方法就是想要跳过求 P(X)P(X)P(X) 的过程，来从一个近似的分布中进行采样，当然这个近似的分布采样要比原分布简单。比如：Rejection Sampling 和 Importance Sampling。
The curse of high dimension. 如果样本空间 X∈Rp\mathcal{X} \in \mathbb{R}^{p}X∈Rp , 每个维度都有 KKK 个状态的话。那么总的样本空间就有 KpK^pKp 的状态。要知道那个状态的概率高，就必须要 bbb 遍历整个样本空间，不然就不知道哪个样本的概率高，如果状态的数量是这样指数型增长的话，全看一遍之后进行采样时不可能的。所以，直接采样的方法是不可行的。

因为采样是困难的，因此之后将讲解一些采样方法：
Rejection Sampling
Importance Sampling
MCMC（ MH、Gibbs）

5.4 采样方法

Rejection Sampling 和 Importance Sampling，都是借助一个 Q(x)Q(x)Q(x) 去逼近目标分布 P(x)P(x)P(x) , 通过从 Q(x)Q(x)Q(x) 中进行采样来达到在 P(x)P(x)P(x) 中采样的目的，而且在 Q(x)Q(x)Q(x) 中采样比较简单。当时如果 Q(x)Q(x)Q(x) 和 P(x)P(x)P(x) 直接的差距太大的话，采样效率会变得很低。而 MCMC 方法，我们主要介绍了 MH Sampling 和 Gibbs Sampling，我们主要是通过构建一个马氏链去逼近目标分布，具体的描述将在下一节中展开描述。

6 Method of MCMC

这一小节主要是对前面的补充，希望可以详细的介绍一下MCMC 原理，将前面的知识点可以顺利的串起来。MCMC 采样中，我们借助了一条马氏链，马氏链的性质，经过若干步以后会收敛到一个平稳分布。马尔可夫链的组成可以大致分成两个部分：

状态空间： {1,2,3,⋯,k}\{1,2,3, \cdots, k\}{1,2,3,⋯,k}
状态转移空间 Q=[Qij]k×kQ=\left[Q_{i j}\right]_{k \times k}Q=[Qij]k×k

马尔可夫链的模型可以被我们表达为：

图1: 马尔可夫链模型抽象图

每一个时间点有一个状态分布，表示当前时间点位于某个状态的概率分布情况,我们表示为 q(t)(x)q^{(t)}(x)q(t)(x) 。如果是在 t=1t=1t=1 的时间节点，状态的概率分布为 q(1)(x)q^{(1)}(x)q(1)(x), 我们可以用下列表来描述：
x123⋯kq(1)(x)q11q12q13⋯q1k\begin{array}{c|ccccc} x & 1 & 2 & 3 & \cdots & k \\ \hline q^{(1)}(x) & q_{1}^{1} & q_{1}^{2} & q_{1}^{3} & \cdots & q_{1}^{k} \end{array}xq(1)(x)1q112q123q13⋯⋯kq1k我们假设在 t=mt=mt=m 时刻之后到达了平稳分布状态，那么我们就可以得到： q(m)=q(m+1)=q(m+2)q^{(m)}=q^{(m+1)}=q^{(m+2)}q(m)=q(m+1)=q(m+2) 这时的平稳分布就是我们想要的目标分布。相邻时间节点之间的状态转移矩阵为：
Q=[Q11Q12⋯Q1kQ21Q22⋯Q2k⋮⋮⋱⋮Qk1Qk2⋯Qkk]k×kQ=\left[ \begin{array}{cccc} Q_{11} & Q_{12} & \cdots & Q_{1 k} \\ Q_{21} & Q_{22} & \cdots & Q_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ Q_{k 1} & Q_{k 2} & \cdots & Q_{k k} \end{array}\right]_{k \times k}Q=⎣⎢⎢⎢⎡Q11Q21⋮Qk1Q12Q22⋮Qk2⋯⋯⋱⋯Q1kQ2k⋮Qkk⎦⎥⎥⎥⎤k×k
状态转移矩阵描述的是，Qij=Q(x2=j∣x1=i)Q_{ij}=Q\left(x^{2}=j | x^{1}=i\right)Qij=Q(x2=j∣x1=i) 。描述的是从一个状态转移到另外一个状态的概率。所以，状态转移矩阵的每一行 iii 表示为目前状态是 iii 时，到其他状态的概率，那么必然有 ∑k=1kQik=1\sum_{k=1}^{k} Q_{i k}=1∑k=1kQik=1。

6.1 Markov Chain 收敛性介绍

在这一小节中，我们将详细的介绍一下，Markov Chain 中状态转移的过程。并将证明在Markov Chain 随着迭代一定会收敛到一个平稳分布。

6.1.1 Markov Chain 状态转移计算

假设在 t+1t+1t+1 时刻，状态是 x=jx=jx=j, 那么它的分布为所有可能转移到这个状态的概率 iii 乘以这个状态的分布 q(t)(x=i)q^{(t)}(x=i)q(t)(x=i) , 我们用公式表达就是：
q(t+1)(x=j)=∑i=1kq(t)(x=i)Qijq^{(t+1)}(x=j)=\sum_{i=1}^{k} q^{(t)}(x=i) Q_{i j}q(t+1)(x=j)=i=1∑kq(t)(x=i)Qij那么，这仅仅是当 x=jx=jx=j 时概率，实际上在 t+1t+1t+1 时刻，可能出现的状态有 k\mathrm{k}k 个，那么 qt+1q^{t+1}qt+1 的分布，就是将转移到各个状态的概率分别计算出来，也就是如下所示：
q(t+1)=[q(t+1)(x=1)q(t+1)(x=2)q(t+1)(x=3)⋯q(t+1)(x=k)]1×kq^{(t+1)}=\left[q^{(t+1)}(x=1) \quad q^{(t+1)}(x=2) \quad q^{(t+1)}(x=3) \quad \cdots \quad q^{(t+1)}(x=k)\right]_{1 \times k}q(t+1)=[q(t+1)(x=1)q(t+1)(x=2)q(t+1)(x=3)⋯q(t+1)(x=k)]1×k而,
q(t+1)(x=j)=∑i=1kq(t)(x=i)Qijq^{(t+1)}(x=j)=\sum_{i=1}^{k} q^{(t)}(x=i) Q_{i j}q(t+1)(x=j)=i=1∑kq(t)(x=i)Qij 那么， q(t+1)q^{(t+1)}q(t+1) 可以被我们表示为：
q(t+1)=[∑i=1kq(t)(x=i)Qi1∑i=1kq(t)(x=i)Qi2⋯∑i=1kq(t)(x=i)Qik]1×k=q(t)⋅Q\begin{aligned} q^{(t+1)} &=\left[\sum_{i=1}^{k} q^{(t)}(x=i) Q_{i 1} \quad \sum_{i=1}^{k} q^{(t)}(x=i) Q_{i 2} \quad \cdots \quad \sum_{i=1}^{k} q^{(t)}(x=i) Q_{i k}\right]_{1 \times k} \\ &=q^{(t)} \cdot Q \end{aligned}q(t+1)=[i=1∑kq(t)(x=i)Qi1i=1∑kq(t)(x=i)Qi2⋯i=1∑kq(t)(x=i)Qik]1×k=q(t)⋅Q其中, q(t)=[q(t)(x=1)q(t)(x=2)q(t)(x=3)⋯q(t)(x=k)]1×kq^{(t)}=\left[q^{(t)}(x=1) \quad q^{(t)}(x=2) \quad q^{(t)}(x=3) \quad \cdots \quad q^{(t)}(x=k)\right]_{1 \times k}q(t)=[q(t)(x=1)q(t)(x=2)q(t)(x=3)⋯q(t)(x=k)]1×k 。那么，通过这个逆推公式，我们可以得到， q(t+1)=q(t)Q=q(t−1)Q2=⋯=q(1)Qtq^{(t+1)}=q^{(t)} Q=q^{(t-1)} Q^{2}=\cdots=q^{(1)} Q^{t}q(t+1)=q(t)Q=q(t−1)Q2=⋯=q(1)Qt 。

通过上述的描述，详细大家都已经详细的了解了 Markov Chain 中，每个时刻点的状态的分布 q(t)q^{(t)}q(t) 的计算方法。既然我们知道了每个时间点的概率分布的计算方法，下一个问题就是我们怎么可以知道一定是收敛的呢?

6.1.2 Markov Chain 收剑性

由于 Q 是一个随机概率矩阵，那么我们可以得到，每个值都是小于 1 的，所以也必然有特征值的绝对值 ≤1\leq1≤1 。为什么呢？我们可以从特征值的几何意义上好好的想一想，特征值代表变换中方向不变的向量的变化尺度。随机矩阵的变化尺度必然是小于等于 1 的。所以，我们可以对概率转移矩阵做特征值分解，分角成对角短阵：
Q=AΛA−1Λ=[λ1λ2⋱λk],∣λi∣≤1(i=1,2,⋯,k)Q=A \Lambda A^{-1} \quad \Lambda=\left[\begin{array}{cccc} \lambda_{1} & & & \\ & \lambda_{2} & & \\ & & \ddots & \\ & & & \lambda_{k} \end{array}\right], \quad\left|\lambda_{i}\right| \leq 1 \quad(i=1,2, \cdots, k)Q=AΛA−1Λ=⎣⎢⎢⎡λ1λ2⋱λk⎦⎥⎥⎤,∣λi∣≤1(i=1,2,⋯,k)随机矩阵的特征值都小于等于1，我们假设只有一个 λi=1\lambda_{i}=1λi=1, 则：
q(t+1)=q(1)(AΛA−1)t=q(1)AΛtA−1q^{(t+1)}=q^{(1)}\left(A \Lambda A^{-1}\right)^{t}=q^{(1)} A \Lambda^{t} A^{-1}q(t+1)=q(1)(AΛA−1)t=q(1)AΛtA−1当 t→∞t \rightarrow \inftyt→∞ 时，必然有：
Λt=[01⋱0]\Lambda^{t}=\left[\begin{array}{cccc} 0 \\ & 1 \\ & & \ddots \\ & & & 0 \end{array}\right]Λt=⎣⎢⎢⎡01⋱0⎦⎥⎥⎤
我们可以假设存在足够大的 M：
s.t. ΛM=[01⋱0]\text {s.t. } \quad \Lambda^{M}=\left[\begin{array}{cccc} 0 \\ & 1 \\ & & \ddots \\ & & & 0 \end{array}\right]s.t. ΛM=⎣⎢⎢⎡01⋱0⎦⎥⎥⎤
所以,
q(m+1)=q(1)AΛmA−1q(m+2)=q(m+1)AΛA−1=q(1)AΛmA−1AΛA−1=q(1)AΛ(m+1)A−1=q(m+1)\begin{aligned} q^{(m+1)} &=q^{(1)} A \Lambda^{m} A^{-1} \\ q^{(m+2)} &=q^{(m+1)} A \Lambda A^{-1} \\ &=q^{(1)} A \Lambda^{m} A^{-1} A \Lambda A^{-1} \\ &=q^{(1)} A \Lambda^{(m+1)} A^{-1} \\ &=q^{(m+1)} \end{aligned}q(m+1)q(m+2)=q(1)AΛmA−1=q(m+1)AΛA−1=q(1)AΛmA−1AΛA−1=q(1)AΛ(m+1)A−1=q(m+1)通过上述的证明，我们成功的证明了 q(m+2)=q(m+1)q^{(m+2)}=q^{(m+1)}q(m+2)=q(m+1)。我们用数学的语言来表述一下，也就是当 t>mt>mt>m 时, q(m+1)=q(m+2)=⋯=q(∞)q^{(m+1)}=q^{(m+2)}=\cdots=q^{(\infty)}q(m+1)=q(m+2)=⋯=q(∞) 。这就是平稳分布，我们成功的证明了 Markov Chain 经过足够大的步数 mmm 之后，一定会收敛到一个平稳分布。于是，这就启发了我们设计一个 Markov Chain, 收敛到我们想要采样的分布 p(x)p(x)p(x) 。那么，怎么我们才能让它收敛呢? 实际上就是由状态转移矩阵 QQQ 所决定的。我们的核心问题就是设计一个合适的状态转移矩阵 Q。

那么，我们要做的就是设计一个 MCMC，利用 Monte Chain 收敛到一个平稳分布 q(x)q(x)q(x), 使得平稳分布 ≈\approx≈ 目标分布 p(x)p(x)p(x) 。也就是当 mmm 足够大的时候, q(m)(x)=q(m+1)(x)=q(m+2)(x)=q(x)q^{(m)}(x)=q^{(m+1)}(x)=q^{(m+2)}(x)=q(x)q(m)(x)=q(m+1)(x)=q(m+2)(x)=q(x)

那么，我们的 Markov Chain 解决了当维度很高的时候， q(x)≈p(x)q(x) \approx p(x)q(x)≈p(x) 找不到的情况，在 MCMC 中不要显示的去找，而是构建一个 Markov Chain 去近似，跳过了直接去寻找的过程。

这里我们介绍一个概念，也就是从开始到收敛到 mmm 的这段时期被我们称为bum-in，中文翻译为燃烧期(个人觉得非常的难听，所以我从来不用中文的表述形式)。也有说法称这个时间 ttt 为Mix-time。当然也不是任何的分布都可以用MCMC 来进行采样。但是它可以有效的避免我们去寻找 q(z)q(z)q(z) 。下面我们将描述一些用MCMC 采样时遇到的困难的地方。

6.2 Existing Problem

虽然，我们可以证明出MCMC 最终可以收敛到一个平稳分布。但是并没有理论来判断MarkovChain 是否进入了平稳分布，也就是不知道Markov Chain 什么时候会收敛。
Mixing Time 过长，这就是有高维所造成的，维度和维度之间的相关性太强了， p(x)p ( x )p(x) 太过于复杂了。理论上 MCMCMCMCMCMC 是可以收敛的，但是如果 p(x)p(x)p(x) 太过于复杂了，我们基本就是认为它是不收敛的。实际上，现在有各种各样的MCMC 算法的变种都是在想办法解决这个Mixing Time 过长的问题。
3 我们希望采到的样本之间的样本相互独立，也就是采到的样本之间的相关性越小越好。

这个有关于样本之间独立性的问题，大家可能不太好理解，这是实际在高维分布中我们采用MCMC来进行采样很有可能造成样本单一，相关性太强的问题。我们我们来举一个Mixture Gaussian Distribution的例子。下图所示是一个Mixture Gaussian Distribution 的例子：

会有一个什么问题呢？就是样本都趋向于一个峰值附近，很有可能会过不了低谷，导致样本都聚集在一个峰值附近。这个问题出现的原因我们可以从能量的角度来解释这个问题。在无向图中，我们常用下列公式来进行表示：
P(X)=1ZP^(X)=1Zexp⁡−E(X)P(X)=\frac{1}{Z} \hat{P}(X)=\frac{1}{Z} \exp ^{-\mathbb{E}(X)}P(X)=Z1P^(X)=Z1exp−E(X)实际上这里的 −E(X)-\mathbb{E}(X)−E(X) 指的就是能量函数，能量和概率是成反比的，概率越大意味着能量越低，能量越低，越难发生跳跃的现象。所以，采样很容易陷入到一个峰值附近。并且，多峰还可以分为均匀和陡峭，陡峭的情况中，能量差实在是太大了，就很难发生跳跃。就像孙悟空翻出如来佛祖的五指山一样，佛祖的维度很高，孙悟空在翻跟头的时候，一直在一个低维里面不同的打转，根本就跳不出来，就是来自佛祖的降维打击。
所以，在高维情况下，很容易发生在一个峰值附近不停的采样，根本就跳不出来，导致采到的样本的多样性低，样本之间的关联性大，独立性低。

7 参考

MCMC随机采样.

版权声明：本文为CSDN博主「AI路上的小白」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/cengjing12/article/details/106586088

13 MCMC(Markov Chain Monte Carlo)相关推荐

Markov Chain Monte Carlo 和 Gibbs Sampling算法
Welcome To My Blog 一.蒙特卡洛模拟蒙特卡洛模拟(Monte Carlo Simulation)是随机模拟的别名,关于随机模拟的一个重要的问题就是:给定一个概率分布p(x),如何生 ...
马尔可夫链蒙特卡罗法（Markov Chain Monte Carlo，MCMC）
文章目录 1. 蒙特卡罗法 2. 马尔可夫链 3. 马尔可夫链蒙特卡罗法 4. Metropolis-Hastings 算法 5. 吉布斯抽样蒙特卡罗法(Monte Carlo method),也称 ...
论文辅助笔记（代码实现）：Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo
1 主要思路回顾具体可见:论文笔记 Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo (ICML 2 ...
Markov Chain Monte Carlo
转载至https://zhuanlan.zhihu.com/p/25610149 [数据分析] Markov Chain Monte Carlo Markov Chain Monte Carlo简称M ...
从概率论到Markov Chain Monte Carlo（MCMC）-- 转
大学本科时代开始学习的概率论,从变着花样从箱子里取不同颜色的球计算概率,到计算各种离散或连续的随机分布期望.方差,再高深点就是利用生成函数求期望和方差,再就是估计理论,包括点估计.极大 ...
【ML】Markov Chain Monte Carlo(MCMC)---Slice sampler(切片采样）和Hierarchical Models(层次模型)
导航 Slice sampler 2D slice sample General Slice Sampler Hierarchical models python Code download Refe ...
R语言与Markov Chain Monte Carlo(MCMC)方法学习笔记(2)
前面已经大致的叙述了MCMC方法.今天来分享一下R中的一个实现MCMC算法的包mcmc. mcmc包的一个核心函数就是metrop,其调用格式为: metrop(obj, initial, nbatc ...
论文笔记 Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo （ICML 2008）
0 摘要低秩矩阵逼近方法是协同过滤中最简单.最有效的方法之一.这类模型通常通过寻找模型参数的MAP估计来拟合数据,这一过程即使在非常大的数据集上也能有效地执行. 然而,除非正则化参数被仔细地调整,否 ...
martingale、markov chain、Monte Carlo、MCMC
文章结构如下: 1: MCMC 1.1 MCMC是什么 1.2 为什么需要MCMC 2: 蒙特卡罗 2.1 引入 2.2 均匀分布,Box-Muller 变换 2.3 拒绝接受采样(Acceptanc ...

13 MCMC(Markov Chain Monte Carlo)