analytical_hmc.py

from numpy import sin, cos, arctan, sqrt, exp, random, pi, linspace
import matplotlib.pyplot as plt

def draw_sample(xold, sigma):

    t = 3.0
    vold = random.normal()
    phi = arctan(-vold / xold * sigma)
    A = vold * sigma * sqrt(xold ** 2 / sigma ** 2 / vold ** 2 + 1)
    xnew = A * cos(t / sigma + phi)
    vnew = -A / sigma * sin(t / sigma + phi)

    E = lambda x: 0.5 * x ** 2 / sigma ** 2
    K = lambda v: 0.5 * v ** 2
    H = lambda x, v: E(x) + K(v)

    p_acc = min(1, exp(-(H(xnew, vnew) - H(xold, vold))))

    if random.random() < p_acc:
        return xnew, True
    else:
        return xold, False


sigma = 2.0
samples = [2.0]
accepted = 0
n_samples = 100000
for _ in range(n_samples):
    new_state, acc = draw_sample(samples[-1], sigma)
    samples.append(new_state)
    accepted += acc


fig, ax = plt.subplots()
ax.hist(samples, bins=40, density=True)
gaussian = lambda x: exp(-0.5 * x ** 2 / sigma ** 2) / sqrt(2 * pi * sigma ** 2)
xspace = linspace(-5, 5, 300)
ax.plot(xspace, list(map(gaussian, xspace)))
plt.show()

print("Acceptante rate:", accepted / n_samples)

imgs.sh

#!/bin/bash

img_list=$(ls -v output*.png)
b=$(<$2)

while read strA <&3 && read strB <&4; do
    rstring="..\/..\/img\/posts\/${strB}"
    echo $rstring
    sed -i "s/${strA}/${rstring}/g" $1
    mv $strA $strB
    # cp $strB ~/projects/tweag/www/app/assets/img/posts/
done 3<<<"$img_list" 4<<<"$b"

# cp $1 ~/projects/tweag/www/app/views/posts/

mcmc_introduction_gibbs.ipynb

mcmc_introduction_HMC.ipynb

mcmc_introduction_MH.ipynb

mcmc_introduction_MH.zip

mcmc_introduction_RE.ipynb

replace_inlinemath.py

#!/usr/bin/env python3

import sys
from itertools import cycle
import re

with open(sys.argv[1]) as ipf:
    lines = ipf.readlines()
# ## replace \{ with \\{ and \} with \\}
lines = [l.replace('\\{', r'\\{') for l in lines]
lines = [l.replace('\\}', r'\\}') for l in lines]

## replace \\ with \\\\
lines = [l.replace(r' \\', r' \\\\') for l in lines]

## replace ^* with ^\*
lines = [l.replace(r'^*', r'^\*') for l in lines]

## alternatingly replace $ with \\( and \\)
## if it's not part of $$
lines2 = []
for line in lines:
    if '$$' in line:
        lines2.append(line)
        continue
    else:
        cycler = cycle((True, False))
        matches = re.finditer('\$', line)
        offset = 0
        for match in matches:
            replacement = '\\\(' if next(cycler) else '\\\)'
            line = line[:match.start()+offset] + replacement + line[match.start()+1+offset:]
            offset += 2
        lines2.append(line)

with open(sys.argv[2]) as ipf:
    header = ipf.readlines()

with open(sys.argv[3], 'w') as opf:
    for line in header + lines2[2:]:
        opf.write(line)

Thanks for your feedback, @MMesch and @feuerbach! Highly appreciated!
For now I mostly worked on the MH part and put it into a separate notebook, where I addressed some of @MMesch's points.

@MMesch: I fear the \pi is indeed somewhat standard - as p often denotes other probabilities such as p_acc.

@feuerbach: This is very interesting - and you're right. The background to this example is that, a few years ago, I did something like that to sample from a different kind of mixture model and back then we used to call it Gibbs sampling. Looking this up, it seems like a very common method to sample from Gaussian (and other) mixture models (where you can easily marginalize out x to obtain p(k)). It seems to be called "Collapsed Gibbs sampling". I will rewrite this part to feature an actual Gibbs sampler. I just need to think of a nice, concise example where you can sample from the conditional distributions without resorting to MCMC. And as for being stuck in one of the mixture components: highly correlated samples often are a problem when using Gibbs sampling and I should discuss this.

I am also thinking of discussing how to assess convergence for MCMC methods, but this is a very difficult topic. It would make the series much more complete, but is also a really ugly rock to look under...

@feuerbach: Just FYI, I implemented the actual Gibbs sampler for this problem and it turns out that

you'd have the same problem of being stuck in one of the mixture components, because p(k|x) would very likely give you the same k as the one that generated that x.

is a slight understatement. If the sampler starts in the mode on the right, there's (on average) a ~1/1e6 chance for it to jump to the mode on the left. That was fun and very instructive and might even serve in a next version of that blog post / notebook as a negative example.