Monte Carlo integration in Python
Monte Carlo Integration is a process of solving integrals having numerous values to integrate upon. The Monte Carlo process uses the theory of large numbers and random sampling to approximate values that are very close to the actual solution of the integral. It works on the average of a function denoted by <f(x)>. Then we can expand <f(x)> as the summation of the values divided by the number of points in the integration and solve the Left-hand side of the equation to approximate the value of the integration at the right-hand side. The derivation is as follows.
Where N = no. of terms used for approximation of the values.
Now we would first compute the integral using the Monte Carlo method numerically and then finally we would visualize the result using a histogram by using the python library matplotlib.
The modules that would be used are:
- Scipy: Used to get random values between the limits of the integral. Can be installed using:
pip install scipy # for windows or pip3 install scipy # for linux and macos
- Numpy: Used to form arrays and store different values. Can be installed using:
pip install numpy # for windows or pip3 install numpy # for linux and macos
- Matplotlib: Used to visualize the histogram and the result of the Monte Carlo integration.
pip install matplotlib # for windows or pip3 install matplotlib # for linux and macos
The integral we are going to evaluate is:
Once we are done installing the modules, we can now start writing code by importing the required modules, scipy for creating random values in a range and NumPy module to create static arrays as default lists in python are relatively slow because of dynamic memory allocation. Then we defined the limits of integration which are 0 and pi (we used np.pi to get the value of pi. Then w create a static numpy array of zeros of size N using the np. zeros(N). Now we iterate over each index of array N and fill it with random values between a and b. Then we create a variable to store sum of the functions of different values of the integral variable. Now we create a function to calculate the sin of a particular value of x. Then we iterate and add up all the values returned by the function of x for different values of x to the variable integral. Then we use the formula derived above to get the results. Finally, we print the results on the final line.
The value calculated by monte carlo integration is 2.0256756150767035.
The value obtained is very close to the actual answer of the integral which is 2.0.
Now if we want to visualize the integration using a histogram, we can do so by using the matplotlib library. Again we import the modules, define the limits of integration and write the sin function for calculating the sin value for a particular value of x. Next, we take an array that has variables representing every beam of the histogram. Then we iterate through N values and repeat the same process of creating a zeros array, filling it with random x values, creating an integral variable adding up all the function values, and getting the answer N times, each answer representing a beam of the histogram. The code is as follows:
Here as we can see, the most probable results according to this graph would be 2.023 or 2.024 which is again quite close to the actual solution of this integral that is 2.0.
The integral we are going to evaluate is:
The implementation would be the same as the previous question, only the limits of integration defined by a and b would change and also the f(x) function that returns the corresponding function value of a specific value would now return x^2 instead of sin (x). The code for this would be:
The value calculated by monte carlo integration is 0.33024026610116575.
The value obtained is very close to the actual answer of the integral which is 0.333.
Now if we want to visualize the integration using a histogram again, we can do so by using the matplotlib library as we did for the previous expect in this case the f(x) functions return x^2 instead of sin(x), and the limits of integration change from 0 to 1. The code is as follows:
Here as we can see, the most probable results according to this graph would be 0.33 which is almost equal (equal in this case, but it generally doesn’t occur) to the actual solution of this integral that is 0.333.