通俗易懂的Monte Carlo积分方法(四)

1.理论目标

​ 利用辛钦大数定律和中心极限定理对$MentoCarlo$方法的收敛性和误差进行定量的描述。

2.收敛性的描述

​ 若随机变量$X_i(i=1,2,...N)$是独立同分布的随机变量，其期望值$Ex<\infty$,由辛钦大数定律计算可以得到如下关系,$\forall ϵ>0:$
l i m N → ∞ P { ∣ 1 N ∑ i = 1 N X i − E x ∣ < ϵ } = 1 lim_{N\rightarrow\infty}P\{|\frac{1}{N}\sum_{i=1}^{N}X_i-Ex|<\epsilon\}=1 limN→∞​P{∣N1​i=1∑N​Xi​−Ex∣<ϵ}=1
​ 如果我们设:
$${\overline{X}}_N=\frac{1}{N}\sum _{i=1}^N{X}_i$$
​ 那么当$N\to \infty$时，${\overline{X}}_N$依概率$1$收敛到$Ex$。这就是其收敛性的数学基础。

3.误差的描述与控制

​ 若随机变量$X_i(i=1,2,...N)$是独立同分布的随机变量，且具有误差${\sigma }^2({\sigma }^2\ne 0)$，且：
$${\sigma }^2={\int }_{-\infty }^{\infty }(x-Ex{)}^{2}f(x)dx<\infty$$
​ 其中：$Ex$为数学期望，$f(x)$为$X_i$的概率密度函数。

​ 如果我们设样本个数为$N$，且${\overline{X}}_N=\frac{1}{N}{\sum }_{i=1}^N{X}_i$，由中心极限定理:
$$li{m}_{n\to \infty }\sum _{i=1}^N{X}_i\sim N(NEx,N{\sigma }^2)$$
​ 得到如下的数学关系:
l i m N → ∞ P { ∣ ∑ i = 1 N X i − N E x N σ 2 ∣ < λ a } = 1 2 π ∫ − λ a λ a e − t 2 2 d t lim_{N \rightarrow \infty}P\{|\frac{\sum_{i=1}^NX_i-NEx}{\sqrt{N\sigma^2}}|<\lambda_a\}=\frac{1}{\sqrt{2\pi}}\int_{-\lambda_a}^{\lambda_a}e^{-\frac{t^2}{2}}dt limN→∞​P{∣Nσ2 ​∑i=1N​Xi​−NEx​∣<λa​}=2π ​1​∫−λa​λa​​e−2t2​dt
​ 化简一下：
l i m N → ∞ P { ∣ X ‾ N − E x ∣ < λ a σ N } = 2 2 π ∫ 0 λ a e − t 2 2 d t = 1 − α lim_{N\rightarrow\infty}P\{|\overline{X}_N-Ex|<\frac{\lambda_a \sigma}{\sqrt{N}}\}=\frac{2}{\sqrt{2\pi}}\int_0^{\lambda_a}e^{-\frac{t^2}{2}}dt=1-\alpha limN→∞​P{∣XN​−Ex∣<N ​λa​σ​}=2π ​2​∫0λa​​e−2t2​dt=1−α
​ 其中，$\alpha$为置信度，当$\alpha =0.5,0.05,0.003$时，${\lambda }_a$分别对应于$0.6745,1.96,3$。

​ 因此我们可以给出误差的定量描述：$|{\overline{X}}_N-Ex|<\frac{{\lambda }_a\sigma }{\sqrt{N}}$以$1-\alpha$的概率成立。而且误差的收敛速度的阶数为$O(N^{-\frac{1}{2}})$。因此到后面随着$N$的增加，误差收敛的速度会越来越慢。

​ 如果我们做一下定义$ϵ$为误差：
$$ϵ=\frac{{\lambda }_a\sigma }{\sqrt{N}}$$
​ 那么要想控制$N$使得误差以$1-\alpha$的概率小于$ϵ$，那么:
$$N>(\frac{{\lambda }_a\sigma }{ϵ}{)}^2$$
​ 或者控制$\sigma$使得误差以$1-\alpha$的概率小于$ϵ$，那么:
$$\sigma <\frac{\sqrt{N}ϵ}{{\lambda }_a}$$
​ 如果要控制$N$的话，那我们就得计算$\sigma$的值，然而$\sigma$的值不太好计算。我们就用下面的式子做$\sigma$的无偏估计：
$$\hat{\sigma }=\sqrt{\frac{1}{N}\sum _{i=1}^N{X}_i^2-(\frac{1}{N}\sum _{i=1}^N{X}_i{)}^2}$$

4.减少误差的技巧

​ 在给定置信度$\alpha$后，$ϵ$由$\sigma$和$N$共同决定，那么增大$N$和减小$\sigma$都可以减小误差$ϵ$。但是，$ϵ$随和$N$却是呈现平方反比的关系。因此单纯的增大$N$效果并不明显，还得想办法减小$\sigma$的值来减小误差。

关于减少$\sigma$的各种技巧，本篇博文就不再赘述。

5.代码的实现

要求计算${\iiint }_{\Omega }{x}^{2}dV$,其中 Ω : { ( x , y , z ) ∣ x 2 + y 2 ≤ z ≤ 1 } \Omega:\{(x,y,z)|\sqrt{x^2+y^2}\leq z\leq 1\} Ω:{(x,y,z)∣x2+y2 ​≤z≤1}。并且给出误差的动态描述。规定$\alpha =0.05$。

#计算Mentor Carlo积分
import random
import matplotlib.pyplot as plt
import math
def Cal_cul(x,y,z):
    if (z<1)&(z>=math.sqrt(x**2+y**2)):
        return x**2
    else:
        return 0

def lambdaTable(alpha):
    if alpha==0.5:
        return 0.6745
    elif alpha == 0.05:
        return 1.96
    else:
        return 3

def int3Cul(a,b,c,Numbers,alpha):
    V = 2*a*2*b*c
    sum1 = 0
    sum2 = 0
    lambda_a = lambdaTable(alpha)
    for i in range(Numbers):
        x = -a+2*a*random.random()
        y = -b+2*b*random.random()
        z = c*random.random()
        sum1 += Cal_cul(x,y,z)
        sum2 += (Cal_cul(x,y,z))**2
    int3Value = V*1/Numbers*sum1
    sigma = math.sqrt(1/Numbers*sum2-(1/Numbers*sum1)**2)
    error = lambda_a*sigma/math.sqrt(Numbers)
    print('在试验次数为{:}的情况下,计算的积分值是{:.3f},有{:.2f}的概率认为误差小于{:.5f}'.format(Numbers,int3Value,1-alpha,error))

if __name__== '__main__':
    int3Cul(1,1,1,10000,0.05)

输出的结果是:

在试验次数为10000的情况下,计算的积分值是0.159,有0.95的概率认为误差小于0.00223

Process finished with exit code 0

实际上我们可以知道上述三重积分理论值是$\frac{\pi }{20}$。而计算的积分值是$0.159$。它们的误差值:
$$|\frac{\pi }{20}-0.159|\dot{=}0.00192<0.00223$$
说明我们的误差预测大概率是正确的。