Monte Carlo integration¶
Monte Carlo methods use random sampling and statistical estimation to solve problems that may be difficult to treat deterministically. Their applications include high-dimensional integration, Markov-chain simulations, statistical physics, optimization by simulated annealing, and many others.
This notebook introduces pseudorandom numbers and Monte Carlo integration with independent samples. We begin with direct sampling and then reduce the statistical error by importance sampling. The next notebook develops the adaptive VEGAS integration algorithm. Later lectures introduce Markov-chain sampling, the Ising model, simulated annealing, and the Wang–Landau algorithm.
Pseudorandom numbers¶
A deterministic computer produces a reproducible sequence of pseudorandom numbers from an internal state. A good pseudorandom-number generator is designed so that the sequence passes demanding statistical tests and has a sufficiently long period. This is entirely adequate for most Monte Carlo calculations, but the choice of generator and the construction of independent streams matter in large simulations.
NumPy's modern interface is Generator:
rng = np.random.default_rng(seed)
The fixed seed makes this notebook reproducible. Omitting it initializes the generator with entropy supplied by the operating system. NumPy's current default bit generator is PCG64. PCG, Wiki.
Before spending a large amount of computer time, it is prudent to verify that the result is stable within its statistical uncertainty and, for especially demanding calculations, to repeat it with an independent seed or another high-quality generator.
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
%matplotlib inline
rng = np.random.default_rng(509)
How does a generator work?¶
A simple example is the linear congruential generator
$$ x_{n+1}=(a x_n+c)\bmod m. $$
Carelessly chosen parameters can produce a short period, poor low-order bits, and conspicuous correlations: points in several dimensions may lie on a relatively small number of hyperplanes.
PCG combines a congruential update of a larger internal state with a permutation of the output bits. The permutation hides the simple algebraic structure of the state evolution and substantially improves the statistical properties of the output. The result is still pseudorandom: the same initial state produces exactly the same sequence.
For parallel calculations, one should create independent streams through a facility such as SeedSequence.spawn, rather than giving every worker the same seed.
Direct Monte Carlo integration¶
Suppose that $X$ is sampled uniformly from a region of volume $V$. Then
$$ I=\int_V f(x)\,dx \approx V\,(\overline f \pm \frac{s_f}{\sqrt N}), $$
where
$$ \overline f=\frac1N\sum_{i=1}^N f(X_i). $$
and $s_f^2$ is the sample variance
$$ s_f^2=\frac{1}{N-1}\sum_{i=1}^N\bigl(f(X_i)-\overline f\bigr)^2. $$
The exponent $N^{-1/2}$ does not explicitly depend on dimension. The variance—and therefore the prefactor—can nevertheless grow strongly with dimension.
Why does the error of the Monte Carlo average scale as $1/\sqrt{N}$? $$\delta f_1+\delta f_2+\cdots+\delta f_N \sim \sqrt N\,\sigma_f,$$ hence $$ \delta\overline f \sim \frac{\sqrt N\,\sigma_f}{N} = \frac{\sigma_f}{\sqrt N}.$$
Why Monte Carlo is useful in high dimension¶
If a tensor-product quadrature uses $n$ points along each of $d$ coordinates, it requires
$$ N=n^d $$
function evaluations. A second-order rule has an illustrative error scaling
$$ n^{-2}=N^{-2/d}. $$
This becomes less favorable than $N^{-1/2}$ when $d>4$. The comparison is only a guide: the actual crossover depends on smoothness, error constants, coordinate transformations, sparse grids, and the variance of the Monte Carlo estimator.
If the desired region is difficult to sample directly, it can be enclosed in a simpler region and the integrand set to zero outside the desired region. Rejected points then increase the variance but do not bias the result.
Example 1: estimating $\pi$¶
We generate random points uniformly inside the square
$$ -1\le x\le1,\qquad -1\le y\le1. $$
The square has area $4$, while the unit circle has area $\pi$. Therefore, the probability that a random point lies inside the circle is
$$ P(\text{inside}) = \frac{\text{area of circle}}{\text{area of square}} = \frac{\pi}{4}. $$
In a numerical experiment, this probability is approximated by the fraction of sampled points inside the circle:
$$ P(\text{inside}) \approx \frac{N_{\text{inside}}}{N_{\text{total}}}. $$
Consequently,
$$ \boxed{ \pi\approx 4\frac{N_{\text{inside}}}{N_{\text{total}}} }. $$
n_samples = 1_000_000
x = rng.uniform(-1.0, 1.0, size=n_samples)
y = rng.uniform(-1.0, 1.0, size=n_samples)
inside = x**2 + y**2 <= 1.0
pi_estimate = 4.0*np.mean(inside)
pi_error = 4.0*np.std(inside, ddof=1)/np.sqrt(n_samples) # ddof just means sdt with N/(N-1) (<f^2>-<f>^2)
print(f"pi = {pi_estimate:.6f} +/- {pi_error:.6f}")
print(f"exact value = {np.pi:.6f}")
pi = 3.139208 +/- 0.001644 exact value = 3.141593
Example 2: a Gaussian over the unit disk¶
We calculate
$$ I=\int_{x^2+y^2\le1}e^{-x^2-y^2}\,dx\,dy =2\pi\int_0^1 r e^{-r^2}\,dr =\pi(1-e^{-1}). $$
r_squared = x**2 + y**2
values = np.exp(-r_squared)*(r_squared <= 1.0)
integral = 4.0*np.mean(values)
error = 4.0*np.std(values, ddof=1)/np.sqrt(n_samples)
exact = np.pi*(1.0 - np.exp(-1.0))
print(f"I = {integral:.6f} +/- {error:.6f}")
print(f"exact value = {exact:.6f}")
I = 1.983391 +/- 0.001221 exact value = 1.985865
Importance sampling¶
Direct sampling is inefficient when the integrand is appreciable only in a small part of the integration region. Let $w(x)$ be a normalized probability density with $w(x)>0$ wherever $f(x)\ne0$. Then
$$ I=\int f(x)\,dx =\int \frac{f(x)}{w(x)}w(x)\,dx =\mathbb E_w\!\left[\frac{f(X)}{w(X)}\right]. $$ Here $\mathbb E_w$ means the expected value, with $X$ sampled according to the probability density $w$. The variance is
$$ \operatorname{Var}_w\!\left(\frac{f}{w}\right) =\int \frac{f(x)^2}{w(x)^2} w(x)\,dx-I^2. $$
The purpose of importance sampling is to choose $w$ so that $f/w$ fluctuates much less than the direct-sampling estimator.
Inverse-transform sampling¶
Define the cumulative distribution
$$ W(x)=\int_{-\infty}^{x}w(t)\,dt. $$
If $U$ is uniform on $[0,1)$ and the inverse is available, then
$$ X=W^{-1}(U) $$
has probability density $w(x)$. Indeed,
$$ \frac{dP}{dx} =\frac{dP}{dW}\frac{dW}{dx} =1\cdot w(x). $$
Exponential sampling¶
For
$$ w(x)=\frac1\lambda e^{-x/\lambda},\qquad x\ge0, $$
the cumulative distribution and its inverse are
$$ W(x)=1-e^{-x/\lambda}, \qquad x=-\lambda\ln(1-W). $$
Therefore \begin{equation} \int_0^\infty f(x)dx=\int_0^1 \frac{f(-\lambda\ln(1-W))}{w(-\lambda\ln(1-W))} dW = \int_0^1 f(-\lambda\ln(1-W))\frac{\lambda dW}{1-W} \end{equation} is easily evaluated with MC if $f(x)$ is exponentially falling function. \begin{eqnarray} \int_0^\infty f(x)dx\rightarrow \langle \lambda \frac{f(-\lambda\ln(1-W))}{1-W}\rangle_{W\;uniform\in[0,1]} \end{eqnarray}
Since $\frac{dP}{dW}=1$ ($W$ is uniformly distributed), the probability for $x$ is $\frac{dP}{dx}=w(x)$, therefore $x$ is exponentially distributed random number.
We could also write \begin{eqnarray} \int_0^\infty f(x)dx\rightarrow \left\langle\frac{f(x)}{\frac{e^{-x/\lambda}}{\lambda}}\right\rangle_{\frac{dP}{dx}=e^{-x/\lambda}/\lambda} \end{eqnarray} or $$ \int_0^\infty f(x)\,dx =\mathbb E_w\!\left[\frac{f(x)}{\frac{e^{-x/\lambda}}{\lambda}}\right] $$ with $w=\frac{e^{-x/\lambda}}{\lambda}$.
NumPy provides this distribution directly through rng.exponential, but deriving the transformation makes the change of probability measure explicit.
n_points = 200_000
exponential_samples = rng.exponential(scale=1.0, size=n_points)
W = rng.uniform(0, 1, size=n_points) # uniform mesh
# Inverse cumulative-distribution transformation.
X = -np.log(1-W)
L=5
x_hist = np.linspace(0.0, L, 80)
plt.hist(X, bins=x_hist, density=True, alpha=0.65, label="transformed samples")
plt.hist(exponential_samples, bins=x_hist, density=True, alpha=0.65,label="samples")
x_plot = np.linspace(0.0, L, 300)
plt.plot(x_plot, np.exp(-x_plot), "r", label=r"$e^{-x}$")
plt.xlabel("x")
plt.ylabel("probability density")
plt.legend()
plt.show()
Example 3: exponential importance sampling¶
Consider
$$ I=\int_0^\infty e^{-x}\frac{\sin x}{x}\,dx =\arctan(1)=\frac{\pi}{4}. $$
With $w(x)=e^{-x}$.
def f(x):
return np.exp(-x)*np.sin(x)/x
batch = 100
faver = 0
Nall = 0
f2 = 0
for i in range(100):
W = rng.uniform(0, 1, size=batch)
X = -np.log(1-W)
faver += sum(f(X)/np.exp(-X))
Nall += batch
f2 += sum((f(X)/np.exp(-X))**2)
f2 /= Nall
faver /= Nall
print(faver, '+-',np.sqrt((f2-faver**2)/Nall),'exact=', 0.785398)
0.7900291006964224 +- 0.00292408184316791 exact= 0.785398
For comparison, direct sampling must replace the infinite interval by a finite one. Taking $L=20$ makes the neglected exponential tail much smaller than the statistical error in this example.
batch = 100
faver = 0
Nall = 0
L=20.
f2 = 0
for i in range(100):
X = rng.uniform(0, L, size=batch)
faver += np.sum(f(X))
Nall += batch
f2 += np.sum(f(X)**2)
f2 /= Nall
faver /= Nall
print(L*faver,'+-', L*np.sqrt((f2-faver**2)/Nall), 'exact=', 0.785398)
0.7745936264669585 +- 0.028060953871350828 exact= 0.785398
Same scaling, different variance¶
Importance sampling usually changes the variance—the prefactor of the error—not the independent-sampling scaling $N^{-1/2}$. We can demonstrate this by repeating both estimators many times and measuring their root-mean-square error.
sample_sizes = np.array([100, 300, 1_000, 3_000, 10_000, 30_000])
n_repetitions = 100
L = 100.0
exact = np.pi/4
direct_rmse = []
importance_rmse = []
for size in sample_sizes:
# Direct sampling: X is uniform on [0, L].
X = rng.uniform(0.0, L, size=(n_repetitions, size)) # two dimensional array (repetitions, size)
direct_trials = L*np.mean(f(X), axis=1)
# Importance sampling: W is uniform on [0, 1], while
# X = -log(1-W) has probability density w(X) = exp(-X).
W = rng.uniform(0.0, 1.0, size=(n_repetitions, size))
X = -np.log(1.0 - W)
importance_trials = np.mean(f(X)/np.exp(-X), axis=1)
# how far from exact results is the approximation?
direct_rmse.append(np.sqrt(np.mean((direct_trials - exact)**2)))
importance_rmse.append(np.sqrt(np.mean((importance_trials - exact)**2)))
direct_rmse = np.asarray(direct_rmse)
importance_rmse = np.asarray(importance_rmse)
reference = importance_rmse[0]*np.sqrt(sample_sizes[0]/sample_sizes)
plt.loglog(sample_sizes, direct_rmse, "o-", label="direct sampling")
plt.loglog(sample_sizes, importance_rmse, "o-", label="importance sampling",)
plt.loglog(sample_sizes, reference, "k--", label=r"$N^{-1/2}$")
plt.xlabel("number of samples N")
plt.ylabel("root-mean-square error")
plt.legend()
plt.show()
Gaussian sampling¶
How to generate Gaussian probability?
The Gaussian integral has no elementary expression; it is written in terms of the error function. Its inverse also requires a special function or a numerical approximation.
The Box–Muller transformation¶
To generate two independent Gaussian variables, draw $U_1,U_2$ uniformly from $[0,1)$ and define
$$ X_1=x_0+\sigma\sqrt{-2\ln(1-U_1)}\cos(2\pi U_2), $$
$$ X_2=x_0+\sigma\sqrt{-2\ln(1-U_1)}\sin(2\pi U_2). $$
The Jacobian gives the joint density
$$ p(X_1,X_2)= \prod_{i=1}^{2} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left[-\frac{(X_i-x_0)^2}{2\sigma^2}\right]. $$
Thus $X_1$ and $X_2$ are independent normal variables with mean $x_0$ and variance $\sigma^2$. In ordinary work one would call rng.normal; here we implement the transformation to see where Gaussian random numbers come from.
Why does the Box–Muller transformation produce Gaussian random numbers?¶
The essential point is that $U_1$ and $U_2$ are independent random variables uniformly distributed on the unit square. Their joint probability density is therefore
$$ p_U(u_1,u_2)=1, \qquad 0\le u_1,u_2<1. $$
The Box–Muller transformation maps these two uniform variables into
$$ X_1 = x_0+\sigma\sqrt{-2\ln(1-U_1)}\cos(2\pi U_2), $$
$$ X_2 = x_0+\sigma\sqrt{-2\ln(1-U_1)}\sin(2\pi U_2). $$
We want to determine the joint probability density $p_X(x_1,x_2)$ of the transformed variables.
Probability must be preserved by the change of variables:
$$ p_X(x_1,x_2)\,dx_1\,dx_2 = p_U(u_1,u_2)\,du_1\,du_2. $$
Therefore,
$$ p_X(x_1,x_2) = p_U(u_1,u_2) \left| \frac{\partial(u_1,u_2)} {\partial(x_1,x_2)} \right|. $$
The Jacobian of the inverse transformation is the reciprocal of the Jacobian of the forward transformation:
$$ \left| \frac{\partial(u_1,u_2)} {\partial(x_1,x_2)} \right| = \frac{1}{ \left| \frac{\partial(x_1,x_2)} {\partial(u_1,u_2)} \right| }. $$
Since $p_U(u_1,u_2)=1$, we obtain
$$ \boxed{ p_X(x_1,x_2) = \frac{1}{ \left| \frac{\partial(x_1,x_2)} {\partial(u_1,u_2)} \right| } }. $$
The forward Jacobian is
$$ \frac{\partial(x_1,x_2)} {\partial(u_1,u_2)} = \begin{vmatrix} \dfrac{\partial x_1}{\partial u_1} & \dfrac{\partial x_1}{\partial u_2} \\[3mm] \dfrac{\partial x_2}{\partial u_1} & \dfrac{\partial x_2}{\partial u_2} \end{vmatrix}. $$
It is a simple excercise to compute this determinant and express the result entirely in terms of $x_1$ and $x_2$.
n_points = 20_000
sigma = 1.0
x0 = 0.0
u1 = rng.random(n_points)
u2 = rng.random(n_points)
radius = sigma*np.sqrt(-2.0*np.log1p(-u1))
x1 = x0 + radius*np.cos(2*np.pi*u2)
x2 = x0 + radius*np.sin(2*np.pi*u2)
gaussian_samples = np.concatenate((x1, x2))
plt.hist(gaussian_samples, bins=100, density=True, alpha=0.65)
x_plot = np.linspace(-5.0, 5.0, 400)
normal_density = np.exp(-(x_plot - x0)**2/(2*sigma**2))/ np.sqrt(2*np.pi*sigma**2)
plt.plot(x_plot, normal_density, "r", label="Gaussian density")
plt.xlabel("x")
plt.ylabel("probability density")
plt.legend()
plt.show()
Example 4: Gaussian importance sampling¶
Consider
$$ I=\int_{-\infty}^{\infty}e^{-x^2}\frac{\sin x}{x}\,dx =\pi\,\operatorname{erf}\!\left(\frac12\right). $$
Choose
$$ w(x)=\frac{e^{-x^2}}{\sqrt\pi}, $$
which is a normal distribution with variance $1/2$. The estimator becomes
$$ \frac{f(x)}{w(x)}=\sqrt\pi\frac{\sin x}{x}. $$
def sinx_over_x(x):
"""Return sin(x)/x, including the limiting value 1 at x=0."""
return np.sinc(x/np.pi)
n_samples = 20_000
x_gaussian = rng.normal(loc=0.0, scale=1/np.sqrt(2), size=n_samples)
gaussian_values = np.sqrt(np.pi)*sinx_over_x(x_gaussian)
gaussian_estimate = np.mean(gaussian_values)
gaussian_error = np.std(gaussian_values, ddof=1)/np.sqrt(n_samples)
exact = np.pi*erf(0.5)
print(f"Gaussian importance sampling: {gaussian_estimate:.6f} +/- {gaussian_error:.6f}")
print(f"exact value: {exact:.6f}")
Gaussian importance sampling: 1.634857 +/- 0.001276 exact value: 1.635199
What is the optimal importance distribution?¶
For a normalized positive density $w$, the variance is
$$ \sigma^2[w]=\int \frac{f(x)^2}{w(x)}\,dx- \left(\int \frac{f(x)}{w(x)} w(x) dx\right)^2. $$
Minimize it subject to $\int w(x)\,dx=1$:
$$ \delta\left[ \int \frac{f(x)^2}{w(x)}\,dx +\lambda\left(\int w(x)\,dx-1\right) \right]=0. $$
Because $\delta w(x)$ is arbitrary,
$$ -\frac{f(x)^2}{w(x)^2}+\lambda=0, $$
and therefore $w_{\mathrm{opt}}(x)\propto |f(x)|$ or
$$ w_{\mathrm{opt}}(x)= \frac{|f(x)|}{\int |f(t)|\,dt}. $$
For a nonnegative integrand this would give zero variance, but its normalization is the unknown integral itself. It is therefore an ideal target rather than usually an available solution. In practice we choose a normalized approximation to $|f|$ and may improve it adaptively. For sign-changing functions, sampling $|f|$ leaves fluctuations of the sign; severe cancellation leads to the Monte Carlo sign problem.
The central limit theorem¶
Let $X_1,X_2,\ldots,X_n$ be independent random variables drawn from the same probability distribution. Suppose that this distribution has mean
$$\langle X_i\rangle =a$$
and variance
$$\operatorname{Var}(X_i)=s^2,$$
where $s$ is the standard deviation.
Consider the sum of $n$ samples,
$$ S_n=X_1+X_2+\cdots+X_n. $$
The mean and variance of this sum are
$$ \langle S_n \rangle = n a, \qquad \operatorname{Var}(S_n)=ns^2. $$
Therefore, the standard deviation of the sum is
$$ \operatorname{std}(S_n)=\sqrt{n}\,s. $$
These relations are exact. The central limit theorem adds that, for sufficiently large $n$, the probability distribution of the sum approaches a Gaussian distribution:
$$ S_n \approx \mathcal{N}(na,ns^2). $$
The original distribution of the $X_i$ does not need to be Gaussian. It must, however, have a finite mean and variance, and the samples must be independent, or at least sufficiently weakly correlated.
Why Monte Carlo error decreases as $1/\sqrt{n}$¶
In Monte Carlo calculations, we are usually interested in the sample average rather than the sum:
$$ \overline{X}_n =\frac{1}{n}\sum_{i=1}^{n}X_i =\frac{S_n}{n}. $$
Its mean is
$$ \langle \overline{X}_n \rangle = a, $$
and its variance is
$$ \operatorname{Var}(\overline{X}_n) =\frac{\operatorname{Var}(S_n)}{n^2} =\frac{ns^2}{n^2} =\frac{s^2}{n}. $$
Consequently, the standard deviation of the sample average is
$$ \operatorname{std}(\overline{X}_n)=\frac{s}{\sqrt{n}}. $$
This standard deviation is the statistical uncertainty of the Monte Carlo estimate. Thus,
$$ \boxed{ \overline{X}_n = a\pm\frac{s}{\sqrt{n}} } $$
for large $n$, where the uncertainty represents approximately one standard deviation.
The important consequence is that Monte Carlo error decreases only as
$$ \boxed{\text{error}\propto\frac{1}{\sqrt{n}}}. $$
To reduce the statistical error by a factor of $10$, we therefore need approximately $100$ times more independent samples.
Consider uniform distribution in the interval $x\in [0,1)$. We know $$\langle x \rangle = 1/2$$ and $$ Var(x) = \langle x^2\rangle - \langle x\rangle^2 = \frac{1}{12} .$$ Then, according to the central limit theorem, the sum of $n$ independent samples approaches a Gaussian distribution with mean $n/2$ and variance $n/12$. The sample average approaches a Gaussian distribution with mean $1/2$ and variance $1/(12n)$.
n_experiments = 50_000
n_draws = 100
samples = rng.random((n_experiments, n_draws))
sample_means = np.mean(samples, axis=1)
plt.hist(sample_means, bins=80, density=True, alpha=0.65)
mean_variance = 1.0/(12*n_draws)
x_plot = np.linspace(0.35, 0.65, 400)
normal_approximation = np.exp(-(x_plot - 0.5)**2/(2*mean_variance))/ np.sqrt(2*np.pi*mean_variance)
plt.plot(x_plot, normal_approximation, "r", label="CLT Gaussian")
plt.title(f"Distribution of the mean of {n_draws} uniform samples")
plt.xlabel("sample mean")
plt.ylabel("probability density")
plt.legend()
plt.show()
Adaptive integration: MISER and VEGAS¶
Stratified sampling divides the integration region into subregions and allocates more points where the estimated contribution to the variance is large. The MISER algorithm applies this idea recursively, choosing subdivisions intended to reduce the total variance.
VEGAS instead learns a separable importance distribution from earlier samples. It is especially effective when the important structure can be represented well by adaptations along the individual coordinate directions. Narrow structures along correlated combinations of coordinates can be more challenging.
The next notebook develops VEGAS in detail. A later lecture treats strongly correlated high-dimensional integrals using Metropolis sampling and an auxiliary normalized distribution.