Monte Carlo and “the power of randomness”

Monte Carlo and “the power of randomness”#

What you will learn

Monte Carlo = Random Sampling to Estimate Integrals
- Replace integration with averaging over random samples.
Expectation Viewpoint: $ $\int p (x) f (x) d x = E_{p} [f (x)] \approx \frac{1}{n} \sum_{i = 1}^{n} f (x_{i})$ $
- sample from $x_{i} \sim p (x)$ , and compute averages of $f (x)$ .
Law of Large Numbers:
- As $n \to \infty$ , sample average converges to the true expectation.
Boltzmann Sampling in Statistical Mechanics:
- Sample from $p (x) \propto e^{- β E (x)}$ to estimate thermodynamic averages: $ $⟨ E ⟩ \approx \frac{1}{n} \sum E (x_{i})$ $
Rejection Sampling:
- Sample $x \sim Uniform$ , accept if $y \leq f (x)$ .
- Useful for sampling from unnormalized distributions.
Markov Chain sampling
- generate random samples in such a way where each sample depends on previous one.
- This allows clevery moving into areas of interest.

Estimating the pi by throwing pebbles on the sand#

The key idea of this technique is that the ratio of the area of the circle to the square area that inscribes it is $π / 4$ , so by counting the fraction of the random points in the square that are inside the circle, we get increasingly good estimates to $π$ .

\frac{V_{c i r c l e}}{V_{s q u a r e}} = \frac{π r^{2}}{(2 r)^{2}} = \frac{π}{4}

circle_pi_estimate(N=100000, r0=100)

../_images/78df3b9e59c96c592d96982a43e741ed7a3c9d4e5e6eb9ce4b6e9411f561d302.png

Shapes more complex than a circle#

I = e^{- x} + e^{- x^{2}} x^{2} c o s (x)^{2} + e^{- 2 x} x^{4} c o s (2 x)^{2}

We will now use the same technique but compute a 1D definite integral from $x_{1}$ to $x_{2}$ by drawing a rectangle to cover the curve with dimensions $x = [x_{1}, x_{2}]$ and $y = [a, b]$ .
The area of the rectangle is simply $A$ =2. The area under the curve is $I$ .
If we choose a point uniformly at random in the rectangle, What’s the probability that the point falls into the region under the curve? It is obviously

p = \frac{n_{i n}}{N} \approx \frac{I}{A}

Thus we can estimate definite integral by drawing N uniform numbers covering range and computing integral as $I = A \frac{n_{i n}}{N}$

Text(0, 0.5, 'f(x)')

../_images/1bf6586fe52153f16b46615f0a221b429af0a450ba286963165d459f142653df.png

Exact result: 2.2898343018663505
MC result: 2.2748

../_images/9113b8386eba1a82e6fc51f7bbc33645b4fee92418b70a15a56133f92054e9a9.png

The Essence of Monte Carlo Simulations#

Suppose we want to evaluate an integral $I$ .

I = \int f (x) d x

A powerful perspective is to reinterpret this integral as the expectation of a function $g (x)$ under some probability distribution $p (x)$ :

I = \int \frac{f (x)}{p (x)} p (x) d x = \int g (x) p (x) d x = E_{p} [g (x)]

In this form, the integral becomes the expected value of $g (x)$ with respect to the distribution $p (x)$ .
To estimate $E_{p} [g]$ , we draw samples $x_{i} \sim p (x)$ and apply the law of large numbers, which guarantees that the sample average converges to the expected value as $n \to \infty$ :

E_{p} [g] \approx \frac{1}{n} \sum_{i = 1}^{n} g (x_{i}), where x_{i} \sim p (x)

Simple 1D applications of MC#

Ordinary Monte Carlo and Uniform Sampling#

A common and intuitive case is when we draw samples uniformly from the interval $[a, b]$ . In this setting, the sampling distribution is constant: $p (x) = \frac{1}{b - a}$ , and the integral simplifies as follows:

I = \int_{a}^{b} f (x) d x = (b - a) \int_{a}^{b} \frac{f (x)}{b - a} d x \approx (b - a) \cdot \frac{1}{n} \sum_{i = 1}^{n} f (x_{i}) = (b - a) \cdot {\bar{f}}_{n}

This gives a clear interpretation of Monte Carlo integration: we approximate the average height of the function $f (x)$ over the interval by randomly sampling points, much like tossing pebbles onto a plot and estimating the shaded area.

MC result 2.277565400031629
Exact result: 2.2898343018663505

Sampling from the Boltzmann Distribution#

Quantities like average energy, heat capacity, or pressure are computed as ensemble averages under the Boltzmann distribution.
The Boltzmann distribution for a system with energy $E (x)$ at inverse temperature $β = 1 / (k_{B} T)$ is:

p (x) = \frac{e^{- β E (x)}}{Z}, where Z = \int e^{- β E (x)} d x

Suppose we are interested in the average energy:

⟨ E ⟩ = \int E (x) p (x) d x

This is an expectation value under the Boltzmann distribution $p (x)$ .
If we can draw samples $x_{i} \sim p (x)$ , we can estimate $⟨ E ⟩$ by:

⟨ E ⟩ \approx \frac{1}{n} \sum_{i = 1}^{n} E (x_{i})

Estimated ⟨E⟩: 0.4896
Exact ⟨E⟩:     0.5000

../_images/192f04f48e4790a08a4c527fb619573121f44c1798dd870d2ff838cd76ea3a3c.png

More dimensions#

Calculate are under 2D Gaussian over [-a, a] and [-b, b] region

f (x, y) = \frac{1}{2 π σ_{x} σ_{y}} e^{- \frac{x^{2}}{2 σ_{x}^{2}} - \frac{y^{2}}{2 σ_{y}^{2}}}

(0.9946076967722628, 0.9959510667787033)

Why Does Monte Carlo Outperform Brute-Force Integration?#

Using i.i.d. random variables in Monte Carlo (MC) allows us to apply the Central Limit Theorem and infer that the mean we are calculating $\bar{g}$ will have variance proportional to N number of samples $σ_{N}^{2} = N σ_{1}^{2}$ where $σ_{1}$ is variance in single steps.
Consequently, the convergence rate of Monte Carlo integration is $O (n^{- 1 / 2})$ , which is notable because it is independent of the number of dimensions of the integral.
This property gives Monte Carlo an edge over numerical integration methods, which have a convergence rate of $O (n^{- d})$ , especially in moderate- to high-dimensional contexts. Even in low-dimensional scenarios, Monte Carlo can be advantageous, particularly when the region of interest within the integration space is small. This allows for targeted sampling in critical areas.

On convergence of MC simulations#

We are often interested in knowing how many iterations it takes for Monte Carlo integration to “converge”. To do this, we would like some estimate of the variance, and it is useful to inspect such plots. One simple way to get confidence intervals for the plot of Monte Carlo estimate against number of iterations is simply to do many such simulations.

For the example, we will try to estimate the function (again)

f (x) = x \cos 71 x + \sin 13 x, 0 \leq x \leq 1

Text(0, 0.5, '$f(x)$')

../_images/a66a36bece31729b34fc73fb871a473efe9c2580e07d1b0bd044b2b9ba6f7303.png

We will vary the sample size from 1 to 100 and calculate the value of $y = \sum x / n$ for 1000 replicates. We then plot the 2.5th and 97.5th percentile of the 1000 values of $y$ to see how the variation in $y$ changes with sample size.

../_images/3a266bb4a3a460db29e240e3c3b61e3c9540bc11c342c17af21573a4c8a14469.png

Importance Sampling#

Suppose we want to evaluate the expectation of a function $h (x)$ under a probability distribution $p (x)$ :

I = \int h (x) p (x) d x

If sampling directly from $p (x)$ is difficult, we can instead introduce an alternative distribution $q (x)$ — one that is easier to sample from — and rewrite the integral as:

I = \int h (x) p (x) d x = \int h (x) \frac{p (x)}{q (x)} q (x) d x

This reformulation allows us to draw samples $y_{i} \sim q (x)$ , and weight them by the importance ratio $\frac{p (y_{i})}{q (y_{i})}$ . The expectation can then be estimated using:

I \approx \frac{1}{n} \sum_{i = 1}^{n} \frac{p (y_{i})}{q (y_{i})} h (y_{i})

This is the essence of importance sampling: reweighting samples from an easier distribution to approximate expectations under a more complex one.

../_images/5cac3c05e8cd99393492f5fe718a7efacd37ab942cf7b0f05e8f0f49957c96a6.png

../_images/38401c616373d10d4f507557cf94cef51c53b0860d0067991115c08a39d92c14.png

Using Monte Carlo to sample probability distributions#

Another application of a simple MC technique is to turn uniformly distributed random numbers into random numbers sampled according to different probability distributions.
The key is to employ rejection criteria; if points go under the curve, they are accepted, hence generating a probability distribution

../_images/acddbada84fd95e119ce4b9959df475d76dd02f6a6b0aa16f6498ddc6b1899cf.png

Markov Chain Monte Carlo#

mcmc_pi(N=10000, delta=0.1)

../_images/edeaa0004fd27886cad27b7aac4d52b7313a1fc70db0e565f5c06aeb548cdb0c.png

3.4307178631051753

Text(0.5, 1.0, 'Mean: 3.2182, Var: 0.9890')

../_images/bd3d510ac69df7509fc2df82f12faf35bfd5274338ff9321c5477013f2ccc657.png

Problems#

MC, the crude version#

Evaluate the following integral $\int_{0}^{\infty} \frac{e^{- x}}{1 + (x - 1)^{2}} d x$ using Monte Carlo methods.

Start by doing a direct monte carlo on uniform interval.
Try an importance sampling approach using en exponential probability distribution.
Find the optimal value of $λ$ that gives the most rapid reduction of variance [Hint: experiment with different values of $λ$ ]

MC integral of 3D and 6D spheres!#

Generalize the MC code above for computing the volume of 3D and 6D spheres.
The analytical results are known: $V_{3 d} = \frac{4}{3} π r^{3}$ and $V_{3 d} = \frac{1}{6} π π^{3} r^{6}$ . So you can check the statistical error made in the simulations.