蒙特卡洛采样【python实例】

---

蒙特卡洛方法（Monte Carlo Simulation）是一种近似推断的方法，通过采样大量粒子的方法来求解期望、均值、面积、积分等问题。

蒙特卡洛对某一种分布的采样方法有直接采样、接受拒绝采样与重要性采样三种：

• 直接采样最简单，但是需要已知累积分布的形式，并且好求逆函数，得到p(x)分布的样本
• 接受拒绝采样是给出一个提议分布，对满足原分布的粒子接受，得到服从p(x)分布的样本
• 重要性采样则是给予每个粒子不同的权重，计算f(x)在概率分布p(x)下的均值

一、均匀分布采样

均匀分布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}$的分布，则可以通过逆变换的方法对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)分布的样本

三、拒绝接受采样

对于累积分布函数比较复杂，不好求反函数的分布，我们可以采用接受-拒绝采样。如下图所示，p(z)是我们希望采样的分布，q(z)是我们提议的分布(proposal distribution)，q(z)分布比较简单，令kq(z)>p(z)，我们首先在kq(z)中按照直接采样的方法采样粒子，接下来判断这个粒子落在途中什么区域，对于落在蓝线以外的粒子予以拒绝，落在蓝线下的粒子接受，最终得到符合p(z)的N个粒子。

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

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

实例

定义拒绝抽样函数：

def Reject_Sampling(p,q_x):
    """对p(x)进行拒绝采样

    Args:
        p (_type_): 目标CDF函数
        q_x (_type_): 选取的简单概率密度函数

    Returns:
        _type_: 采样样本
    """
    size = 1e+05   #初始采样点
    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):
    fx=np.exp(-(x**2)/2)*(np.sin(6+x)**2+3*(np.cos(x)**2)*(np.sin(4*x)**2)+1)
    return fx
  
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)分布下的期望，即$E(f(x))=\int f(x)p(x)dx$，当p(x)很复杂时，不解析，积分不好求时，可以通过重要性采样来计算。当$f(x)=x$，则可以算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)$ 也许本身就是定义在一个分布之上的, 我们定义这个分布为 $\pi (x)$, 我们无法直接从 $\pi (x)$ 上进行采样, 所以另辟蹊径重新找到一个更加简明的分布 $p(x)$, 从它进行取样, 希望间接地求出 $f(x)$ 在分布 $\pi (x)$ 下的期望。

（2.1） $\pi (x)归一化$

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

但是这个期望的值我们无法直接得到, 因此我们需要借助 $p(x)$ 来进行采样, $p(x)$ 可以选取简单的分布，比如设p(x)为均匀分布，当我们在 $p(x)$ 上采样 { x 1 , x 2 , … , x n } \left\{x_{1}, x_{2}, \ldots, x_{n}\right\} {x1?,x2?,…,xn?} (即$x_i$服从p(x)分布）后可以估计 $f$ 在分布 $p(x)$ 下的期望为：
$$E[f]={\int }_{x}p(x)f(x)dx\approx \frac{1}{N}\sum _{i=1}^{N}f\left(x_i\right)。$$接着我们可以对式子进行改写, 即： $\pi (x)f(x)=p(x)\frac{\pi (x)}{p(x)}f(x)$, 所以我们可以得到:
$$E[f]={\int }_{x}p(x)\frac{\pi (x)}{p(x)}f(x)dx$$这个式子我们可以看作是函数 $\frac{\pi (x)}{p(x)}f(x)$ 定义在分布 $p(x)$ 上的期望, 当我们在 $p(x)$ 上采样 { x 1 , x 2 , … , x n } \left\{x_{1}, x_{2}, \ldots, x_{n}\right\} {x1?,x2?,…,xn?} (服从p(x)分布)，可以估计 $f$ 的期望:

$$\begin{array}{rl}E[f]& =\frac{1}{N}\sum _{i=1}^N\frac{\pi \left(x_i\right)}{p\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{\pi \left(x_i\right)}{p\left(x_i\right)}$就是重要性权重。

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

$$\begin{array}{rl}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}$$而分子分母可分别得到：
$$\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}$$则最终可得E(f(x)):
$$\begin{array}{r}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}$$其中$W_i=\frac{p(x_i)}{q(x_i)}$.

（3）实例

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
def f(x):   #概率密度函数
    fx=np.exp(-(x**2)/2)*(np.sin(6+x)**2+3*(np.cos(x)**2)*(np.sin(4*x)**2)+1)
    return fx
  
q_x = stats.norm(0, 1) #取q(x)为正太分布
 
# 从 q(x) 进行重要性采样
def Importance_Sampling(f):
    """对f(x)进行重要性采样
    Args:
        f (_type_): 传入f(x)函数

    Returns:
        _type_: 返回采样结果
    """
    EX_list = []
    W=0
    for i in range(n):
        
        x_i = np.random.normal(0, 1)     #x为正太分布的采样样本
        w_i=f(x_i) / (1*q_x.pdf(x_i))    #权重  
        value = x_i*w_i
        W+=w_i
        EX_list.append(value)

    W=W/n
    return EX_list/W        #归一化,并返回采样结果


IS=Importance_Sampling(f)

IS_Mean=np.mean(IS)   

print("E(X)=",IS_Mean)

参考资料

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