Search code examples
pythonalgorithmmontecarlo

Montecarlo integration in D dimension in Python

I'm tring to solve a D dimensional integral by Monte Carlo Integration:

enter image description here

The idea is to generate N point and calculate the aria below te curve as:

enter image description here

In order to do this i implemented this Python code:

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 output are:

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 .

Comparing the exact results of integral with the results of Monte Carlo integration, we can see that the Monte Carlo integration failed.

Where is the error?

Thanks in advance.

  • John Snowden
Solution

  • Well, why do you need this "below curve" crap?

    You're integrating over hypercube, just compute mean value of the function and be done.

    E.g., in 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)

    will print something like

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

    and for 6D case

    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)

    I've got

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