Monte Carlo Methods

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 77 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.

There are many other distributions we can use.

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 π#

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