蒙特卡洛采样【拒绝采样、重要性采样】

---

蒙特卡洛采样（Monte Carlo Simulation），蒙特卡洛方法又称“随机抽样方法”，和一般数值计算方法有本质区别的计算方法，是一种近似推断的方法，属于试验数学的分支，利用随机数进行统计试验，以求得统计特征，比如通过采样大量粒子的方法来求解期望、均值、面积、积分等问题。

本文介绍直接采样、接受拒绝采样与重要性采样三种蒙特卡洛采样方法：

• 直接采样最简单，但是需要已知累积分布函数(CDF)f(x)的形式，并且好求逆函数，得到服从分布f(x)的样本.
• 接受拒绝采样是给出一个提议分布$p(x)$(概率密度函数，PDF)，对满足原分布的粒子进行采样，得到服从p(x)（PDF)分布的样本.
• 重要性采样则是给予每个粒子不同的权重，计算f(x)在概率密度函数为p(x)下的均值:$\mathbb{E}[f(x)]=\int f(x)p(x)dx$.

一、均匀分布采样

首先我们介绍均匀分布采样，这是其他采样方法的基础。均匀分布Uniform(0,1)的样本是相对容易生成的，主要原理是线性同余随机数生成器。 通过线性同余发生器(LCG)可以生成伪随机数，我们用确定性算法生成[0,1]之间的伪随机数序列后， 这些序列的各种统计指标和均匀分布Uniform(0,1)的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

线性同余随机数生成器如下: $$x_{n+1}=(a{x}_n+c)\mathrm{mod}m$$式中a，c，m是数学推导出的合适的常数。这种算法产生的下一个随机数完全依赖当前的随机数，所以当随机数序列足够大的时候，随机数也会出现重复子序列的情况，这也是为什么说计算机产生的随机数是伪随机数。

总的来说，服从均匀分布采样的粒子我们可以很容易获得，均匀分布采样也是后面几种采样的基础。后文我们都假设能够直接从计算机获得服从均匀分布采样的粒子。

二、直接采样

直接采样的方法是根据累积分布函数进行采样，得到服从某个概率分布的样本。对一个已知概率密度函数(比如均匀分布）与累积概率密度函数的概率分布，我们可以直接从**累积分布函数（CDF）**进行采样。

原理：

实例

设$x$服从累积分布函数为$F(x)=1-e^{-x}$(可验证是单调不减，且积分为1的函数）的分布，则可以通过逆变换的方法对$F(x)$直接采样，产生服从F(X)分布的样本X.
$$\begin{gathered} 令y=1-e^{-x}\to e^{-x}=1-y \\ 两边求对数可得:x=-ln(1-y) \end{gathered}$$则$F^{-1}(x)=-ln(1-x)$，令$x_i$为均匀分布样本，则$X_i=-ln(1-x_i)$为服从累积密度函数为F(x)分布的样本.

三、拒绝接受采样

拒绝接受采样的目的仍然是得到服从某个概率分布的样本，不过这种方法是直接利用概率密度函数(PDF)得到样本。如下图所示，$p(x)$是我们希望采样的分布，$q(x)$是我们提议的分布(proposal distribution)，$q(x)$分布比较简单，令$kq(x)>p(x)$，我们首先在$kq(x)$中按照直接采样的方法采样粒子，接下来判断这个粒子落在途中什么区域，对于落在蓝线以外的粒子予以拒绝，落在蓝线下的粒子接受，最终得到符合$p(x)$的N个粒子。

拒绝接受采样的基本步骤：

1. 生成服从q(x)的样本$\to x_i$
2. 生成服从均匀分布U(0,1)的样本——$u_i$
3. 当$q(x_i)\cdot u_i<p(x_i)$,也就是二维点落在蓝线以下，此时接受$X_k=x_i$
4. 最终得到的$X_k$为服从p(x)的样本

实例

我们希望采样得到服从如下分布的样本
$$p(x)=e^{-\frac{x^2}{2}}\left({\sin }^2(6+x)+3{\cos }^2(x){\sin }^2(4x)+1\right)$$
这个函数的形式比较复杂，不容易积分得到CDF的形式，也就更不好对CDF求反函数，所以我们采用拒绝-接受采样，直接对PDF进行采样.

定义拒绝抽样函数：

def Reject_Sampling(p,q_x):
    """对p(x)进行拒绝采样
    Args:
        p : 目标PDF函数
        q_x : 选取的简单概率密度函数
        
    Returns:
        list: 采样样本
    """
    size = 1e06   #初始采样点数目
    k=10          #系数
    sample=[]
    for _ in range(int(size)):      #不建议用for循环，速度慢，这样写比较好理解
        xi=np.random.normal(0,1)
        qi=k*q_x.pdf(xi)   #乘上系数
        ui=np.random.uniform(0,1)
        pi=p(xi)

        #判断
        if qi*ui<pi:
            sample.append(xi)

    return sample 

实验并作图：

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

def f(x):
    px=np.exp(-(x**2)/2)*(np.sin(6+x)**2+3*(np.cos(x)**2)*(np.sin(4*x)**2)+1)
    return px
  
q_x = stats.norm(0, 1)
RS=Reject_Sampling(f,q_x)   #抽样数据集

fig,ax=plt.subplots(figsize=(10,5))
ax.hist(RS,bins=2000,density = True, stacked=True)  #画出抽样分布
ax.set_xlabel('x', fontsize=15)
ax.set_ylabel("f(x)", fontsize=15)
ax.set_title("Rejection Sampling", fontsize=20)

我们发现采样的样本的分布图和目标函数分布图基本一致。

四、重要性采样

1 目的

重要性采样的目的：求一个函数f(x)在p(x)分布下的期望，即
$$\mathbb{E}[f(x)]=\int f(x)p(x)dx$$
当$p(x)$很复杂、不解析、积分不好求时，可以通过重要性采样来计算。比如当$f(x)=x$时，我们可以算采样样本在$p(x)$分布下的期望.

2 原理

首先, 当我们想要求一个函数 $f(x)$ 在区间 $[a,b]$ 上的积分 ${\int }_a^{b}f(x)dx$ 时有可能会面临一个问题, 那就是积分曲线难以解析, 无法直接求积分。这时候我们可以采用一种估计的方式, 即在区 间 $[a,b]$ 上进行采样: { x 1 , x 2 … , x n } \left\{x_{1}, x_{2} \ldots, x_{n}\right\} {x1​,x2​…,xn​}, 值为 { f ( x 1 ) , f ( x 2 ) , … , f ( x n ) } \left\{f\left(x_{1}\right), f\left(x_{2}\right), \ldots, f\left(x_{n}\right)\right\} {f(x1​),f(x2​),…,f(xn​)} .

如果采样是均匀的, 即如下图所示:

那么显然可以得到这样的估计:
$${\int }_a^{b}f(x)dx\approx \frac{b-a}{N}\sum _{i=1}^{N}f\left(x_i\right)，\text{\hspace{1em}(1)}$$

在这里 $\frac{b-a}{N}$ 可以看作是上面小长方形的底部的 “宽”, 而 $f\left(x_i\right)$ 则是坚直的 “长”。

上述的估计方法随着取样数的增长而越发精确，那么有什么方法能够在一定的抽样数量基础上来增加准确度，减少方差呢？比如x样本数量取10000，那么显然在f(x)比较大的地方，有更多的$x_i$，近似的积分更精确，如下图所示，在圆圈部分$x_i$样本应该更多，所以$x_i$就不用按均匀分布进行采样。

并且原函数 $f(x)$ 也许本身就是定义在一个分布之上的, 我们定义这个分布为 $p(x)$, 这个分布可能比较复杂，我们无法直接从 $p(x)$ 上进行采样, 那么我们可以另辟蹊径重新找到一个更加简明的分布 $q(x)$, 从它进行采样, 希望间接地求出 $f(x)$ 在分布 $p(x)$ 下的期望，这就是重要性采样的主要思想.

搞清楚了这一点我们可以继续分析了，下面我们从$p(x)$是否归一化两种情况分别进行讨论.

（2.1） 若$p(x)$归一化

假设$p(x)$已经归一化，即$\int p(x)=1$，那么我们知道函数 $f(x)$ 在概率分布 $p(x)$ 下的期望为:
$$\mathbb{E}[f(x)]={\int }_{x}p(x)f(x)dx\text{\hspace{1em}(2)}$$

当$p(x)$比较复杂时，这个期望的值我们无法直接得到, 因此我们需要借助 $q(x)$ 来进行采样, $q(x)$ 可以选取简单的分布，比如设$q(x)$为均匀分布，当我们在 $q(x)$ 上采样得到 { x 1 , x 2 , … , x n } \left\{x_{1}, x_{2}, \ldots, x_{n}\right\} {x1​,x2​,…,xn​} (即$x_i$服从$q(x)$分布），由公式（1）我们可以得到 $f$ 在 $q(x)$ 下的期望为：
$$\mathbb{E}[f(x)]={\int }_{x}q(x)f(x)dx\approx \frac{1}{N}\sum _{i=1}^{N}f\left(x_i\right).\text{\hspace{1em}(3)}$$
上面这个式子就简单很多了，只要我们得到$x_i$然后代入$f(x)$然后求和就行了，而且均匀分布的样本$x_i$很容易获得。接着我们来考虑原问题，对式(2)进行改写, 即： $p(x)f(x)=q(x)\frac{p(x)}{q(x)}f(x)$, 我们可以得到:
$$\mathbb{E}[f(x)]={\int }_{x}q(x)\frac{p(x)}{q(x)}f(x)dx$$
这个式子我们可以看作是函数 $\frac{p(x)}{q(x)}f(x)$ 定义在分布 $q(x)$ 上的期望, 当我们在 $q(x)$ 上采样 { x 1 , x 2 , … , x n } \left\{x_{1}, x_{2}, \ldots, x_{n}\right\} {x1​,x2​,…,xn​} (服从$q(x)$分布)，由（3）式，我们可以得到 $f$ 的期望估计如下

$$\begin{array}{rl}\mathbb{E}[f(x)]& \approx \frac{1}{N}\sum _{i=1}^N\frac{p\left(x_i\right)}{q\left(x_i\right)}f\left(x_i\right)\\ & =\frac{1}{N}\sum _{i=1}^N{w}_{i}f\left(x_i\right)\end{array}\text{\hspace{1em}(4)}$$
在这里 $w_i=\frac{p\left(x_i\right)}{q\left(x_i\right)}$就是重要性权重。

（2.2）若$p(x)没有归一化$

上面的讨论是假设$p(x)$已经完成归一化了，也就是$\int p(x)=1$,假如$p(x)$没有归一化，那么我们可以先对$p(x)$进行归一化：
$$\begin{array}{rl}\mathbb{E}[f(x)]& =\int f(x)\frac{p(x)}{\int p(x)dx}dx\\ & =\frac{\int f(x)p(x)dx}{\int p(x)dx}\\ & =\frac{\int f(x)\frac{p(x)}{q(x)}q(x)dx}{\int \frac{p(x)}{q(x)}q(x)dx}.\end{array}$$
而分子分母可分别得到如下两式（下面两式约等于都利用$q(x)$是均匀分布的假设）：
$$\begin{array}{rl}\int f(x)\frac{p(x)}{q(x)}q(x)dx& \approx \frac{1}{n}\sum _{i=1}^n{W}_{i}f\left(x_i\right),\\ \int \frac{p(x)}{q(x)}q(x)dx& \approx \frac{1}{n}\sum _{i=1}^n{W}_i.\end{array}$$
其中$W_i=\frac{p(x_i)}{q(x_i)}$，则最终可得$\mathbb{E}[f(x)]$:
$$\begin{array}{r}\mathbb{E}[f(x)]\approx \sum _{i=1}^n{w}_{i}f\left(x_i\right),w_i=\frac{W_i}{{\sum }_{i=1}^n{W}_i}\end{array}\text{\hspace{1em}(5)}$$

3 实例

我们希望求如下函数
$$p(x)=e^{-\frac{x^2}{2}}\left({\sin }^2(6+x)+3{\cos }^2(x){\sin }^2(4x)+1\right),$$
的期望$\mathbb{E}[p(x)]={\int }_{-\infty }^{\infty }xp(x)dx$，这个函数的形式比较复杂，积分有点痛苦，所以我们尝试用重要性采样来得到这个期望值，注意在这个问题中$f(x)=x$，$f(x)$服从$p(x)$分布.

这个问题$p(x)$的形式我们知道，那么只要有服从均匀分布的样本$x_i$，那么$p(x_i)$我们是直接能算出来的，$p(x_i)$有了，重要性权重$w_i$也就有了，那么利用式$(5)$期望的估计也就能算了，下面是对应的Python代码实现.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

def p(x):
    """ 概率密度函数，作为被积函数 """
    px = np.exp(-(x**2)/2) * (np.sin(6 + x)**2 + 3 * (np.cos(x)**2) * (np.sin(4 * x)**2) + 1)
    # px = 1 / 2 * np.exp(- (x-1) ** 2 /2) # 正太分布可以用来验证，均值是否为1
    return px

def f(x):
    return x

q_x = stats.norm(0, 1)  # 这里我们取 q(x) 为标准正态分布



def importance_sampling(f,p, n=1000):
    """ 使用重要性采样来估计函数的期望。
    Args:
        f (function): 目标函数 f(x)
        p (function): f服从 p(x) 的分布,p(x)为概率密度函数
        n (int): 采样的数量，默认为 1000。

    Returns:
        float: 估计的期望值 E[f(X)]。
    """
    total_weighted_value = 0
    total_weight = 0
    for _ in range(n):
        x_i = np.random.normal(0, 1)    # 从 q(x) 中采样
        W_i = p(x_i) / q_x.pdf(x_i)     # 计算权重 W_i = p(x_i) / q(x_i)
        total_weighted_value += f(x_i) * W_i
        total_weight += W_i

    expected_value = total_weighted_value / total_weight
    return expected_value


estimated_mean = importance_sampling(f, p, 100000)
print("E(X) =", estimated_mean)

最终我们得到期望值为$\mathbb{E}[p(x)]=-0.0343$，我们还可以用上面拒绝接受采样的结果验证，拒绝接受采样是得到服从$p(x)$分布的样本，那么我们可以用样本均值来估计期望，即：
$$\mu =\frac{x_1+x_2+\cdots +x_n}{n}$$

mean = np.mean(RS)
print("mean:",mean)

我们得到$\mu =-0.03178$，可以发现两个的结果一致（采样有随机性，每次结果都有细微差别）。我们还可以用正太分布来验证，比如代码里我写的$p(x)=\frac{1}{2}{e}^{-\frac{(x-1{)}^2}{2}}$，显然这个正太分布的期望为1，我们用重要性采样得到的结果为$\mu =1.0020$，可以看到基本上可以精确到小数点后三位，相对误差为0.206%.

4 方差分析

虽然我们可以把 $p$ 换成任何的 $q$，但是在实现上，$p$ 和 $q$ 的差距不能太大。差距太大，会有一些问题。虽然通过上面的推导我们知道
$${\mathbb{E}}_{x\sim p}[f(x)]={\mathbb{E}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right],\text{\hspace{1em}(6)}$$
但如果不是计算期望值，而是计算方差，${\mathrm{Var}}_{x\sim p}[f(x)]$ 和 ${\mathrm{Var}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]$ 是不一样的。两个随机变量的平均值相同，并不代表它们的方差相同。 我们可以将 $f(x)$ 和 $f(x)\frac{p(x)}{q(x)}$ 代入方差的公式
$$\mathrm{Var}[X]=\mathbb{E}[X^2]-(\mathbb{E}[X]{)}^2,$$
可得
$${\mathrm{Var}}_{x\sim p}[f(x)]={\mathbb{E}}_{x\sim p}[f(x{)}^2]-{\left({\mathbb{E}}_{x\sim p}[f(x)]\right)}^2,\text{\hspace{1em}(7)}$$
以及（下面的转换注意利用公式(6)）
$$\begin{array}{rl}{\mathrm{Var}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]& ={\mathbb{E}}_{x\sim q}\left[{\left(f(x)\frac{p(x)}{q(x)}\right)}^2\right]-{\left({\mathbb{E}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]\right)}^2\\ & ={\mathbb{E}}_{x\sim q}\left[\frac{p(x)}{q(x)}f(x{)}^2\frac{p(x)}{q(x)}\right]-{\left({\mathbb{E}}_{x\sim p}[f(x)]\right)}^2\\ & ={\mathbb{E}}_{x\sim p}\left[f(x{)}^2\frac{p(x)}{q(x)}\right]-{\left({\mathbb{E}}_{x\sim p}[f(x)]\right)}^2\end{array}\text{\hspace{1em}(8)}$$
我们可以看到 ${\mathrm{Var}}_{x\sim p}[f(x)]$ 和 ${\mathrm{Var}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]$ 的差别在于第二项是不同的，${\mathrm{Var}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]$ 的第一项多乘了 $\frac{p(x)}{q(x)}$。如果 $\frac{p(x)}{q(x)}$ 差距很大，$f(x)\frac{p(x)}{q(x)}$ 的方差就会很大。所以理论上它们的期望值一样，也就是说，我们只要对分布 $p$ 来采样足够多次，对分布 $q$ 来采样足够多次，得到的结果会是一样的。但是如果我们采样的次数不够多，因为它们的方差差距是很大的，所以我们就有可能得到差距非常大的结果。

例如，当 $p(x)$ 和 $q(x)$ 差距很大时，就会有问题。如下图所示(图源¹），假设蓝线是 $p(x)$ 的分布，绿线是 $q(x)$ 的分布，红线是 $f(x)$。

如果我们要计算 $f(x)$ 的期望值：

1. 从分布 $p(x)$ 做采样，显然 ${\mathbb{E}}_{x\sim p}[f(x)]$ 是负的。这是因为左边区域 $p(x)$ 的概率很高，所以采样会到这个区域，而 $f(x)$ 在这个区域是负的，所以理论上这一项算出来会是负的。
2. 接下来我们改成从 $q(x)$ 采样，因为 $q(x)$ 在右边区域的概率比较高，所以如果我们采样的点不够多，可能只会采样到右侧。如果我们只采样到右侧，可能 ${\mathbb{E}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]$ 是正的。我们这边采样到这些点，去计算它们的 $f(x)\frac{p(x)}{q(x)}$ 都是正的。我们采样到这些点都是正的，取期望值以后也都是正的，这是因为采样的次数不够多。

当然如果我们采样足够多，左边虽然概率很低，但也有可能被采样到。假设我们好不容易采样到左边的点，因为左边的点的 $p(x)$ 和 $q(x)$ 是差很多的，这边 $p(x)$ 很大，$q(x)$ 很小。$f(x)$ 好不容易终于采样到一个负的，这个负的就会被乘上一个非常大的权重，这样就可以平衡刚才那边一直采样到正的值的情况。最终我们算出这一项的期望值，终究还是负的。 但前提是我们要采样足够多次，这件事情才会发生。 所以当采样次数不够多时，${\mathbb{E}}_{x\sim p}[f(x)]$ 与 ${\mathbb{E}}_{x\sim q}\left[f(x)\frac{p(x)}{q(x)}\right]$ 可能就有很大的差距。这就是重要性采样的问题。

五、参考资料

1. 重要性采样
2. 蒙特卡洛方法
3. 一文看懂蒙特卡洛方法

---

1. 王琦，杨毅远，江季，Easy RL：强化学习教程，人民邮电出版社，https://github.com/datawhalechina/easy-rl. ↩ ↩︎