Python 中 D 维的蒙特卡洛积分

Question

我正在尝试通过 Monte Carlo 求解 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 .

将积分的确切结果与Monte Carlo积分的结果进行比较，可以看出Monte Carlo积分失败。

哪里出错了？

提前致谢。

约翰·斯诺登

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})

Python 中 D 维的蒙特卡洛积分

Montecarlo integration in D dimension in Python

python

algorithm

montecarlo