通俗易懂的Monte Carlo积分方法(二)
Monte Carlo积分的计算(期望法)
-
Monte Carlo算法的期望法计算的数学基础:
辛钦大数定律:
如 果 X i 是 独 立 的 随 机 变 量 , 且 E X i 是 相 应 的 数 学 期 望 , ∀ ϵ > 0 , ( i = 1 , 2 , . . . N ) , 则 有 : l i m n → ∞ P { ∣ 1 n ∑ i = 1 n X i − 1 n ∑ i = 1 n E X i ∣ < ϵ } = 1 如果X_i 是独立的随机变量,且EX_i是相应的数学期望,\\\forall\epsilon>0,(i=1,2,...N),则有:\\lim_{n\rightarrow\infty}P\{|\frac{1}{n}\sum_{i=1}^nX_i-\frac{1}{n}\sum_{i=1}^nEX_i|<\epsilon\}=1 如果Xi是独立的随机变量,且EXi是相应的数学期望,∀ϵ>0,(i=1,2,...N),则有:limn→∞P{∣n1∑i=1nXi−n1∑i=1nEXi∣<ϵ}=1该定律指出了当样本数量足够大的时候,样本的算术平均值可以无限逼近数学期望
-
Mentor Carlo期望法的计算
主要是用随机数序列做随机模拟的数值方法。精度比投点法要高的很多。现在我们直接从多重积分来进行考虑
考虑一个 k k k维积分:
I = ∫ Ω . . . ∫ f ( x 1 , x 2 , x 3 , . . . , x k ) d x 1 d x 2 d x 3 . . . d x k I=\int_{\Omega}...\int f(x_1,x_2,x_3,...,x_k)dx_1dx_2dx_3...dx_k I=∫Ω...∫f(x1,x2,x3,...,xk)dx1dx2dx3...dxk
Ω \Omega Ω是一个 k k k维积分区域。考虑选取 Ω \Omega Ω上的一个概率密度: g ( x 1 , x 2 , . . . , x k ) , Ω ⊆ V ( V : { ( x 1 , x 2 , . . x k ) ∣ a i ≤ x i ≤ b i , i = 1 , 2 , . . . k } ) g(x_1,x_2,...,x_k),\Omega\subseteq V\\(V:\{(x_1,x_2,..x_k)|a_i\leq x_i\leq b_i,i=1,2,...k\}) g(x1,x2,...,xk),Ω⊆V(V:{(x1,x2,..xk)∣ai≤xi≤bi,i=1,2,...k})其中 g ( x 1 , x 2 , . . . , x k ) g(x_1,x_2,...,x_k) g(x1,x2,...,xk)必须满足以下条件,而且选取的随机变量 ( x 1 , x 2 , . . . , x k ) (x_1,x_2,...,x_k) (x1,x2,...,xk)必须服从该概率密度 g ( x 1 , x 2 , . . . , x k ) g(x_1,x_2,...,x_k) g(x1,x2,...,xk)的随机分布:
-
g ( x 1 , x 2 , . . . x k ) ≥ 0 g(x_1,x_2,...x_k)\geq 0 g(x1,x2,...xk)≥0
-
( x 1 , x 2 , . . . x k ) ∈ Ω (x_1,x_2,...x_k)\in\Omega (x1,x2,...xk)∈Ω
-
∫ Ω . . . ∫ g ( x 1 , x 2 , x 3 , . . . , x k ) d x 1 d x 2 d x 3 . . . d x k = 1 \int_{\Omega}...\int g(x_1,x_2,x_3,...,x_k)dx_1dx_2dx_3...dx_k=1 ∫Ω...∫g(x1,x2,x3,...,xk)dx1dx2dx3...dxk=1
但是,为了使得我们的随机变量是在一个规范的 k k k维矩形体内产生的,我们有必要将积分 I I I进行一下变换,变换成规则的区域 V V V上的积分。变换方法如下:
I = ∫ V . . . ∫ f ∗ ( x 1 , x 2 , . . . , x k ) d x 1 d x 2 . . . d x k I=\int_{V}...\int f^*(x_1,x_2,...,x_k)dx_1dx_2...dx_k I=∫V...∫f∗(x1,x2,...,xk)dx1dx2...dxk
其中:
f ∗ ( x 1 , x 2 , . . . , x k ) = { f ( x 1 , x 2 , . . . , x k ) if ( x 1 , x 2 , . . x k ) ∈ Ω 0 ( x 1 , x 2 , . . x k ) ∉ Ω f^*(x_1,x_2,...,x_k)=\begin{cases} f(x_1,x_2,...,x_k) & \text{if $(x_1,x_2,..x_k)\in \Omega$}\\ 0 &\text{$(x_1,x_2,..x_k)\notin\Omega$} \end{cases} f∗(x1,x2,...,xk)={f(x1,x2,...,xk)0if (x1,x2,..xk)∈Ω(x1,x2,..xk)∈/Ω在 V V V上的随机变量的密度函数为 g ∗ ( x 1 , x 2 , . . . x k ) g^*(x_1,x_2,...x_k) g∗(x1,x2,...xk),满足以下条件 g ∗ ( x 1 , x 2 , . . . x k ) > 0 , ( x 1 , x 2 , . . . , x k ) ∈ V , ∫ V . . . ∫ g ∗ ( x 1 , x 2 , . . . x k ) d x 1 d x 2 . . . d x k = 1 g^*(x_1,x_2,...x_k)>0,(x_1,x_2,...,x_k)\in V,\int_V...\int g^*(x_1,x_2,...x_k)dx_1dx_2...dx_k =1 g∗(x1,x2,...xk)>0,(x1,x2,...,xk)∈V,∫V...∫g∗(x1,x2,...xk)dx1dx2...dxk=1
I = ∫ V . . . ∫ f ∗ ( x 1 , x 2 , . . . x k ) d x 1 d x 2 . . . d x k = ∫ V . . . ∫ f ∗ ( x 1 , x 2 , . . . x k ) g ∗ ( x 1 , x 2 , . . . x k ) g ∗ ( x 1 , x 2 , . . . x k ) d x 1 d x 2 . . . d x k = E f ∗ ( x 1 , x 2 , . . . x k ) g ∗ ( x 1 , x 2 , . . . x k ) I=\int_V...\int f^*(x_1,x_2,...x_k)dx_1dx_2...dx_k\\ =\int_V...\int \frac{f^*(x_1,x_2,...x_k)}{g^*(x_1,x_2,...x_k)}g^*(x_1,x_2,...x_k)dx_1dx_2...dx_k\\ =E\frac{f^*(x_1,x_2,...x_k)}{g^*(x_1,x_2,...x_k)} I=∫V...∫f∗(x1,x2,...xk)dx1dx2...dxk=∫V...∫g∗(x1,x2,...xk)f∗(x1,x2,...xk)g∗(x1,x2,...xk)dx1dx2...dxk=Eg∗(x1,x2,...xk)f∗(x1,x2,...xk)
由辛钦大数定律有:
∀ ϵ > 0 , l i m N → ∞ P { ∣ 1 N ∑ i = 1 N f ∗ ( x i 1 , x i 2 , . . . x i k ) g ∗ ( x i 1 , x i 2 , . . . x i k ) − E { f ∗ ( x 1 , x 2 , . . . x k ) g ∗ ( x 1 , x 2 , . . . x k ) } ∣ < ϵ } = 1 \forall \epsilon >0,lim_{N\rightarrow\infty}P\{|\frac{1}{N}\sum_{i=1}^{N}\frac{f^*(x_{i1},x_{i2},...x_{ik})}{g^*(x_{i1},x_{i2},...x_{ik})}-E\{\frac{f^*(x_1,x_2,...x_k)}{g^*(x_1,x_2,...x_k)}\}|<\epsilon\}=1 ∀ϵ>0,limN→∞P{∣N1i=1∑Ng∗(xi1,xi2,...xik)f∗(xi1,xi2,...xik)−E{g∗(x1,x2,...xk)f∗(x1,x2,...xk)}∣<ϵ}=1
其中: ξ i = ( x i 1 , x i 2 , . . . , x i k ) \xi_i=(x_{i1},x_{i2},...,x_{ik}) ξi=(xi1,xi2,...,xik)是在 V V V中以概率密度函数为 g ∗ ( x 1 , x 2 , . . . x k ) g^*(x_1,x_2,...x_k) g∗(x1,x2,...xk)的随机一次取样,总共要在 V V V中要取样 N N N次。得到:
I = 1 N ∑ i = 1 N f ∗ ( x i 1 , x i 2 , . . . x i k ) g ∗ ( x i 1 , x i 2 , . . . x i k ) I=\frac{1}{N}\sum_{i=1}^{N}\frac{f^*(x_{i1},x_{i2},...x_{ik})}{g^*(x_{i1},x_{i2},...x_{ik})} I=N1i=1∑Ng∗(xi1,xi2,...xik)f∗(xi1,xi2,...xik) -
-
Python的代码实现
例:计算 ∫ ∫ D x 2 + y 2 d σ \int\int_{D}x^2+y^2d\sigma ∫∫Dx2+y2dσ。其中: σ = { ( x , y ) ∣ x 2 + y 2 = 1 , 0 ≤ x ≤ 1 , 0 ≤ y ≤ 1 } \\ \sigma=\{(x,y)|x^2+y^2=1,0 \leq x \leq1, 0\leq y\leq 1\} σ={(x,y)∣x2+y2=1,0≤x≤1,0≤y≤1}
#用蒙特卡洛计算二重积分 def Cal_cul(x,y): if (x**2+y**2)<1: value = x**2+y**2 else: value = 0 return value import random,cmath random.seed(10) numbers = 10000 int_value = 0 sum = 0 for i in range(numbers): x_value = random.random() y_value = random.random() sum += Cal_cul(x_value,y_value) int_value = sum/numbers print("最后的值是:{:.4f} ".format(int_value))
理论值是: π 8 \frac{\pi}{8} 8π,计算值是: 0.3912 0.3912 0.3912。非常接近。