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)

Reads nicely @simeoncarstens !

Here are a few thoughts that I had while reading the beginning (will continue another day). Some of them are resolved by reading it several times, but maybe it is a good hint for where a reader might hang a bit:

A sufficient (but not neccessary), convenient condition for a Markov chain to have a distribution $p(x)$ as its invariant distribution is called detailed balance: $p(x)T(y|x) = p(y)T(x|y)$.

what is y in the equation?

A distribution $\pi$ of the state space of a Markov Chain

Calling the stationary distribution eigenvector \pi confuses me because I associate it so much with a number. But maybe it is the standard. In that case never mind ...

So all that is fun, but back to sampling arbitrary probability distributions. If we could design our transition kernel in such a way that the next state is already drawn from $\pi$, we would be done, as our Markov Chain would... well... immediately sample from $\pi$. Unfortunately, to do this, we need to be able to sample from $\pi$, which we can't. Otherwise you wouldn't be reading this, right?

I feel this is a nice paragraph but it could be easier to understand without prior knowledge: You don't really introduce what "design our transition kernel" means. Also, how can I "draw from \pi" I thought I only draw a new state based on the row in T of the previous step?

n $q(x_{i+1}^*|x_i)$, from

the ^* is difficult to see and read. Probably better if it has different support. Otherwise, can you use a different symbol?

probability $p_\mathrm{acc}(x_{i+1}=x_{i+1}^|x_{i+1}^, x_i)$

For this transition kernel to obey detailed balance, we need to choose $p_\mathrm{acc}$ correctly. One possibility is the Metropolis acceptance criterion:

how can we know when it is set correctly? Is this difficult to grasp intuitively?

I don't think the method you use for the mixture of Gaussians can be called Gibbs sampling. As you say, in Gibbs sampling you'd have to draw from p(k|x), but it does not equal p(k) because k and x are not independent. You'd have to use the Bayes theorem to calculate p(k|x).

Your method is to draw hierarchically: first draw k from the marginal distribution, then draw x from the conditional distribution, whereas in Gibbs sampling you draw from both conditional distributions, each time using the value from the previous step.

Your method works, but it works for a different reason than the reason Gibbs sampling works. I don't think it has a name.

Edit: moreover, I think your method works well precisely because you are not using GIbbs sampling. If you did, 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.

Finished reading — very nice! I corrected a few typos here: https://gist.github.com/feuerbach/5ccaeee166b45be13a8e375f980f405a.

Maybe one day I'll make a similar post except everything is done in Stan :)

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.