# -*- coding: utf-8 -*- """ Created on Fri Oct 2 12:48:44 2020 @author: jbobowsk """ # This tutorial will numerically evaluate a definite integral using the # Monte Carlo "f-average" method. import numpy as np import matplotlib.pyplot as plt # In all Monte Carlo simulations it is necessary to generate random or # pseudo-random numbers. The Python command 'np.random.uniform()' # from the NumPy module generates random numbers uniformly distributed # between zero and one. print(np.random.uniform()) # To confirm that the numbers generated are uniformly distribted on # the interval [0, 1], we can generate a long list of these numbers and # plot a histogram. 'np.random.uniform(a, b, n)' returns n random number # uniformly distributed between a and b. nos = np.random.uniform(0, 1, int(1e4)) # Here's the histogram. nbins = 50 plt.hist(nos, nbins, color = 'c', edgecolor = 'k') plt.xlabel('x'); plt.ylabel('counts'); # This tutorial will first attempt to evaluate an integral for which the # solution is easily obtained. This approach has been taken purposely so # that we can confirm that our numerical techniques are reliable. The # function that we will integrate is a simply polynomial. Below, the # function is plotted on the interval x = [0, 1] and the exact value of the # integral is evaluated over the same interval. We are going to use this # polynomial over and over again in this tutorial, so we define it as a # function. fcn = lambda x : (-65536*x**8 + 262144*x**7 - 409600*x**6 + 311296*x**5 -\ 114688*x**4 + 16384*x**3)/27 xx = np.arange(0, 1, 0.001) plt.figure() plt.plot(xx, fcn(xx), 'k-', linewidth = 1.2) plt.xlabel('x') plt.ylabel('f(x)') # Here's how we can use the SymPy module to do symbolic integration in # Python. This is an indefinite integral. import sympy as sym x = sym.Symbol('x') fint = sym.integrate(fcn(x), x) print(fint) # Here is the definite integral of our function over 0 < x < 1. fint = sym.integrate(fcn(x), (x, 0, 1)) print(fint) print('The exact integral of f(x) on 0 < x < 1 is', sym.N(fint)) # The f-average method requires that we repeatedly generate x_i random # numbers (n times) uniformly distributed between the integration interval # (0 to 1 in this example). For each x_i we calculate the corresponding # f(x_i) and sum all of these to find f_total. f_average is obtained # from f_average = f_total/n. An estimate of the integral is given by # the integration interval (b-a)*f_average. In the plot below, we show # each of the f(x_i) and f_average which is our estimate of the integral. # The square area below the red line is approximately equal to the area # beneath the blue curve which is the integral that we're trying to evaluate. n = int(1e4) ftot = 0 xList = [] fList = [] for i in range(n): xx = np.random.uniform() xList = xList + [xx] fList = fList + [fcn(xx)] ftot = ftot + fcn(xx) favg = ftot/n print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', favg) # Plot the function as a line and the add the hits and misses as scatter # plots. plt.figure() plt.plot(xList, fList, 'bo', fillstyle = 'none') plt.plot(xList, favg*np.ones(n), 'r-', linewidth = 3) plt.xlabel('x') plt.ylabel('f(x)') # If you don't need to keep track of the actual coordinates of each # (x, f(x)), then you can get away with a more compact loop keeps a running # total of all the f(x) values. n = int(1e4) ftot = 0 for i in range(n): xx = np.random.uniform() ftot = ftot + fcn(xx) favg = ftot/n print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', favg) # Now we will numberically approximate the integral using n = 10,000 one # thousand times and plot the resulting distribution of our determination of # the integral. The distribution is expected to be Gaussian. ## To keep track of how long this block of codes takes to complete, I'll use the ## "datetime.now()" function to display the time before and after the loop. from datetime import datetime print(datetime.now()) n = int(1e4) m = int(1e3) intList = [] for j in range(m): ftot = 0 if j % int(m/10) == 0: print(int(j/10), '%', sep = '') for i in range(n): xx = np.random.uniform() ftot = ftot + fcn(xx) intList = intList + [ftot/n] print(datetime.now()) plt.figure() nbins = 20 plt.hist(intList, nbins, color = 'c', edgecolor = 'k') plt.xlabel('integral') plt.ylabel('counts') # The width of the distribution can be calculated from the standard # deviation of the list of 1000 approximations of the integral and is an # estimate in the uncertainty of our determination of the definite integral. # The standard deviation can be calculated easily using the statistics module. import statistics statistics.stdev(intList) # Finally, the last thing we'll attempt to do is to understand how our # uncertainly (standard deviation) depends on the n. So far, all of our # calculations have used n = 1e4. Now we'll determine the standard # deviation for values of n that range from 100 to 1e5. This # block of code took my laptop about 18 minutes to complete... # If you've got a lot time to kill, trying running the same loop in Maple! print(datetime.now()) nList = [100, 200, 500, 1000, 2000, 5000, 1e4, 2e4, 5e4, 1e5] sigmaList = [] for k in range(len(nList)): intList = [] print(k+1, 'of', len(nList)) n = int(nList[k]) for j in range(m): ftot = 0 for i in range(n): xx = np.random.uniform() ftot = ftot + fcn(xx) intList = intList + [ftot/n] sigmaList = sigmaList + [statistics.stdev(intList)] print(datetime.now()) # Below we plot the uncertainty in the numerical integral estimation as # a function of n (the number of trials in the Monte Carlo simulation). # As expected, the uncertainty decreases as the number of trials increases. # The sigma values are proportional to 1/sqrt(n). plt.figure() plt.plot(nList, sigmaList, 'bo', fillstyle = 'none') plt.xlabel('n') plt.ylabel(r'$\sigma$ (std)') plt.axis((-5e3, 1.2e5, 0, 0.04)) # To get a better appreciation of the dependence of sigma on n, below we # plot sigma as a function of 1/sqrt(n) and observe the linear relationship # between the two. Beautiful! All of this generated from uniformly # distributed random numbers! Take a moment to reflect on what we've # accomplished. We've used Monte Carlo simulations to study the behaviour # of Monte Carlo simulations! The objective of a Monte Carlo calculation # is always to study the characteristics of some system (often a physical # system) by simulating data using random numbers. plt.figure() plt.plot(1/np.sqrt(nList), sigmaList, 'bo', fillstyle = 'none') plt.xlabel(r'$1/\sqrt{n}$') plt.ylabel(r'$\sigma$ (std)') plt.axis((0, .11, 0, 0.04))