通俗易懂的MonteCarlo积分方法(五)

通俗易懂的MonteCarlo积分方法(五)

​ 这次主要目的是想办法提高$MonteCarlo$的计算精度。

​ 如果我们要求解一个定积分:
$$J={\int }_a^{b}f(x)dx$$
​ 在求解这个定积分之前我们先假设我们可以产生在$[a,b]$上的具有给定密度函数$p(x)$的独立随机变量$\zeta$。

​ 那么这个$p(x)$在$[a,b]$上就满足以下条件:
$$\begin{gathered} p(x)>0 \\ {\int }_a^{b}p(x)dx=1 \end{gathered}$$
​ 如果我们设定随机变量$\eta =\frac{f(x)}{p(x)}$，那么$\eta$的均值如下:
$$E\eta ={\int }_a^b\frac{f(x)}{p(x)}p(x)dx=J$$
​ 那么$\eta$的方差如下:
$$D\eta =E{\eta }^2-(E\eta {)}^2={\int }_a^b\frac{f^2(x)}{p(x)}dx-J^2$$
​ 由我们之前讨论的辛钦大数定理，我们可以用$\widehat{J_N}=\frac{1}{N}{\sum }_{i=1}^N{\eta }_i=\frac{1}{N}{\sum }_{i=1}^N\frac{f({\zeta }_i)}{p({\zeta }_i)}$来近似描述$J$。

​ 我们由中心极限定理可知:
l i m n → ∞ P { ∣ J N ^ − J ∣ < 3 D η N } = 0.997 lim_{n\rightarrow \infty}P\{|\hat{J_N}-J|<\frac{3\sqrt{D\eta}}{\sqrt{N}}\} = 0.997 limn→∞​P{∣JN​^​−J∣<N ​3Dη ​​}=0.997
​ 由中心极限定理，误差$\delta ={\hat{J}}_N-J$随着$N$的增大逐渐减小，但是误差却是与$\frac{1}{\sqrt{N}}$成反比。这说明在$N$特别大时候，再增大$N$时数值误差就会减少的比较小。

​ 因此误差随着$N$的增加时，数值误差的改进就会特别小。这就是$MonteCarlo$方法最大的弊端。

​ 现在我们再考虑在提高精度，减少误差的另一种办法，就是通过：选取合适的$p(x)$来使得$\sqrt{D\eta }$尽可能的小。

​ 我们由柯西不等式可知:
$$({\int }_a^b|f(x)|dx{)}^2\le {\int }_a^b\frac{f^2(x)}{p(x)}dx{\int }_a^{b}p(x)dx={\int }_a^b\frac{f^2(x)}{p(x)}dx=E{\eta }^2=D\eta +J^2$$
​ 当且仅当$p(x)=\frac{|f(x)|}{{\int }_a^b|f(x)|dx}$时上述等式成立，此时$D\eta =0$成立。

​ 但是我们这个时候选取得$p(x)$不现实。因为我们要计算${\int }_a^b|f(x)|dx$，相比于${\int }_a^{b}f(x)dx$,本身就需要巨大得计算量。

​ 因此我们可以选择一个合适的算法来求较优的$p(x)$。比如说，我们可以考虑$p(x)$在$|f(x)|$的极值附近增加较快，而在不明显的地方取$p(x)$近似为常数。

​ 假设我们的$p(x)$已经选定了，由于$J$我们不知道时多少，因此我们用统计的方法估计$D\eta$。估计方法如下：
$$D\eta \approx {\hat{D}}_{\eta ,N}=\frac{1}{N-1}\sum _{i=1}^N(\frac{f({\zeta }_i)}{p({\zeta }_i)}{)}^2-{\hat{J}}_N^2$$
​ 而相关的${\hat{J}}_N$和${\hat{D}}_{\eta ,N}$可以由下面的算法决定：
$$\{\begin{array}{l}{\hat{J}}_0=0,{\hat{D}}_0=0\\ {\hat{J}}_{N+1}={\hat{J}}_N-\frac{1}{N+1}(\frac{f({\zeta }_{N+1})}{p({\zeta }_{N+1})}-{\hat{J}}_N)\\ {\hat{D}}_{N+1}={\hat{D}}_N-\frac{1}{N+1}(\frac{f^2({\zeta }_{N+1})}{p^2({\zeta }_{N+1})}-{\hat{D}}_N)\\ {\hat{D}}_{\eta ,N}={\hat{D}}_N-{\hat{J}}_N^2\end{array}$$
​ 设定一个偏差$ϵ$,这个算法在以下条件时有满足$0.997$的概率达到结束。
$$3\sqrt{\frac{D_{\eta ,N}}{N}}\le ϵ$$
​ 以上的算法的关键时产生服从我们设定的在$[a,b]$上任意分布$p(x)$的随机数。这是一个很复杂的问题，一般是用$MCMC$来产生分布为$p(x)$随机数。这个算法的$python$代码的实现过程下次再描述。