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

通俗易懂的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 XiEXiϵ>0,(i=1,2,...N),limnP{n1i=1nXin1i=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)aixibi,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,limNP{N1i=1Ng(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=1Ng(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,0x1,0y1}

    #用蒙特卡洛计算二重积分
    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。非常接近。

©️2020 CSDN 皮肤主题: 1024 设计师:上身试试 返回首页