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)\))

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.

Uniform distribution

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

Listing 91 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 distrbution 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.

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 92 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;
}