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 $V$ of volume $|V|$. Then
$$ I=\int_V f(x)\,dx =|V|\,\mathbb E[f(X)] \approx |V|\,(\overline f \pm \frac{s_f}{\sqrt N}), $$
where
$$ \overline f=\frac1N\sum_{i=1}^N f(X_i). $$
For independent samples and finite variance, the estimated standard error is
$$ 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.
Intuitive understanding: $$\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)
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]. $$
The variance is
$$ \operatorname{Var}_w\!\left(\frac{f}{w}\right) =\int \frac{f(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)}{e^{-X/\lambda}/\lambda}\right]. $$
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.random(n_points) # uniform mesh
# Inverse cumulative-distribution transformation.
X = -np.log(1-W)
L=5
x_plot = np.linspace(0.0, L, 300)
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")
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.random(batch)
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.7877969081567248 +- 0.0029581700619980607 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=100.
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.8159931791725863 +- 0.06603726047683531 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))
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)
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: 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.
n_points = 200_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))
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.hist(gaussian_samples, bins=100, density=True, alpha=0.65)
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}. $$
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} "
f"+/- {gaussian_error:.6f}"
)
print(f"exact value: {exact:.6f}")
Gaussian importance sampling: 1.633041 +/- 0.001283 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-I^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)= \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¶
For independent, identically distributed variables with finite mean $\mu$ and variance $\sigma^2$, the sample mean approaches a Gaussian distribution as the sample size $n$ grows:
$$ \overline X\approx \mathcal N\!\left(\mu,\frac{\sigma^2}{n}\right). $$
For the uniform distribution on $[0,1)$,
$$ \mu=\frac12, \qquad \sigma^2=\frac1{12}. $$
The theorem explains the familiar $1/\sqrt n$ statistical uncertainty. For correlated Markov-chain samples, the effective number of independent samples is smaller; that issue will be discussed in the Markov-chain lecture.
n_experiments = 50_000
n_draws = 100
samples = rng.random((n_experiments, n_draws))
sample_means = np.mean(samples, axis=1)
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.hist(sample_means, bins=80, density=True, alpha=0.65)
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.