如何科学地采样

发表于 更新于

这是你的分布（滑稽

如何只采集1k个点，就描述出这个分布的特点呢？

为什么要采样？

最直接的目的是用离散的数据 近似 代替连续（或原先就是离散）的数据，比如说要计算一个很复杂的函数f的积分，假设找出这个函数原函数F是十分困难的，那么可以用近似的值代替真实的值。一个简单而有效的方式是在定义域内按每隔0.01的距离计算f的值，然后用小矩形面积——0.01*f(x)表示这段范围内函数的积分，最后整个函数的面积用众多小矩形的面积之和代替。计算函数积分确实是采样的一个重要应用。

采样可以很大程度减轻计算难度，不过怎么采样还是有很深的学问的，采样的结果的分布要能近似真实分布。

换句话说，概率密度大的地方要能多采集点，而概率密度小的地方少采集点。

比方说上面这个例子里，同样采集1w个点，左边的采样效果就更符合真实分布。

怎么采样？

概率分布采样

假设现在有一个分布，已知该分布的概率密度函数（pdf），且可以计算出它的分布函数（cdf）。然后用计算机生成一个均匀分布$Y \sim U(0,1)$，在Y上的采样结果通过cdf映射回X，就得到了该分布的均匀采样。采样的思想用下面的图形象地表示

最后公式一下：
$$
Y \sim U(0,1) \
X = cdf^{−1}(Y)
$$

不过万一pdf太复杂以至于cdf无法计算呢？那么就无法用概率分布采样了。

拒绝采样

请一个帮手

拒绝采样通过引入一个已知pdf和cdf的分布，来迂回计算原先难以采样的分布，这个被引入的分布就称为“参考分布”，不妨记作$q(x)$，而待采样记作$p(x)$。拒绝采样还需要另一个系数$k$，并保证在需要采样的区间内，$\forall x, kq(x) \ge p(x)$。

由于$q(x)$的cdf已知，所以可以很轻易得到$q(x)$的合理采样，假设得到的样本为$X={x_1,\dots,x_n}$。但此时属于$X$的点的分布并不符合$q(x)$的分布，那么接下来只需要遍历每个点$x_i$，同时每次扔一枚能产生0到1之间均匀分布随机数的“骰子”，将“骰子”得到的数u和接受率$\alpha = q(x_i)/p(x_i)$作比较，如果满足$u \le \alpha$就接受这个样本点，否则丢弃即可。

拒接采样之所以合理，其物理含义也是很直观的，因为乘上系数$k$后，q(x)很大于等于p(x)，根据q(x)分布采集得到的样本点集$X$近似曲线q(x)下的面积，而接受率又是比较q(x)和p(x)面积的比，所以经过接受率筛选后剩下的点集就近似p(x)曲线下的面积了。

不过拒接采样还是有很多限制的，比如q(x)一定要尽可能接近p(x)，否则采样的效率会很低。其次k的选取也是要恰到好处，在满足拒绝采样等式的要求下也要尽可能小。

重要性采样

也是请一个帮手

重要性采样的用武之地是计算一个函数f关于分布p的期望，即$E(f(x)), x \sim p$。并且这个p不仅贼难采样，而且真实分布比较偏，换句话说，即使采集到了一些样本点，但是这些样本点大多对应于f取值较小的部位，就像下面这张图所示的

为了解决这个问题，还是引入一个简单下手的参考分布q，并且希望通过q改善p“偏离群众”的问题。

集体做法是这样的，上面说到要计算期望$E(f(x))$，详细写开来就是

$$E(f(x))=\int_x{f(x)p(x)dx}$$

引入q后做一下置换变换

$$E(f(x))=\int_x{f(x)\frac{p(x)}{q(x)}q(x)dx}$$

经过这么一个操作后，原来求f关于p的期望变成了求$f(x)\frac{p(x)}{q(x)}$关于q的期望了！因为q的采样容易，并且和上图中展示的一样，原本对f来说不容易采集到的却又十分重要的部分，受$q(x)$的影响变得十分容易得到了。而这样一来，所求的期望的偏差就不会因为数据稀疏偏差太远。

话说重要性采样在一些论文中出现的频率还是很大的，确实是重要的一种采样

MCMC

本文的核心角色出现了，MCMC全称markov-chain monte carlo，含义是在马尔可夫链上应用随机化的方法。马尔可夫链可以由一个转移矩阵描述，当给定一个初始状态，然后这个初始状态开始在转移矩阵的作用下随机切换状态。也就是说，MCMC将马尔可夫链视为工具，用随机化的手段模拟状态的转移。

说到马尔可夫，有一个很重要的性质是不得不说的，那就是它的“平稳分布”。假如随意给定一个初始状态，然后开始状态转移，接着对每个状态出现次数的频率做统计。这个过程可以用矩阵乘法描述，假设初始状态为$A_0 \in R^{n \times 1}$，转移矩阵为$M \in R^{n \times n}$，那么每一次转移可以看作是这样一个变换

$$
A_{t+1} = MA_t
$$

如果从初始状态开始计算的话，那么t时刻的状态可以表示为

$$
A_t = M^{t}A_0
$$

而随着t的不断不断增大，分布也会近似趋于平稳，就像下面这张图表示的那样，纵使初始分布在空间的各个角落，随着迭代的继续，最终分布会趋向于稳定。

也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。

只要马尔可夫链收敛到平稳分布π后，就存在一个被称为“细致平衡”的等式，

$$
\pi(x)P(y|x) = \pi(y)P(x|y)
$$

等式中的$P(y│x)$就是马尔可夫转移矩阵中的元素，表示从x到y的转移概率。这个等式很好地描述了马尔可夫平稳的本质。

说完了马尔可夫链模型的性质，下面要开始引入正题了，即MCMC是怎么把马尔可夫链模型和蒙特卡洛方法结合起来，并采集出和原分布一致的样本分布的。

M-H（Metropolis-Hastings）

M-H是MCMC的一种实现。回顾一下，上面这个细致平衡的等式是怎么的来的呢？第一步，也是最关键的是找到了一个满足马尔可夫性质的转移矩阵M，于是才有了之后的平稳分布，也就最终能达到细致平衡条件。

然而，在“对原分布采样”这个需求下，我们只有平稳分布（也就是原分布本身），却没有满足该平稳分布的状态转移矩阵。所以，如何构造这样的一个矩阵呢？毕竟不是所有的矩阵都能推导出想要的平稳分布。

M-H的做法十分大胆，它拿一个和目标转移矩阵一样大的随机马尔可夫转移矩阵Q作为转移矩阵的一部分。显然此时$\pi(x)Q(y|x) \neq \pi(y)Q(x|y)$。而后，为了让等式相等，引入类似拒绝采样里提到的接受率这样一个概念。令$\alpha(x,y) = min⁡(1,\frac{\pi(x)Q(y|x)}{\pi(y)Q(x|y)})$，于是一个细致平衡就被构造出来了

$$
\pi(x)Q(y|x)\alpha(x,y) = \pi(y)Q(x|y)\alpha(y,x)
$$

因为Q是随机的，也就是任意满足马尔可夫性质的转移矩阵都能work，那么简单起见，甚至可以把Q设置为对称矩阵，此时接受率可以简化为：$\alpha(x,y) = min⁡(1, \pi(x) / \pi(y))$。M-H算法的步骤如下：

  1. 输入需要采样的平稳分布$\pi$，转移的次数n；
  2. 随机构造一个马尔可夫转移矩阵Q；
  3. 挑选一个初始状态$x_0$；
  4. for t = 1 to n：
    1. 采样$x^∗ \sim Q(x|x_{t−1})$；
    2. 从均匀分布中采样$u \sim U(0,1)$；
    3. if u < $\alpha(x_{t−1},x^∗)$：接受转移，令$x_t = x^∗$；
    4. else: 不转移，令$x_t = x_{t−1}$；
  5. 返回采集到的一系列x；

似乎很令人惊讶，大名鼎鼎的MCMC竟然用程序来实现简单得不能再简单，但是采样的效果确实十分优秀的，让我们看一下M-H算法对一个beta(5,12)分布的采样结果：

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta

# 转移矩阵用正态分布
N = 5000
xs = [i for i in range(N)]
xs[0] = random.random()
p = beta(5, 12)
for i in range(1, N):
    x_star = random.uniform(0,1)    # 从U[0,1]中采样
    alpha = min(1, p.pdf(x_star) / p.pdf(xs[i-1]))
    u = random.uniform(0,1)
    xs[i] = x_star if u < alpha else xs[i-1]

plt.hist(xs, 50, density=1)
x = np.linspace(0, 1, 200)
plt.plot(x, [p.pdf(i) for i in x])
plt.show()

Gibbs

Gibbs适用于更高维的分布的采样，并且用一句话概括它的思想就是：固定其他维度上的随机变量，那么剩下那一维上的随机变量仍然满足细致平衡。
首先还是需要寻找合适的细致平衡条件。假设现在要采集二维联合分布$\pi(x_1, x_2)$，该分布上有两个点$A(x_1^1), x_2^1)$和B(x_1^1, x_2^2)，可以看到它们的第一个维度是相同的，此时$\pi(x_1^1,x_2^1) \neq \pi(x_1^1, x_2^2)$。这个是显而易见的，那么如何让其相等呢？还是用构造的方式。因为$\pi(x_1, x_2) = \pi(x_1)\pi(x_2|x_1)$，所以只要等式左右两边相互乘以“对方有而自己没有”的那一项，等式就成立了，即

$$
\pi(x_1^1, x_2^1)\pi(x_2^2|x_1^1) = \pi(x_1^1,x_2^2)\pi(x_2^1|x_1^1)
$$

即

$$
\pi(A)\pi(x_2^2|x_1^1) = \pi(B)\pi(x_2^1|x_1^1)
$$

因为点A和点B的关系是某一个维度相同，经过上面的构造得到等式的过程，我们可以发现，在二维空间里的两个点只要一个维度相同，那么这两个点的相互转化满足细致平衡。

就比如A和B，A和C之间存在细致平衡转化关系，而B和C之间就不存在。

那么Gibbs如何采样呢？以二维采样为例

  1. 给定一个平稳分布$\pi(x_1,x_2)$，给定转移次数n；
  2. 随机初始化一个点$X_0(x_1^0,x_2^0)$；
  3. for t = 1 to n：
    1. 固定$x_1$，采样$x_2$，得到$x_2^{t+1} \sim \pi(x_2|x_1^t)$；
    2. 固定$x_2$，采样$x_1$，得到$x_1^{t+1} \sim \pi(x_1|x_2^t)$；
    3. 得到新采样点$X_{t+1}(x_1^{t+1}, x_2^{t+1})$；
  4. 返回采集到的一系列X；

也是很简单的过程。

参考