Open In Colab

Monte Carlo Integration

Introduction

In this tutorial, you will explore the use of the Monte Carlo method to compute integrals. This will begin by illustrating the Metropolis method to sample from an arbitrary probability distribution. Then you will gain some experience in estimating the value of integrals with the Monte Carlo method, including the use of importance sampling.

In the accompanying notebook MC-Lennard-Jones.ipynb, you will perform a Monte Carlo simulation in the \(NVT\) ensemble for a shifted-truncated Lennard-Jones fluid. The two notebooks are independent of each other, but some basic ideas are introduced here, so we recommend that you start with this one.

A third notebook, MC-Pressure.ipynb, deals with \(NPT\) simulations of the Lennard-Jones fluid. This will be covered in a later workshop.

Preliminaries

Start by importing some useful Python modules, setting up the random number generator, and loading a plotting style (feel free to change this as you wish).

[ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad
from numpy.random import default_rng
rng = default_rng()
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis')
plt.rc('legend',frameon=True,framealpha=1.0)
plt.rc('hist',bins=100) # Default number of bins to use in histograms

Metropolis sampling

In this section, you will examine the use of the Metropolis method to sample points from an arbitrary probability distribution. The code in the following cells selects values for the variable \(x\) from the probability distribution \begin{equation*} P(x) = 12 \, \left( x - \tfrac{1}{2} \right)^{2} \end{equation*} where \(x \in [0,1]\). This function is normalized on that range, although this is not necessary in order to do the sampling. In fact, we are going to omit the normalization factor from the function definition, simply to make this point. This is important when we come to do practical Metropolis Monte Carlo: in the canonical ensemble for many atoms, for instance, the normalization factor for the configurational distribution function is very hard to calculate (although, hopefully, you know what it is called?).

[ ]:

def P(x):
    """Un-normalized probability function. Argument x may be scalar or NumPy array."""
    return (x-0.5)**2

This implementation of the Metropolis method starts from an initial point \(x_0 \in [0,1]\), and generates trial points randomly and uniformly in that range. Suppose that \(0,1,2,\ldots, k\) trials have been carried out so far. The next trial point \(x_t\) is accepted or rejected by comparing \(P(x_t)\) with \(P(x_k)\):

  • if \(P(x_t) \ge P(x_k)\), accept the trial;

  • if \(P(x_t) < P(x_k)\), accept the trial with probability \(P(x_t)/P(x_k)\), otherwise reject.

Both the above tests can be combined in the following way, using a random number \(0<u<1\) sampled from the standard uniform distribution:

  • If \(\dfrac{P(x_t)}{P(x_k)} \ge u\) accept the trial, otherwise reject.

Accepting the trial means \(x_{k+1}=x_t\). Rejecting the trial means \(x_{k+1}=x_k\).

Hopefully, it is clear that, in a large set of values of \(x\) sampled like this, those with high \(P(x)\) will occur more often than those with low \(P(x)\). It can be shown that the expected frequency of occurrence of \(x\) is indeed proportional to \(P(x)\).

For convenience, a function metropolis_1d is defined to do this. The argument list consists of

  • the starting point, x0,

  • the range of sampling xrange=(xmin,xmax)=(0.0,1.0) in this case,

  • the sampling function prob which will be \(P(x)\) when we call the function,

  • the number of points to be sampled, n.

The function builds up a list of sampled points, and returns it as a NumPy array.

[ ]:

def metropolis_1d ( x0, xrange, prob, n ):
    """Carries out 1D sampling of specified probability function.

    Arguments
    ---------
    x0 : float, scalar
        starting point, must lie within xrange
    xrange : tuple of scalar floats (xmin,xmax)
        specified range
    prob : function
        specified probability function prob(x) of variable x
    n : int, scalar
        specified number of samples

    Returns
    -------
    float, NumPy 1d array X of length n
        contains the sampled x points
    """

    assert xrange[0] < x0 < xrange[1], 'x0 must be in xrange'
    assert prob(x0) > 0,               'x0 must have prob(x0)>0'
    xk = x0
    X = []
    for _ in range(n):
        xt = rng.uniform(*xrange) # uniform random number in xrange
        if prob(xt)/prob(xk) >= rng.uniform():
            xk = xt
        X.append(xk)
    X = np.array(X)
    return X

This is how we use the function: call it with the desired arguments, including the probability function that we specified earlier.

[ ]:

n      = 10000
xrange = (0.0,1.0)
x0     = rng.uniform(*xrange) # random initial value in range
X      = metropolis_1d(x0,xrange,P,n)

The next cell plots a probability histogram of the sampled values of \(x\), along with the exact \(P(x)\). The histogram is normalized numerically by selecting the density=True option. In plotting the exact \(P(x)\), the correct normalizing factor is now included, since we want to compare visually with the simulated distribution.

[ ]:

fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$P(x)$')
ax.set_xlim(xrange)
ax.hist(X,range=xrange,density=True,label='sampled distribution')
x = np.linspace(*xrange,501)
Pnormed = 12*P(x)
ax.plot(x,Pnormed,label='exact distribution')
ax.legend()
plt.tight_layout()

From the sample, we can estimate averages and moments for this distribution, such as the mean value \(\langle x\rangle\) and the variance \(\langle x^2\rangle-\langle x\rangle^2\), calculated in the following cell. Again, we emphasize that this has been done without normalizing \(P(x)\). Longer runs would give a more precise estimate. For comparison, the exact values are 0.5 and 0.15 respectively, which you can confirm algebraically by calculating integrals over \(P(x)\), but this will involve knowing the normalizing factor!

[ ]:

print(f'Mean value of x = {X.mean():10.5f}')
print(f'Variance of x   = {X.var():10.5f}')

Now it’s your turn. Your task is to sample, in a similar way, points from the distribution \begin{equation*} P(x) = \Bigl(\frac{2}{L}\Bigr) \sin^{2} \Big( \frac{3 \pi x}{L} \Big) \end{equation*} in the range \(x \in [0,L]\) where \(L=5\). This function is defined in the next cell. Note that, in order to re-use the metropolis_1d routine, P is still a function of the single argument x, and the value of L is simply inherited from outside (not necessarily good practice in general!). Again, the normalization factor \((2/L)\) is omitted from the definition, but must be remembered when the exact function is plotted.

[ ]:

def P(x):
    """Un-normalized probability function. Argument x may be scalar or NumPy array."""
    return np.sin(3*np.pi*x/L)**2

Over to you. Add statements to the next cell, defining xrange and x0, and using this function in a suitable call to metropolis_1d. Put the results, once more, in the array X.

[ ]:

n = 10000
L = 5.0
# Insert your code here

The following cell compares the histogram of your X with the exactly known function. They should agree quite well, if the previous cell is correct.

[ ]:

fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$P(x)$')
ax.set_xlim(xrange)
ax.hist(X,range=xrange,density=True,label='sampled distribution')
x = np.linspace(*xrange,501)
Pnormed = (2/L)*P(x)
ax.plot(x,Pnormed,label='exact distribution')
ax.legend()
plt.tight_layout()

Now we aim to sample the two-dimensional distribution \begin{equation*} P(x,y) = \Bigl(\frac{4}{L_x L_y}\Bigr) \sin^{2} \Big( \frac{2 \pi x}{L_{x}} \Big) \, \sin^{2} \Big( \frac{3 \pi y}{L_{y}} \Big) \end{equation*} over the range \(x\in [0,L_x]\), \(y\in [0,L_y]\), where \(L_{x}=3\) and \(L_{y}=5\). Again, the values of Lx and Ly are inherited from outside by the function given below.

[ ]:

def P(x,y):
    """Un-normalized probability function. Arguments x, y may be scalars or NumPy arrays (of the same shape)."""
    px = np.sin(2*np.pi*x/Lx)**2
    py = np.sin(3*np.pi*y/Ly)**2
    return px*py

The approach is to write a new version of metropolis_1d, called metropolis_2d, for the 2D case. The skeleton of this function is given below, and it is up to you to insert the main loop. Two empty lists X and Y are defined, ready to contain the \(x\) and \(y\) coordinates of the sampled points. At each step, the current point is \((x_k,y_k)\). For the next sample, a new trial point \((x_t,y_t)\) should be chosen, randomly and uniformly within the specified area and a decision taken as to whether to accept it or reject it, based on the ratio \(P(x_t,y_t)/P(x_k,y_k)\). Acceptance means that \((x_k,y_k)\) is replaced by \((x_t,y_t)\); rejection means that \((x_k,y_k)\) is left unchanged. Then these coordinates are appended to the lists, and the process moves on to the next step. The function converts the lists to arrays at the end, before returning.

[ ]:

def metropolis_2d ( x0, y0, xrange, yrange, prob, n ):
    """Carries out 2D sampling of specified probability function.

    Arguments
    ---------
    x0, y0 : two floats, scalar
        starting point, must lie within xrange, yrange
    xrange : tuple of scalar floats (xmin,xmax)
        specified range in x
    yrange : tuple of scalar floats (ymin,ymax)
        specified range in y
    prob : function
        specified probability function prob(x,y) of 2 variables x, y
    n : int, scalar
        specified number of samples

    Returns
    -------
    a tuple of two NumPy 1d arrays X, Y, of floats, length n
        X contains the x coordinates of the sampled points
        Y contains the y coordinates of the sampled points
    """

    assert xrange[0] < x0 < xrange[1], 'x0 must be in xrange'
    assert yrange[0] < y0 < yrange[1], 'y0 must be in yrange'
    assert prob(x0,y0) > 0,  '(x0,y0) must have prob(x0,y0)>0'
    xk, yk = x0, y0
    X, Y = [], []
    # Insert your code here
    X, Y = np.array(X), np.array(Y)
    return X, Y

Once you have inserted the necessary statements, run the program in the following cell. Notice that the number of samples has been increased to \(1000000\), because the dimensionality has gone up. This will take several seconds to run, so while you are testing you might like to temporarily reduce that number.

[ ]:

n      = 1000000
Lx, Ly = 3.0, 5.0
xrange = (0.0,Lx)
yrange = (0.0,Ly)
x0, y0 = rng.uniform(*xrange), rng.uniform(*yrange)
X, Y   = metropolis_2d(x0,y0,xrange,yrange,P,n)

Once you are happy, and have carried out the \(1000000\)-sample run, you can compare with with the theoretical result. For our purposes it is sufficient to produce side-by-side plots of the sample histogram (using hist2d from matplotlib), and contours for the theory (using contourf from matplotlib). Hopefully they look similar.

[ ]:

fig, (ax1,ax2) = plt.subplots(1,2,sharey=True,figsize=(7,5))
ax1.set_xlabel(r'$x$')
ax1.set_ylabel(r'$y$')
ax2.set_xlabel(r'$x$')
ax1.hist2d(X,Y,bins=100,range=[xrange,yrange],density=True)
x   = np.linspace(*xrange,301)
y   = np.linspace(*yrange,501)
x,y = np.meshgrid(x,y)
Pnorm = (2/Lx)*(2/Ly)*P(x,y)
ax2.contourf(x,y,Pnorm)
plt.tight_layout()

Monte Carlo integration

In this section, the Monte Carlo method is used to numerically integrate a simple function: \begin{equation*} I = \int_{0}^{L} dx \, F(x) \qquad\text{where}\quad F(x) = e^{-x} \sin x \end{equation*} and \(L=100\). The exact result is \begin{equation*} I = \tfrac{1}{2} \left[ 1 - e^{-L}(\cos L + \sin L) \right] \approx \tfrac{1}{2} \end{equation*} and for practical purposes, replacing \(L=100\) by \(L\rightarrow\infty\) would make no difference, as the integrand is vanishingly small for large \(x\). However, as will become apparent, for most of this exercise a finite value of \(L\) is needed. Note, in passing, that simple quadrature, rather than Monte Carlo, is almost always preferable for numerical estimation of low-dimensional integrals such as this; it is just for illustration.

Uniform sampling Monte Carlo

The first approach is to estimate the integral using the straightforward Monte Carlo method, uniformly sampling \(x \in [0,L]\). Formally, the integral is being rewritten as \begin{equation*} I = L \times \left[ \frac{1}{L} \int_{0}^{L} dx \, F(x) \right] = L \times \langle F(x) \rangle = L \times \langle e^{-x} \sin x \rangle \end{equation*} where \(\langle\cdots\rangle\) represents an average over the uniform distribution \begin{equation*} P(x)= \begin{cases} 1/L & 0 < x < L \\ 0 & \text{otherwise} \end{cases} \end{equation*} This is implemented in the next cells. This code actually performs s MC simulations, each consisting of n sampled points. The s estimates of the integral from the simulations are stored in the array Iuni. The spread of results will give a guide to the (im)precision of a typical result from a single simulation. The direct sampling of \(x\) values can be coded up very compactly.

[ ]:

def F(x):
    """Function to be integrated"""
    return np.exp(-x)*np.sin(x)

[ ]:

L = 100.0
xrange = (0.0,L)
s = 1000
n = 1000

Iuni = []
for _ in range(s):
    X = rng.uniform(*xrange,size=n) # Uniformly sampled x values in xrange
    I = L*np.average(F(X))          # Estimate of integral from above formula
    Iuni.append(I)                  # Append this estimate to list
Iuni = np.array(Iuni)

# From these simulation estimates, calculate a mean value and standard deviation
print(f'Uniform sampling: I = {Iuni.mean():8.5f}',u'\u00B1',f'{Iuni.std():8.5f}')
# Compare exact value
I = 0.5*(1.0-np.exp(-L)*(np.sin(L)+np.cos(L)))
print(f'Exact value:      I = {I:8.5f}')

The distribution of estimates (and hence, the uncertainty in any one of them) can be visualized by creating a histogram of these values, and the following code does this.

[ ]:

fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel('Integral')
ax.set_ylabel('Histogram')
ax.set_xlim((0.0,1.0)) # Sensible likely range of estimates of integral
ax.hist(Iuni,range=(0.0,1.0),label='uniform sampling') # Histogram of estimates of integral
ax.axvline(x=I,c='k',label='exact value') # Exact value of integral
ax.legend()
plt.tight_layout()

Importance sampling

In real life, we might do just one Monte Carlo estimate of such an integral, so the spread of results here is a little worrying, even if it is centred on the correct answer. This method can be improved by introducing importance sampling. Rather than uniformly sampling the entire possible range of \(x\), attention is focused on regions where the value of the integrand, \(e^{-x}\sin x\), is significant. In this case, it might be worth sampling \(x\) from the normalized probability distribution \(P(x) = e^{-x}\). It is convenient (but not essential) to take \(L\rightarrow\infty\). Formally, the integral is being rewritten as \begin{equation*} I = \int_{0}^{\infty} dx \, P(x) \frac{F(x)}{P(x)} = \left\langle \frac{F(x)}{P(x)} \right\rangle_P \end{equation*} where \(\langle\cdots\rangle_P\) indicates an average over the distribution \(P(x)\). The uniform sampling discussed in the previous section is just a special case of this formula, with \(P(x)=1/L\) over the applicable range \([0,L]\) in that case.

Note that, in this type of integration, we do need to know the properly normalized \(P(x)\) in order to calculate the quantity being averaged, \(F(x)/P(x)\)! (In most Monte Carlo simulations, we are typically only calculating ratios of integrals of this kind, so the normalizing factor once again is not needed).

Finally, note that (in this case) the factor \(e^{-x}\) within \(F(x)\) happens to cancel exactly with \(P(x)\), so we end up just calculating \(\langle\sin x\rangle_P\), but this is not essential for the method to work.

For this simple case, the points can be sampled directly: a suitable random number generator is even built in to the NumPy library. This is implemented in the code below, which otherwise looks very similar to the previous cell, and stores the results in the array Iexp. The plot compares the spread of results with Iuni.

[ ]:

def P(x):
    """Probability function. Argument x may be scalar or NumPy array."""
    return np.exp(-x)
def FoverP(x):
    """Function to be integrated divided by probability function"""
    return F(x)/P(x)

[ ]:

s    = 1000
n    = 1000
Iexp = []
for _ in range(s):
    X = rng.standard_exponential(size=n) # Exponentially sampled x values
    I = np.average(FoverP(X))            # Estimate of integral
    Iexp.append(I)                       # Append to list
Iexp = np.array(Iexp)

# From these estimates, calculate a mean value and standard deviation
print(f'Exponential sampling: I = {Iexp.mean():8.5f}',u'\u00B1',f'{Iexp.std():8.5f}')
# Compare exact value
I = 0.5
print(f'Exact value:          I = {I:8.5f}')

[ ]:

fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel('Integral')
ax.set_ylabel('Histogram')
ax.set_xlim((0.0,1.0)) # Sensible likely range of estimates of integral
ax.hist(Iuni,range=(0.0,1.0),label='uniform sampling')
ax.hist(Iexp,range=(0.0,1.0),alpha=0.8,label='exponential sampling')
ax.axvline(x=I,color='k',label='exact value') # Exact value of integral
ax.legend()
plt.tight_layout()

Markov chain random walk

The exponential (importance) sampling should be significantly better, in the sense that the numerical results have a much narrower distribution around the exact value. In the more general case, however, it is not possible to sample the points directly. One alternative approach is to construct a Markov chain random walk.

In this context, a chain is a sequence of coordinates \(x_0,x_1,x_2,\ldots, x_k, \ldots\) (we specialize to 1D for simplicity here). In a Markov chain, the probability of generating the next value \(x_{k+1}\) depends solely on the current value \(x_k\), not on the prior history of the chain. The metropolis_1d code presented earlier produces a Markov chain: the trial value \(x_t\) (chosen randomly and uniformly in the range \(x_\text{min}<x_t<x_\text{max}\)), is accepted (\(x_{k+1}=x_t\)) or rejected (\(x_{k+1}=x_k\)) by comparing \(P(x_t)\) with \(P(x_k)\). So we could use this for importance sampling based on the \(P(x)=e^{-x}\) probability function, or any other probability distribution.

In this section, we are going to refine this a little, and generate uniformly sampled trial values \(x_t\) close to \(x_k\), such that \(x_k-\Delta < x_t < x_k+\Delta\), where \(\Delta\) is a relatively small number (the maximum step size). This step-by-step process is often called a random walk. An advantage of the method is that the probability of acceptance of such trial values can usually be made higher by reducing \(\Delta\), since then \(x_t\approx x_k\) and so \(P(x_t)\approx P(x_k)\). Provided some care is taken with the method of selecting the trial move, the same Metropolis criterion may be applied, to generate the desired distribution \(P(x)\). A finite length system is used, bounded by \(x_\text{min}\) and \(x_\text{max}\) which are stored in xrange; trials are rejected if they fall outside this range, as well as if they fail the Metropolis criterion.

[ ]:

def random_walk_1d ( x0, xrange, delta, prob, n ):
    """Carries out 1D sampling by random walk of specified function.

    Arguments
    ---------
    x0 : float, scalar
        starting point, must lie within xrange
    xrange : tuple of scalar floats (xmin,xmax)
        specified range
    delta : float, scalar
        maximum displacement
    prob : function
        specified probability function prob(x) of variable x
    n : int, scalar
        specified number of samples

    Returns
    -------
    float, NumPy 1d array X of length n
        contains the sampled x points
    """

    assert xrange[0] < x0 < xrange[1], 'x0 must be in xrange'
    assert prob(x0) > 0,               'x0 must have prob(x0)>0'
    xk = x0
    X = []
    for _ in range(n):
        xt = xk + rng.uniform(-delta,delta) # Random walk step
        if xrange[0] < xt < xrange[1]:
            if prob(xt)/prob(xk) >= rng.uniform():
                xk = xt
        X.append(xk)
    X = np.array(X)
    return X

The random walk is implemented in the cells below, again attempting to sample from the \(P(x)=\exp(-x)\) distribution, to integrate \(F(x)=\exp(-x)\sin x\). We re-use the functions P(x) and FoverP(x) defined above.

Each simulation consists, once more, of n=1000 sampled points, within the range \(x\in[0,L]\) where \(L = 100\). The starting point x0 is chosen at the midpoint of the range \([0,L]\). Do you think that this is a good idea? Maybe we should re-examinine this afterwards. Once more the program performs s=1000 MC simulations, so the distribution of answers gives an idea of precision. In the plot, comparison is made with the other two methods: uniform sampling and importance sampling directly from the exponential distribution.

[ ]:

s       = 1000
n       = 1000
delta   = 2.5
L       = 100.0
xrange  = (0.0,L)
x0      = 0.5*L

Iwlk = []
for _ in range(s):
    X = random_walk_1d(x0,xrange,delta,P,n) # Exponentially sampled x values
    I = np.average(FoverP(X))               # Estimate of integral
    Iwlk.append(I)                          # Append to list
Iwlk = np.array(Iwlk)

# From these estimates, calculate a mean value and standard deviation
print(f'Random walk: I = {Iwlk.mean():8.5f}',u'\u00B1',f'{Iwlk.std():8.5f}')
# Compare exact value
I = 0.5*(1.0-np.exp(-L)*(np.sin(L)+np.cos(L)))
print(f'Exact value: I = {I:8.5f}')

[ ]:

fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel('Integral')
ax.set_ylabel('Histogram')
ax.set_xlim((0.0,1.0))
ax.hist(Iuni,range=(0.0,1.0),label='uniform sampling')
ax.hist(Iexp,range=(0.0,1.0),label='exponential sampling')
ax.hist(Iwlk,range=(0.0,1.0),alpha=0.8,label='random walk')
ax.axvline(x=I,color='k',label='exact value')
ax.legend()
plt.tight_layout()

There should be some good points, and some bad points, in the above plot. Spend some time thinking about the questions below, particularly the first two, perhaps experimenting with the code in the cells above.

  1. How does the precision of the random walk method compare with that of the uniform MC method, and the direct exponential sampling method? (Precision relates to the spread in results). How about the accuracy? (This relates to how close the results are to the correct answer).

  2. Concerning the systematic inaccuracy, could the starting position for the random walk be chosen in a better way? Make a more sensible choice, and re-run the random walk.

  3. How do the precision and accuracy of the estimate of the integral vary with the number of points sampled in the simulation?

  4. How does the precision of the random walk method vary with the delta parameter? (We expect this to be connected with the fraction of trials accepted, which is not calculated above).

Comments

You may wish to insert a few comments regarding the above questions in this cell.

Two-dimensional importance sampling

If time permits, use one of the above approaches to numerically estimate the value of the integral \begin{equation*} I = \int_{-\infty}^{\infty} dx \int_{-\infty}^{\infty} dy \, F(x,y) \quad \text{where} \quad F(x,y) = \exp \bigl[ -(x^{2}+y^{2})^{2}+x^{2} \bigr] . \end{equation*} The value is \(I\approx 3.94476\), obtained by standard quadrature, as we can see in the following cell which uses dblquad, imported from the scipy.integrate sub-package at the top of this notebook.

[ ]:

def F(x,y):
    """Function to be integrated."""
    return np.exp(-(x**2+y**2)**2+x**2)

I, Ierr = dblquad ( F, -np.inf, np.inf, -np.inf, np.inf )
print(f'Quadrature value = {I:8.5f}')
print(f'Estimated error  = {Ierr:8.1e}')

Do this by importance sampling, using the double-Gaussian probability density \begin{equation*} P(x,y) = p(x) p(y) \quad\text{where}\quad p(x) = \frac{1}{\sqrt{\pi}} e^{-x^{2}} \end{equation*} This means that we are evaluating the integral as follows \begin{equation*} I = \int_{-\infty}^{\infty} dx \int_{-\infty}^{\infty} dy \, P(x,y) \, \frac{F(x,y)}{P(x,y)} = \left\langle \frac{F(x,y)}{P(x,y)}\right\rangle_P \end{equation*} where \(\langle\cdots\rangle_P\) denotes an average over \(P\). Three approaches spring to mind. If time is short, feel free to postpone thinking about these until after the workshop!

  1. Recognize \(p(x)\) as the normal distribution with mean \(\mu=0\) and variance \(\sigma^2=\frac{1}{2}\), i.e. standard deviation \(\sigma=\sqrt{\frac{1}{2}}\). There is a random number function for this in the NumPy library, rng.normal. You can use it to sample \(P(x,y)\) directly.

  2. Do Metropolis sampling, in a square box \(-L < x,y < L\), accepting or rejecting moves according to \(P(x,y)\). The function metropolis_2d above, can do this, without any changes. You’ll need to choose \(L\) large enough that the integrand is essentially zero at the boundary: \(L=2.5\) should be OK.

  3. Do a Markov chain random walk in 2D, adapting the earlier random_walk_1d function into a 2D version random_walk_2d, and otherwise following the general scheme of the Metropolis sampler.

For convenience the necessary additional functions are defined in the following cell. To get an idea of the reliability of the answer, again conduct s=1000 independent simulations, and plot the distribution of results. However, for a 2D example, the number of samples should be increased: n=5000 may be a good starting point. Note: running the Metropolis or random walk sampler may take a few minutes.

[ ]:

def P(x,y):
    """Normalized probability density."""
    return np.exp(-x**2-y**2)/np.pi
def FoverP(x,y):
    """Integrand F divided by probability density P."""
    return F(x,y) / P(x,y)

In the following cells, if time permits or after the workshop, try any or all of the above 3 approaches, comparing the distribution of answers with the quadrature value just calculated.

[ ]:



[ ]:

This concludes the notebook.

[ ]: