采样（三）：重要性采样与接受拒绝采样

本文介绍重要性采样与接受拒绝采样，并给出使用Laplace分布作为参考分布，采样正态分布的例子。

重要性采样（importance sampling）

在概率分布$p(x)$下计算$f(x)$均值，

$\begin{align} \int_{\Omega} p(x)f(x)dx &= \mathbb{E}_{x\sim p(x)}[f(x)] \newline &\approx \frac{1}{n} \sum_{i=1}^{n} f(x_{i}) \end{align}$

其中$x_{i} \sim p(x)$，然而在很多情况下$p(x)$并不是简单的分布，无法像均匀分布、正态分布一样直接采样。这种情况下重要性采样引入一个方便采样的分布$q(x)$，把均值计算改写为，

$\begin{align} \int_{\Omega} p(x)f(x) dx &= \mathbb{E}_{x\sim p(x)}[f(x)] \newline &= \int_{\Omega} q(x) \frac{p(x)}{q(x)} f(x) dx \newline &= \mathbb{E}_{x \sim q(x)} \Big[ \frac{p(x)}{q(x)} f(x) \Big] \newline &\approx \frac{1}{n}\sum_{i=1}^{n} \frac{p(x_i)}{q(x_i)} f(x_i) \end{align}$

其中$x_{i} \sim q(x)$比较容易，而$\frac{p(x_i)}{q(x_i)} f(x_i)$的计算没有任何难度。注意到$q(x_i) \ne 0$。$ \frac{p(x_i)}{q(x_i)}$可以看着是样本的权重，这正是“重要性”来源。

$p(\boldsymbol{x})$为不容易采样的分布，计算$f(\boldsymbol{x})$的均值有，

$\begin{align} \mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}[f(\boldsymbol{x})] &= \sum_{\boldsymbol{x}} p(\boldsymbol{x}) f(\boldsymbol{x}) \newline &= \sum_{\boldsymbol{x}} q(\boldsymbol{x})\frac{p(\boldsymbol{x})}{q(\boldsymbol{x})}f(\boldsymbol{x}) \newline &= \mathbb{E}_{\boldsymbol{x}\sim q(\boldsymbol{x})}\left[\frac{p(\boldsymbol{x})} {q(\boldsymbol{x})}f(\boldsymbol{x})\right] \newline &= \mathbb{E}_{\boldsymbol{x}\sim q(\boldsymbol{x})}\left[w(\boldsymbol{x}) \cdot f(\boldsymbol{x})\right] \end{align}$

对每个样本加权，

$w(\boldsymbol{x}) = \frac{p(\boldsymbol{x})}{q(\boldsymbol{x})}$

可以验证，权重符合归一性，

$\mathbb{E}_{\boldsymbol{x} \sim q(\boldsymbol{x})} \left[ w(\boldsymbol{x}) \right] = \sum_{\boldsymbol{x}} q(\boldsymbol{x}) \frac{p(\boldsymbol{x})}{q(\boldsymbol{x})} = \sum_{\boldsymbol{x}} p(\boldsymbol{x}) = 1$

接受-拒绝采样（accept-reject sampling）

假设有概率分布$p(x)$，称为目标分布，需要从中采样往往比较难，于是引入较简单的易于采样的分布$q(x)$，称为参考分布或建议分布，并定义如下函数，

$\alpha(x) = \frac{p(x)}{q(x)} \in [0, 1]$

现在采样样本$x_{i} \sim q(x)$，采样随机数$\varepsilon_{i} \sim U[0,1]$，如果$\varepsilon_i \le \alpha(x_{i})$，则接受样本$x_i$否则拒绝该样本，继续以上流程采样新样本，如此下去获得的样本集$x_1, x_2, \dots, x_n$，服从目标分布$p(x)$。

$q(x)$为建议分布，$p(x)$为目标分布，满足如下关系，

$M \cdot q(x) \ge p(x)$

即$M\times q(x)$要包络$p(x)$。

$\begin{align} 1 &= \int p(x) dx \\ &\approx \sum_{i=1}^{n} \alpha(x_i) q(x_{i}) \end{align}$

因此，归一化后可获得$p(x)$的采样集，

$\frac{ \alpha(x_{i})q(x_{i})}{\displaystyle \sum_{i=1}^{n} \alpha(x_{i})q(x_{i})}$

那么接受概率是多少呢？直觉上这明显取决于$M \cdot q(x)$ 包络 $p(x)$ 的质量，考虑到接受样本需要满足，

$\varepsilon_i \le \frac{p(x_i)}{M \times q(x_i)}, x_i \sim q(x), \varepsilon_i \sim U[0, 1]$

注意下式是从$q(x)$中采样，接受概率为，

$\begin{align} P(\text{接受}) &= P(U \lt \frac{p(x)}{M\cdot q(x)}) \newline &= \int q(x) \frac{p(x)}{M \cdot q(x)} dx \newline &= \frac{1}{M} \end{align}$

因此，高质量的采样，需要 $M$ 尽可能小，否则从$q(x)$采样出来的样本大部分都会被拒绝掉，提高了程序运行时的空转时间。

利用贝叶斯公式可以知道最后的采样分布（接受样本条件下得到的分布）为，

$q(x|\text{接受}) = \frac{q(x) \times P(\text{接受}| X = x)}{P(\text{接受}) } = p(x)$

接下来我们给出一个例子，使用Laplace分布作为参考分布，采样正态分布。

拉普拉斯分布采样正态分布

以拉普拉斯分布为参考分布，采样目标分布为正态分布$p(x)$

$f(x;\mu ,\sigma )=\frac {1}{\sqrt {2\pi }\sigma } e^{-\frac {(x-\mu )^{2}}{2\sigma ^{2}}} \Big|_{\mu = 0}$

拉普拉斯分布的概率密度$q(x)$，

$f(x\mid \mu ,b)={\frac {1}{2b}}\exp \left(-{\frac {|x-\mu |}{b}}\right)\, \Big|_{\mu = 0}$

均值和方差分别为$\mu, 2b^{2}$。

计算$M \cdot q(x)$ 包络 $p(x)$ 中$M$的最小值，

$\frac {1}{\sqrt {2\pi }\sigma } \exp\left({-\frac {x^{2}}{2\sigma ^{2}}} \right) \le M \times {\frac {1}{2b}}\exp \left(-{\frac {|x |}{b}}\right)\,$

取特殊值，解得包络值的下界，

$M \ge \frac{2b}{\sqrt{2\pi} } \cdot \exp(\frac{1}{2b^2})$

取等号获得最优包络的$M$值。

均匀分布采样正态分布

类似地，以均匀分布为参考分布，

$\frac {1}{\sqrt {2\pi }\sigma } \exp\left({-\frac {x^{2}}{2\sigma ^{2}}} \right) \le M \times \frac{1}{b-a}$

这里正态分布参数$\sigma=1$，认为$x \in (-5\sigma , 5\sigma)$足够（正态分布的$k\sigma$法则），解得

$M \ge \frac{10}{\sqrt{2 \pi}}$

实现

Python实现如下，

import numpy as np

def accept_reject_sampling(pf, qf, q_gen, M):
    """
    pf:目标采样概率密度函数
    qf:参考分布概率密度函数
    q_gen:参考分布采样生成器
    M:使得pf(x)<=qf(x)*M，qf(x)*M称为pf(x)的包络函数
    """
    for x in q_gen():
        u = np.random.uniform()
        if u < pf(x) / (M * qf(x)):
            yield x