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

Monte Carlo 方法计算的理论基础

1.理论目标
2.收敛性的描述
3.误差的描述与控制
4.减少误差的技巧
5.代码的实现

1.理论目标

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

2.收敛性的描述

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

3.误差的描述与控制

若随机变量Xi(i=1,2,...N)X_i(i=1,2,...N)Xi(i=1,2,...N)是独立同分布的随机变量，且具有误差σ2(σ2≠0)\sigma^2(\sigma^2\ne0)σ2(σ2=0)，且：
σ2=∫−∞∞(x−Ex)2f(x)dx<∞\sigma^2=\int_{-\infty}^{\infty}(x-Ex)^2f(x)dx<\infty σ2=∫−∞∞(x−Ex)2f(x)dx<∞
其中：ExExEx为数学期望，f(x)f(x)f(x)为XiX_iXi的概率密度函数。

如果我们设样本个数为NNN，且X‾N=1N∑i=1NXi\overline{X}_N=\frac{1}{N}\sum_{i=1}^NX_iXN=N1∑i=1NXi，由中心极限定理:
limn→∞∑i=1NXi∼N(NEx,Nσ2)lim_{n \rightarrow \infty}\sum_{i=1}^NX_i \sim N(NEx,N\sigma^2) limn→∞i=1∑NXi∼N(NEx,Nσ2)
得到如下的数学关系:
limN→∞P{∣∑i=1NXi−NExNσ2∣<λa}=12π∫−λaλae−t22dtlim_{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=1NXi−NEx∣<λa}=2π1∫−λaλae−2t2dt
化简一下：
limN→∞P{∣X‾N−Ex∣<λaσN}=22π∫0λae−t22dt=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λae−2t2dt=1−α
其中，α\alphaα为置信度，当α=0.5,0.05,0.003\alpha=0.5,0.05,0.003α=0.5,0.05,0.003时，λa\lambda_aλa分别对应于0.6745,1.96,30.6745,1.96,30.6745,1.96,3。

因此我们可以给出误差的定量描述：∣X‾N−Ex∣<λaσN|\overline{X}_N-Ex|<\frac{\lambda_a\sigma}{\sqrt{N}}∣XN−Ex∣<Nλaσ以1−α1-\alpha1−α的概率成立。而且误差的收敛速度的阶数为O(N−12)O(N^{-\frac{1}{2}})O(N−21)。因此到后面随着NNN的增加，误差收敛的速度会越来越慢。

如果我们做一下定义ϵ\epsilonϵ为误差：
ϵ=λaσN\epsilon = \frac{\lambda_a\sigma}{\sqrt{N}} ϵ=Nλaσ
那么要想控制NNN使得误差以1−α1-\alpha1−α的概率小于ϵ\epsilonϵ，那么:
N>(λaσϵ)2N>(\frac{\lambda_a\sigma}{\epsilon})^2 N>(ϵλaσ)2
或者控制σ\sigmaσ使得误差以1−α1-\alpha1−α的概率小于ϵ\epsilonϵ，那么:
σ<Nϵλa\sigma<\frac{\sqrt{N}\epsilon}{\lambda_a} σ<λaNϵ
如果要控制NNN的话，那我们就得计算σ\sigmaσ的值，然而σ\sigmaσ的值不太好计算。我们就用下面的式子做σ\sigmaσ的无偏估计：
σ^=1N∑i=1NXi2−(1N∑i=1NXi)2\hat{\sigma} = \sqrt{\frac{1}{N}\sum^N_{i=1}X_i^2-(\frac{1}{N}\sum_{i=1}^NX_i)^2} σ^=N1i=1∑NXi2−(N1i=1∑NXi)2

4.减少误差的技巧

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

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

5.代码的实现

要求计算∭Ωx2dV\iiint_{\Omega}x^2dV∭Ωx2dV,其中Ω:{(x,y,z)∣x2+y2≤z≤1}\Omega:\{(x,y,z)|\sqrt{x^2+y^2}\leq z\leq 1\}Ω:{(x,y,z)∣x2+y2≤z≤1}。并且给出误差的动态描述。规定α=0.05\alpha = 0.05α=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**2else:return 0def lambdaTable(alpha):if alpha==0.5:return 0.6745elif alpha == 0.05:return 1.96else:return 3def int3Cul(a,b,c,Numbers,alpha):V = 2*a*2*b*csum1 = 0sum2 = 0lambda_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))**2int3Value = V*1/Numbers*sum1sigma = 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.00223Process finished with exit code 0

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