蒙特卡洛积分和重要性采样(Importance Sampling)

一、蒙特卡洛积分

• 蒙特卡洛积分概述：简而言之蒙特卡洛积分就是，在求定积分时，如果找不到被积函数的原函数，无法使用经典牛顿-莱布尼茨积分法得到定积分结果的。而蒙特卡洛积分方法利用一个随机变量对被积函数进行采样，并将采样值进行一定的处理可以得到定积分的一个近似值，当采样数量很高时，得到的近似值可以很好的近似原积分的结果。这样一来，我们就不用去求原函数的形式，就能求得积分的近似结果。

• 补充一些基础性公式：

  • 假设一连续型随机变量$X$的样本空间为$D$，其概率密度分布函数为$p(x)$，则其数学期望为：$$E(X)={\int }_{D}xp(x)dx$$
  • 若另一连续随机变量$Y$满足$Y=f(X)$，则$Y$的数学期望为：$$E(Y)={\int }_{D}f(x)p(x)dx$$

• 蒙特卡洛积分方法基础形式

  • 现在 假设我们要计算一个定积分：$$A={\int }_a^{b}f(x)dx$$根据牛顿-莱布尼茨公式我们可以得到：$$A={\int }_a^{b}f(x)dx=F(b)-F(a)$$其中$F(x)$是$f(x)$的一个原函数。
  • 如果我们不知道或者无法求得原函数，我们该怎么计算这个定积分呢。那么就需要借助蒙特卡洛积分(Monte Carlo Integration)方法：
    • 首先我们可以在区间$[a,b]$上进行均匀采样得到：{X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​},样本对于的函数值为：{f(X1),…,f(XN)}\{f(X_1),…,f(X_N)\}{f(X1​),…,f(XN​)}
    • 然后再求和得到：$$F_N\approx \frac{b-a}{N}\sum _{i=1}^{N}f(X_i)$$
    • 用$F_N$作为A近似估值。
    • 这个和定积分的定义(黎曼积分)非常相似，只是在定积分的定义中{X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​}是不均匀采样得到而是对区间[a,b]均匀划分得到：{x1=a,…,xN=b}\{x_1=a,…,x_N=b\}{x1​=a,…,xN​=b}。根据定积分的定义如果划分的次数N趋于无穷大的时候$F_N=A$.如下图所示：
      

• 蒙特卡洛积分方法的正确性进行进一步分析：

  • {X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​}是通过均匀分布采样得到的,则$X_i$也是随机变量并且服从均匀分布，即：$X_i$~$U(a,b)$

  • 而$F_N$是{X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​}的函数那么$F_N$是一个样本统计量也是一个随机变量。

  • 那么现在我们计算一下$F_N$的数学期望：$$\begin{gathered} E[F_N]=E[\frac{b-a}{N}\sum _{i=1}^{N}f(X_i)]=\frac{b-a}{N}\sum _{i=1}^{N}E[f(X_i)] \\ =\frac{b-a}{N}\sum _{i=1}^N{\int }_a^{b}f(x)p(x)dx=\frac{b-a}{N}\sum _{i=1}^N{\int }_a^{b}f(x)\frac{1}{b-a}dx \\ =\frac{b-a}{N}\frac{1}{b-a}\sum _{i=1}^N{\int }_a^{b}f(x)dx=\frac{1}{N}\sum _{i=1}^N{\int }_a^{b}f(x)dx \\ ={\int }_a^{b}f(x)dx=A \end{gathered}$$其中$p(x)$是$[a,b]$上均匀分布的概率密度函数

  • 由上面的推到可知$F_N$是一个样本统计量并且其期望为$A$，根据大数定律，$F_N$依概率收敛与$A$，即随着$N$的增加，理论上$F_N$将越逼近$A$的值。也可以称$F_N$是$A$的无偏估计量。所以$F_N$可以视作$A$的一个近似值。

• 蒙特卡洛积分方法做一个一般性的推广：

  • 上面的推导过程是在区间$[a,b]$上进行均匀采样得到的结论，那么如果我们在区间上按照概率密度函数$p(x)$进行采样得到{X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​}上面的结论是否成立呢。答案是肯定的。下面我们就来推导这种情况：
  • 首先我们按照概率密度函数$p(x)$区间$[a,b]$上进行采样得到样本集{X1,…,XN}\{X_1,…,X_N\}{X1​,…,XN​}
  • 再构造新的$F_N$函数：$$F_N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}$$
  • 用$F_N$作为定积分$A$的近似估计
  • 计算一些$F_N$的数学期望：$$\begin{gathered} E[F_N]=E[\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}]=\frac{1}{N}\sum _{i=1}^{N}E[\frac{f(X_i)}{p(X_i)}] \\ =\frac{1}{N}\sum _{i=1}^N{\int }_a^b\frac{f(X_i)}{p(X_i)}p(X_i)dx=\frac{1}{N}\sum _{i=1}^N{\int }_a^{b}f(x)dx \\ ={\int }_a^{b}f(x)dx=A \end{gathered}$$
  • 到这里我们发现其实前推导的那种情况是上面这种情况$p(x)$为均匀分布的一种特殊情况：
    • 若$p(x)$是$[a,b]$上均匀分布，则它的表达式为：$$p(x)=\{\begin{array}{lr}1& 0\le x\le 1\\ 0& 其他\end{array}$$
    • 则$F_N(x)$的表达式为：$$\begin{gathered} F_N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}=\frac{1}{N}\sum _{i=1}^N(b-a)f(X_i) \\ =\frac{b-a}{N}\sum _{i=1}^{N}f(X_i) \end{gathered}$$

• 经过上面的堆到，我们可以得到蒙特卡洛积分方法如下：
  $${\int }_{D}f(x)dx=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}$$

• 蒙特卡洛积分法的收敛性和收敛速度分析

  • 求$F_N$的方差:$${\sigma }^2(F_N)={\sigma }^2[\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}]=\frac{1}{N^2}\sum _i^N{\sigma }^2[\frac{f(X_i)}{p(X_i)}]$$令：$$Y=\frac{f(X)}{p(Y)}$$则：$${\sigma }^2(F_N)=\frac{1}{N^2}\sum _i^N{\sigma }^2[Y_i]=\frac{1}{N^2}(N{\sigma }^2[Y])=\frac{1}{N}{\sigma }^2[Y]$$则：$$\sigma [F_N]=\frac{1}{\sqrt{N}}\sigma [Y]$$
  • 由此可知当采样个数$N$趋向无穷大时$\sigma [F_N]=0$，即$F_N$趋向其均价$A$；这也证明蒙特卡洛积分法的收敛性。
  • 由方差公式可以知道蒙特卡洛积分法收敛速度和稳定性由$Y$决定，如果$Y_i=\frac{f(X_i)}{p(X_i)}$ 因不同$X_i$的取值变化地越剧烈，就会造成$Y$的方差$\sigma [Y]$较大，则会造成估计值的收敛速度越慢, 这就告诉我们，若$p(x)$与$f(x)$越接近(即若他们的形状越接近)，则收敛更快更稳定。最理想的情况是$p(x)=f(x)$,则$Y=1$并且$\sigma [Y]=0$.
  • 若$p(x)$与$f(x)$很接近(甚至$p(x)=f(x)$,那f(x)大的地方$p(x)$也会更大，那么我们更加$p(x)$函数进行采样时就会对这个个地方采样更多，即对重要的地方采样更多。如下图所示：
    
    很明显在圆形区域的函数值对积分的贡献比方形区域要大很多(即圆形区域比比方形区域更重要)，所以我们可以在抽样的时候以更大的概率抽取圆形区域的样本，这就是重要性采样的思想.

• 利用蒙特卡洛方法，我们可以得到任意一个积分的结果，但是可能得不到精确值，我们得到的只是一个对理论值的估计，估计值与理论值之间的误差可以通过增加样本数来减小，但收敛速率仅为$O(\sqrt{N})$，也就是说，若想将误差降为现在的一半，我们需要再多计算4倍的计算量才可以达到。即便如此，原始的蒙特卡洛积分方法也不失为是一种经典有效的方法。

二、重要性采样

• 假设$X$为连续型随机变量，其概率密度函数为$\pi (x)$,$Y$也是一个随机变量，$Y$与$X$关系为$Y=f(X)$,现在我们需要求$Y$的期望：$$E_{\pi }(Y)={\int }_x\pi (x)f(x)dx$$
• 如果$\pi (x)$比较难求或者无法求解，我们需要寻找其他方法来估计$E_{\pi }(Y)$的值。那么就可以使用蒙特卡洛积分法来估计该积分。假设我们找到一个容易采样且与$\pi (x)f(x)$比较接近的分布$p(x)$,然后对$p(x)$进行采样得到样本序列：{xi,…,xN}\{x_i,…,x_N\}{xi​,…,xN​},根据蒙特卡洛积分法可得：$$E_{\pi }(Y)={\int }_x\pi (x)f(x)dx=\frac{1}{N}\sum _{i=1}^N\frac{\pi (x_i)}{p(x_i)}f(x_i)dx$$
• $\frac{\pi (x_i)}{p(x_i)}$就是就是重要性权重