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

---

蒙特卡洛方法（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$，则可以算$\mathbb{E}[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=\frac{b-a}{N}\sum _{i=1}^{N}f\left(x_i\right)$$在这里 $\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)$ 下的期望。

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

搞清楚了这一点我们可以继续分析了。首先我们知道函数 $f(x)$ 在概率分布 $p(x)$ 下的期望为:
$$\mathbb{E}[f(x)]={\int }_{x}p(x)f(x)dx(1)$$

但是这个期望的值我们无法直接得到, 因此我们需要借助 $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)分布）后，那么我们可以估计 $f$ 在 $\mathbf{q}(\mathbf{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).$$
上面这个式子就简单很多了，只要我们得到$x_i$然后代入$f(x)$然后求和就行了，而且均匀分布的样本$x_i$很容易获得。接着我们来考虑原问题，对式(1)进行改写, 即： $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)分布)，可以估计 $f$ 的期望:

$$\begin{array}{rl}\mathbb{E}[f(x)]& =\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}$$
在这里 $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}$$

（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$.

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%.

五、参考资料

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