Skip to article frontmatterSkip to article content

Chapter 12: The Monte Carlo Method

Learning Objectives

By the end of this lecture, you should be able to

  1. Understand the basic principles of the Monte Carlo method.
  2. Apply Monte Carlo sampling to estimate integrals.
  3. Explain the concept of importance sampling and its benefits.

Introduction to Monte Carlo Method

The Monte Carlo method is a powerful computational technique used to solve problems by relying on random sampling. To illustrate its basic idea, let’s start with a simple example: evaluating the integral of f(x)=x2 f(x) = x^2 over the interval [0,1][0, 1]

I=01x2dxI = \int_0^1 x^2 \, dx

Relating Integrals to Averages

Recall that the average value of a function f(x) f(x) over an interval [a,b][a, b] is given by

f(x)=1baabf(x)dx\langle f(x) \rangle = \frac{1}{b - a} \int_a^b f(x) \, dx

In this case, we can write the integral as the average of x2 x^2 over [0,1][0, 1]

I=x2=01x2dxI = \langle x^2 \rangle = \int_0^1 x^2 \, dx

Monte Carlo Estimation

Instead of solving this integral analytically, we can approximate it by randomly sampling points from the interval [0,1][0, 1] and evaluating x2 x^2 at those points. This is the essence of the Monte Carlo method.

Suppose we generate N N random numbers xi x_i uniformly distributed between 0 and 1. We can approximate the integral as the average of xi2 x_i^2

I1Ni=1Nxi2I \approx \frac{1}{N} \sum_{i=1}^N x_i^2

This gives us an estimate for the integral. The larger the number of random samples N N , the closer the estimate gets to the true value.

Python Implementation

The following Python code demonstrates how to estimate the integral using Monte Carlo sampling.

import numpy as np

# Set random seed for reproducibility
np.random.seed(42)

# Number of random samples
N = 1000

# Generate N random numbers uniformly distributed between 0 and 1
x = np.random.rand(N)

# Estimate the average of x^2 over [0, 1]
average_x_squared = np.mean(x**2)

# Estimate the variance of the estimator
variance = np.var(x**2)

# Print the result with a descriptive message
print(f"The estimated average of x^2 over [0, 1] with {N} samples is: {average_x_squared:.6f}")
print(f"The estimated variance of the estimator is: {variance:.6f}")
print(f"The exact value of the integral is: {1/3:.6f}")
The estimated average of x^2 over [0, 1] with 1000 samples is: 0.325610
The estimated variance of the estimator is: 0.089466
The exact value of the integral is: 0.333333

Importance Sampling

In the previous example, we used a uniform distribution to generate random numbers for our Monte Carlo simulation. However, we can often improve the efficiency of these simulations by employing a different probability distribution—this approach is known as importance sampling.

Motivation for Importance Sampling

Let’s reconsider the goal of estimating the integral

I=01x2dxI = \int_0^1 x^2 \, dx

Using a uniform distribution, the Monte Carlo estimator would directly sample from this integral. However, efficiency can be improved by selecting a distribution that matches the behavior of the integrand, i.e., a distribution that emphasizes regions where x2x^2 is larger. This results in a smaller variance of the estimator.

Rewriting the Integral

We can rewrite the integral as

I=01x2g(x)g(x)dxI = \int_0^1 \frac{x^2}{g(x)} g(x) \, dx

where g(x)g(x) is any probability distribution (with support over [0, 1]) that can be used to sample points more efficiently. The idea is to choose a g(x)g(x) that is similar to the shape of the function x2x^2. This way, the sampling is concentrated in regions where the function x2x^2 contributes more to the integral.

Choosing a Suitable g(x)g(x)

A good candidate for g(x)g(x) is the Beta distribution with parameters α=3\alpha = 3 and β=1\beta = 1. The probability density function of the Beta distribution is

g(x)=xα1(1x)β1B(α,β)g(x) = \frac{x^{\alpha - 1} (1 - x)^{\beta - 1}}{B(\alpha, \beta)}

where B(α,β)B(\alpha, \beta) is the Beta function, and in this case, it biases the samples toward larger values of xx, aligning well with the shape of x2x^2.

Monte Carlo Estimator with Importance Sampling

Using the Beta distribution, the integral can be approximated as

I1Ni=1Nxi2g(xi)I \approx \frac{1}{N} \sum_{i=1}^N \frac{x_i^2}{g(x_i)}

where xix_i are samples drawn from the Beta distribution.

Implementation in Python

The following code demonstrates how to estimate the integral of x2x^2 using importance sampling.

import numpy as np
from scipy.stats import beta

# Set random seed for reproducibility
np.random.seed(42)

# Number of random samples
N = 100

# Parameters of the Beta distribution (biasing towards larger x)
alpha, beta_param = 3, 1

# Generate N random samples from the Beta distribution
x = np.random.beta(alpha, beta_param, size=N)

# Evaluate the integrand at the random numbers (x^2)
f_x = x**2

# Evaluate the importance sampling weights
g_x = beta.pdf(x, alpha, beta_param)  # PDF of the Beta distribution

# Since the target is the uniform distribution on [0, 1], its PDF is constant 1
# So the weights become 1/g(x)
weights = 1 / g_x

# Estimate the integral using importance sampling
estimated_integral = np.mean(f_x * weights)

# Estimate the variance of the estimator
variance = np.var(f_x * weights)

# Print the result with a descriptive message
print(f"The estimated value of the integral of x^2 over [0, 1] using importance sampling is: {estimated_integral:.6f}")
print(f"The estimated variance of the estimator is: {variance}")
print(f"The exact value of the integral is: {1/3:.6f}")
The estimated value of the integral of x^2 over [0, 1] using importance sampling is: 0.333333
The estimated variance of the estimator is: 1.4544622940012404e-32
The exact value of the integral is: 0.333333