Monte Carlo Integration using Python

 Monte Carlo Integration

Monte Carlo integration is a  powerful method of computing the value of complex integrals using probabilistic techniques. This technique uses random numbers to compute the definite integral of a function. Here we are going to use the python programs written in the previous post to generate pseudorandom numbers and approximate the value of the definite integral of a function.

Consider a function to be integrated, as shown below:

That is we need to evaluate the definite integral $\int_{a}^{b} f(x)\mathrm{d}x$. The integral is just the area under the curve. The width of the interval $(b- a)$ times the average value of the function is also the value of the integral, that is,

$\int_{a}^{b}f(x)\mathrm{d}x = (b - a)f_{average} = (b-a)\langle f \rangle$

So if we had some independent way of calculating the average value of the integrand, then we could evaluate the integral. This is where random numbers come in.

Imagine that we had a list of random numbers, $x_j$ uniformly distributed between $a$ and $b$. To calculate the function average, we simply evaluate $f(x_j)$ at each of the randomly selected points, add all of them together and divide the sum by the number of points.

$$\langle f \rangle _N = \frac{1}{N}\sum_{i=1}^{N}f(x_i)$$

As the number of points used in calculating the average increases, $\langle f \rangle_N$ approaches the true value, $\langle f \rangle$. Therefore, as a numerical approximation, we can write

$$\int_{a}^{b}f(x)\mathrm{d}x = \frac{(b-a)}{N}\sum_{i=1}^{N}f(x_i)$$

Alternatively, we can look at the Monte Carlo integration method in the following way:

• Find some value $M$ such that $f(x) < M$ over the integral $[a,b]$.
• Select random numbers $x$ and $y$ from the intervals $[a,b]$ and $[0,M]$ respectively.
• Determine if $y > f(x)$ or $y \leq f(x)$.
•  Repeat the process $N$ times, keeping track of the number of times $y \leq f(x)$ or under the curve (successes); call the total number of successes $S$.
  

The estimated probability of success is

$$\frac{S}{N} = \frac{Area \ under \ curve}{Total \ area \ inside \ rectangle} = \frac{\int_{a}^{b}f{x}\mathrm{d}x}{M(b-1)}$$

Thus the value of the integral can be calculated from the formula

$$\int_{a}^{b}f(x)\mathrm{d}x = M(b-a)\frac{S}{N}$$

The python program that implements the Monte Carlo method of integration is given below. Here the random numbers used in the integration is generated using Lehmer's method mention in the previous post.

Program

from sympy import *

def pseudo(a,m,x0):

    list1 = []

    for i in range(100):

        x0 = (a*x0)%m

        list1.append(x0)

    return(list1)

def integrate(f,var):

    f = sympify(f)

    var = Symbol(var)

    rand_values1 = pseudo(61,101,18)

    rand_values2 = pseudo(63,107,18)

    f_values1 = []

    f_values2 = []

    x = []

    y = []

    for i in rand_values1:

        x.append(a + (((b-1)/101)*i))

    i = a

    #to find maximum M

    while i <= b:

        f_values1.append(f.subs({var:i}))

        i += 0.01

    #to check whether y lies below or above the curve

    for i in x:

        f_values2.append(f.subs({var:i}))

    M = int(max(f_values1)) + 1

    for i in rand_values2:

        y.append((M/107)*(i))

    count = 0

    for i in range(len(y)):

        if y[i] <= f_values2[i]:

            count += 1

    return (b-a)*M*(count/100)

f = input("Enter the function: ")

var = input("Enter the variable: ")

a = int(input("Enter the lower limit: "))

b = int(input("Enter the upper limit: "))

print("The value of the definite integral is", integrate(f,var))

The value of some definite integrals obtained using the above program and their actual values are given below in the table: