Python中D维的Montecarlo集成

Question

我试着用蒙特卡罗积分法求解一个D维积分：

​

​

其思想是生成N个点并计算te曲线下面的aria，如下所示：

​

​

为了做到这一点，我实现了以下Python代码：

import numpy as np
from sympy import symbols, integrate     



def f(x,D):                             
  return D*(x**2)



for i in range(1, 9):                    

  x = symbols('x')                      

  print("The exact mathematical value of the integral with D egual", i, "is:", integrate(f(x,i),(x, 0,1)).evalf(2), "\n") 

print("************************************************************************* \n")



N = 10**4

for j in range(1,9):

  ans = 0

  n_tot = N

  n_below_curve = 0


  for i in range(N):

    x0=np.random.uniform(0,1)
    y0=np.random.uniform(0,1)


    if (f(x0,j) <= y0):

      n_below_curve += 1

  ans = ( n_below_curve / n_tot ) * (1*1)



  print("The result of integral with D egual to", j, "is:", ans, ".\n")

产出如下：

The exact mathematical value of the integral with D egual 1 is: 0.33 

The exact mathematical value of the integral with D egual 2 is: 0.67 

The exact mathematical value of the integral with D egual 3 is: 1.0 

The exact mathematical value of the integral with D egual 4 is: 1.3 

The exact mathematical value of the integral with D egual 5 is: 1.7 

The exact mathematical value of the integral with D egual 6 is: 2.0 

The exact mathematical value of the integral with D egual 7 is: 2.3 

The exact mathematical value of the integral with D egual 8 is: 2.7 

************************************************************************* 

The result of integral with D egual to 1 is: 0.6635 .

The result of integral with D egual to 2 is: 0.4681 .

The result of integral with D egual to 3 is: 0.3823 .

The result of integral with D egual to 4 is: 0.3321 .

The result of integral with D egual to 5 is: 0.2978 .

The result of integral with D egual to 6 is: 0.269 .

The result of integral with D egual to 7 is: 0.252 .

The result of integral with D egual to 8 is: 0.2372 .

将积分的精确结果与蒙特卡罗积分的结果进行比较，可以看出蒙特卡罗积分失败了。

错误在哪里？

提前谢谢。

• 约翰·斯诺登

Answer 1

你为什么要这个“曲线下”的废话？

在超立方体上进行积分，只需计算函数的平均值就可以了。

例如，在3D

import numpy as np
from scipy import integrate

rng = np.random.default_rng()

D = 3

N = 100000

I = 0.0 # accumulator
for k in range(0, N):
    pt = rng.random(D) # single point sampled
    I += np.sum(pt*pt) # x0^2 + x1^2 + ...

print(I/N) # mean value

def func(x0, x1, x2):
    return x0*x0 + x1*x1 + x2*x2

R = integrate.nquad(func, ((0,1), (0,1), (0,1)), full_output=True)
print(R)

会打印出这样的东西

1.0010147193589627
(1.0, 2.5808878251226036e-14, {'neval': 9261})

对于6D案件

def func(x0, x1, x2, x3, x4, x5):
    return x0*x0 + x1*x1 + x2*x2 + x3*x3 + x4*x4 + x5*x5

R = integrate.nquad(func, ((0,1), (0,1), (0,1), (0,1), (0,1), (0,1)), full_output=True)

我有

1.9997059362936607
(2.0, 5.89710805049393e-14, {'neval': 85766121})