48  Metropolis-Hastings

using Plots, LaTeXStrings

Suppose we want to approximate the integral

\begin{align} \int_{\mathbb R^D} f(x) \mu(x) \mathrm{d}x \end{align}

over a high dimensional space (D \gg 1). For example, \mu is the probability density function of some probability distribution \mathbb P and the integral above is the expectation of the random variable f(X) where X \sim \mathbb P.

The Metropolis-Hastings (MH) algorithm is a way of sampling the distribution to obtain x_1,\dots,x_L which allows one to approximate the expectation using the empirical average

\begin{align} \int_{\mathbb R^N} f(x) \mu(x) \mathrm{d}x % \approx \frac{1}{L} \sum_{i=1}^L f(x_i). \end{align}

Quadrature rules (studied last semester) suffer from the Curse of Dimensionality: if one uses n grid points per dimension, then the total number of grid points is n^D (when n=10 (which is small) and D=100 (which isn’t big for e.g. computational materials science where D = 3N with N being the number of particles, and is very small for e.g. machine learning) n^D = 10^100 is larger than the number of particles in the universe).

Another feature of the MH algorithm is that, we don’t actually need to know the normalisation of the probablility distribution to be able to sample from it - can you see why from the code below?

a, b = -5, 5
μ(x) = exp(-x^2) * (2 + sin(x/2) + sin(5x)) / (2√π)

# approximate ∑_i δ( ⋅ - x_i ) with Gaussians
h = 0.2
δ(x) = exp(-(x/h)^2) / (h√π)

# proposal distribution
σ = .25
g() = σ * randn()

# number of iterations
L = 2500
S = [0.5]          

anim = @animate for i ∈ 1:L
    x = S[end]
    x_prop = x + g()

    # Metropolis 
    α = min(1.0, μ(x_prop) / μ(x))
    if rand() < α
        push!(S, x_prop)
        current_color = :green 
    else
        push!(S, x)
        current_color = :red 
    end

    # empirical density 
    μL(val) = sum(δ(val - s) for s ∈ S) / length(S)

    plot(μ, a, b;
        title="Iteration $i",
        label=L"Target $μ(x)$", lw=3, 
        color=:black, ylims=(-.1, 1.))
    plot!(μL, a, b;
        label=L"Empirical $μ_L(x)$, $\sigma=%$σ$", 
        lw=2, color=:blue, fill=(0, 0.2, :blue))
    scatter!([S[end]], [0]; 
        color=current_color, primary=false, markersize=6)
end

gif(anim, fps=10)

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\tmp.gif