MCMC sampling for dummies

When I give talks about probabilistic programming and Bayesian statistics, I usually gloss over the details of how inference is actually performed, treating it as a black box essentially. The beauty of probabilistic programming is that you actually don't have to understand how the inference works in order to build models, but it certainly helps. When I presented a new Bayesian model to Quantopian's CEO, Fawce, who wasn't trained in Bayesian stats but is eager to understand it, he started to ask about the part I usually gloss over: "Thomas, how does the inference actually work? How do we get these magical samples from the posterior?". Now I could have said: "Well that's easy, MCMC generates samples from the posterior distribution by constructing a reversible Markov-chain that has as its equilibrium distribution the target posterior distribution. That statement is correct, but is it useful? My pet peeve with how math and stats are taught is that no one ever tells you about the intuition behind the concepts (which is usually quite simple) but only hands you some scary math. This is certainly the way I was taught and I had to spend countless hours banging my head against the wall until that euraka moment came about. Usually things weren't as scary or seemingly complex once I deciphered what it meant. This blog post is an attempt at trying to explain the intuition behind MCMC sampling (specifically, the Metropolis algorithm). Critically, we'll be using code examples rather than formulas or math-speak. Eventually you'll need that but I personally think it's better to start with the an example and build the intuition before you move on to the math. We have P(\theta x), the probability of our model parameters \theta given the data x and thus our quantity of interest. To compute this we multiply the prior P(\theta) (what we think about \theta before we have seen any data) and the likelihood P(x \theta), i.e. how we think our data is distributed. This nominator is pretty easy to solve for. However, lets take a closer look at the denominator. P(x) which is also called the evidence (i.e. the evidence that the data x was generated by this model). This is the key difficulty with Bayes formula -- while the formula looks innocent enough, for even slightly non-trivial models you just can't compute the posterior in a closed-form way. Now we might say "OK, if we can't solve something, could we try to approximate it?