Monte Carlo Methods

Monte Carlo methods involve using random numbers to produce different realizations of a system and then building up statistics.

Pseudo-random numbers

Generally, producing truly random numbers is hard and we most often use pseudo-random numbers. These should have the following properties:

The are fast to generate
Each number is not correlated with the one that came before it (actually, there are a number of statistical tests they should satisfy).
They uniformly sample the range of interest (usually \([0, 1)\))

Note

True random number generators usually use hardware, noise, or some natural source of statistical randomness

There are a few popular methods for generating pseudo-random numbers, including:

linear congruential generators : these are easy to code up and fast but don’t give the best statistical properties.
Mersenne twister : this method is one of the most widely used now and generates good quality random numbers. It is more complicated to code up.

The C++ standard library has these methods and more for us to use in the random library.

Caution

Some (historical) pseudo-random number generators provided by compiler vendors have been “truly horrible”, such as the IBM RANDU generator.

Uniform distribution

Here’s a simple example of a uniform random number distribution:

Listing 114 test_random.cpp

#include <iostream>
#include <random>

int main() {

    // a uniform random number generator that will be used as our seed
    std::random_device rd;

    // the Mersenne Twister -- we initialize it with a random seed
    std::mt19937 generator(rd());

    // our random number distribution
    std::uniform_real_distribution<double> uniform(0.0, 1.0);

    for (int n = 0; n < 10; ++n) {
        // sample our distribution with our generator
        std::cout << uniform(generator) << std::endl;
    }

}

Some notes:

In C++ <random>, the generator and distribution are separate concepts. Here we pick the Mersenne Twister std::mt19937 generator and use it with the uniform distribution std::uniform_real_distribution.
To initialize the generator, we use std::random_device—this is a true random number generator (depending on system support) and will provide a seed as the starting point for our pseudo-random number generator.

Notice that each time we run it we get a different sequence of random numbers.

Tip

We could alternately use std::default_random_engine by changing the line:

std::mt19937 generator(rd());

to

std::default_random_engine generator(rd());

However, this is implementation dependent, so you may not be sure what engine it is using.

Other distributions

We can create other distributions. Imagine that \(q(z)\) is the probably of getting a number between \(z\) and \(z +dz\). For the uniform distribution, \(q(z) = 1\). We want to find an \(x(z)\) that maps the uniformly-distributed numbers into another distribution. Since the total probability is \(1\), we require:

\[p(x) dx = q(z) dz\]

where \(p(x)\) will be the new distribution function. Integrating this allows us to map the uniform distribution to other distributions.

Tip

The cppreference site has a nice example of a normal distributon that also introduces std::map to hold “key:value” pairs.

Seeding

Sometimes we want reproducible random numbers—i.e. we want the sequence of numbers we get to be random but we want to get the same sequence each time we run. To accomplish this, we can explicitly feed a seed into our random generator:

std::mt19937 generator(12345);

To do better, you can use std::seed_seq (https://en.cppreference.com/w/cpp/numeric/random/seed_seq) to build a better seed:

std::seed_seq seed{1, 2, 3, 4, 5};
std::mt19937 generator(seed);

Warning

For applications that require security, you need to be more careful. But for physics simulations where we are not worried about someone being able to guess our sequence, seeding allows for reproducibility, and is good for testing.

Example: computing \(\pi\)

A fun way to compute \(\pi\) is to pick pairs of random numbers in the range \([0, 1)\)—these will all fall inside a unit square, but we also check to see if they fall inside a circle of diameter 1 inscribed inside the unit square.

We can write

\[\frac{A_\mathrm{circle}}{A_\mathrm{square}} = \frac{\pi (d/2)^2}{d^2} = \frac{\pi}{4}\]

We then pick a large number (\(N\)) of random number pairs. We approximate the area of the circle and square as simply the number of random pairs that fall inside each shape. If we call the number of random numbers that fall inside the circle \(N_\mathrm{circle}\) and note that all \(N\) fall inside the square, so we can approximate \(\pi\) as:

\[\pi \approx 4 \frac{N_\mathrm{circ}}{N}\]

Here’s an implementation of our computing \(\pi\):

Listing 115 pi.cpp

#include <iostream>
#include <random>
#include <cmath>
const int N = 1000000;

int main() {

    std::random_device rd;
    std::mt19937 generator(rd());
    std::uniform_real_distribution<double> uniform(0.0, 1.0);

    int N_circle{0};

    for (int n = 0; n < N; ++n) {

        double x = uniform(generator);
        double y = uniform(generator);

        double r = std::sqrt(std::pow(x - 0.5, 2) + std::pow(y - 0.5, 2));

        if (r <= 0.5) {
            N_circle++;
        }

    }

    std::cout << "pi = " << 4 * static_cast<double>(N_circle) / N << std::endl;
}

Example: hypersphere volume

Consider a \(d\)-dimensional hypersphere with radius 1 defined on a domain \([-1, 1]^d\). The volume is:

\[V_d(R) = \frac{\pi^{d/2}}{\Gamma \left (\tfrac{d}{2} + 1\right )} R^d\]

We can integrate this using Monte Carlo methods.

Mean value Monte Carlo integration works by sampling the domain that the function is defined on to estimate the mean. In one-dimension, we know:

\[\langle f \rangle = \frac{1}{b - a} \int_a^b f(x) dx = \frac{I}{b -a}\]

where \(I\) is the integral we desire:

\[I = \int_a^b f(x) dx\]

This means we can approximate \(I\) by doing \(N\) samples as:

\[I \approx \frac{b - a}{N} \sum_{i = 1}^{N} f(x_i)\]

where \(x_i\) is a randomly selected point in \([a, b]\).

This extends to an arbitrary dimension easily, as:

\[I = \int_\omega f(x_0, x_1, \ldots x_{d-1}) dx_0 dx_1, \ldots dx_{d-1} \approx \frac{V}{N} \sum_{i = 1}^{N} f({\bf x}_i)\]

where \({\bf x}_i\) is a randomly-selected point in the \(d\)-dimensional domain.

For \(d > 3\), Monte Carlo integration is usually less computationally-expensive than doing the trapezoid rule in each dimension.

Here’s an example of computing the volume of the \(d\)-dimensional hypersphere using mean-value Monte Carlo integration:

Listing 116 hypersphere.cpp

// compute the volume of a D-dimensional hypersphere
// via mean-value Monte-Carlo integration

#include <iostream>
#include <iomanip>
#include <random>
#include <cmath>
#include <numbers>

constexpr int N_MAX{1000000};

double integrand(std::vector<double> x) {

    // inner_product() here will compute x . x
    double radius = std::sqrt(std::inner_product(x.cbegin(), x.cend(), x.cbegin(), 0.0));

    if (radius <= 1.0) {
        return 1.0;
    }

    return 0.0;
}

double analytic(const double D) {
    return std::pow(std::numbers::pi, D/2.0) / std::tgamma(D/2.0 + 1.0);
}

int main() {

    std::random_device rd;
    std::mt19937 generator(rd());

    std::uniform_real_distribution<double> distribution(-1.0, 1.0);

    double D{-1};
    std::cout << "Enter dimensionality of the hypersphere: ";
    std::cin >> D;

    // we'll hold D random numbers, representing {r_0, r_1, ..., r_{D-1}}

    std::vector<double> r(D, 0.0);

    int N = 10;
    while (N <= N_MAX) {

        double V{0.0};

        for (int n = 0; n < N; ++n) {

            // get D random numbers between -1 and 1

            for (auto &e : r) {
                e = distribution(generator);
            }

            V += integrand(r);
        }

        V *= std::pow(2.0, D);

        std::cout << "samples: " << std::setw(7) << N << ";  "
                  << "volume: " << V / N << std::endl;

        N *= 10;

    }

    std::cout << "analytic result: " << analytic(D) << std::endl;
}

Tip

This requires compiling with C++20 support due to the use of the <numbers> header.

Here’s the output for \(d = 5\):

Enter dimensionality of the hypersphere: 5
samples:      10; volume = 12.8
samples:     100; volume = 8.32
samples:    1000; volume = 5.504
samples:   10000; volume = 5.2416
samples:  100000; volume = 5.31168
samples: 1000000; volume = 5.26419
analytic result: 5.26379

Note

There are better ways to do Monte Carlo integration, including importance sampling.