通俗易懂的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 如果Xi​是独立的随机变量，且EXi​是相应的数学期望，∀ϵ>0,(i=1,2,...N),则有：limn→∞​P{∣n1​∑i=1n​Xi​−n1​∑i=1n​EXi​∣<ϵ}=1

  该定律指出了当样本数量足够大的时候，样本的算术平均值可以无限逼近数学期望

• Mentor Carlo期望法的计算

  主要是用随机数序列做随机模拟的数值方法。精度比投点法要高的很多。现在我们直接从多重积分来进行考虑

  考虑一个$k$维积分：
  $$I={\int }_{\Omega }...\int f(x_1,x_2,x_3,...,x_k)d{x}_{1}d{x}_{2}d{x}_3...d{x}_k$$
  $\Omega$是一个$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)$必须满足以下条件，而且选取的随机变量$(x_1,x_2,...,x_k)$必须服从该概率密度$g(x_1,x_2,...,x_k)$的随机分布：

  • $g(x_1,x_2,...x_k)\ge 0$

  • $(x_1,x_2,...x_k)\in \Omega$

  • ${\int }_{\Omega }...\int g(x_1,x_2,x_3,...,x_k)d{x}_{1}d{x}_{2}d{x}_3...d{x}_k=1$

  但是，为了使得我们的随机变量是在一个规范的$k$维矩形体内产生的，我们有必要将积分$I$进行一下变换，变换成规则的区域$V$上的积分。变换方法如下：
  $$I={\int }_V...\int f^{\ast }(x_1,x_2,...,x_k)d{x}_{1}d{x}_2...d{x}_k$$
  其中：
  $$f^{\ast }(x_1,x_2,...,x_k)=\{\begin{array}{ll}f(x_1,x_2,...,x_k)& \text{if (x1,x2,..xk)∈Ω}\\ 0& \text{(x1,x2,..xk)∈/Ω}\end{array}$$

  在$V$上的随机变量的密度函数为$g^{\ast }(x_1,x_2,...x_k)$，满足以下条件$$g^{\ast }(x_1,x_2,...x_k)>0,(x_1,x_2,...,x_k)\in V,{\int }_V...\int g^{\ast }(x_1,x_2,...x_k)d{x}_{1}d{x}_2...d{x}_k=1$$
  $$\begin{gathered} I={\int }_V...\int f^{\ast }(x_1,x_2,...x_k)d{x}_{1}d{x}_2...d{x}_k \\ ={\int }_V...\int \frac{f^{\ast }(x_1,x_2,...x_k)}{g^{\ast }(x_1,x_2,...x_k)}{g}^{\ast }(x_1,x_2,...x_k)d{x}_{1}d{x}_2...d{x}_k \\ =E\frac{f^{\ast }(x_1,x_2,...x_k)}{g^{\ast }(x_1,x_2,...x_k)} \end{gathered}$$
  由辛钦大数定律有：
  ∀ ϵ > 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{∣N1​i=1∑N​g∗(xi1​,xi2​,...xik​)f∗(xi1​,xi2​,...xik​)​−E{g∗(x1​,x2​,...xk​)f∗(x1​,x2​,...xk​)​}∣<ϵ}=1
  其中：${\xi }_i=(x_{i1},x_{i2},...,x_{ik})$是在$V$中以概率密度函数为$g^{\ast }(x_1,x_2,...x_k)$的随机一次取样，总共要在$V$中要取样$N$次。得到：
  $$I=\frac{1}{N}\sum _{i=1}^N\frac{f^{\ast }(x_{i1},x_{i2},...x_{ik})}{g^{\ast }(x_{i1},x_{i2},...x_{ik})}$$

• Python的代码实现

  例：计算$\int {\int }_D{x}^2+y^{2}d\sigma$。其中: σ = { ( 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))
  

  理论值是：$\frac{\pi }{8}$，计算值是：$0.3912$。非常接近。