MCMC & the art of Sampling without Sampling

Story, Intuition & the gentle Math behind the greatest algorithm of the 20th century

“It is surprising that our simple algorithm worked so well” (it) “simulates what Nature does in looking for the equilibrium of complex systems. And Nature often does its job quite efficiently” — Cerny, after observing the effectiveness of MCMC on the traveling salesman problem.

A raging debate happened about a century back between two leading Math luminaries in Russia. It was on the matter of whether free will existed or not. The debate produced an idea that would fuel the progress of technology — ranging from developing the atom bomb to powering Google search, understanding genes to seeding the foundations of AI Language Models.

Nekrasov was bent on showing that free will existed & he was keen to use Maths to justify his theory. The other Mathematician named Markov was a purist. He did not like mixing Math with philosophy and religion & was bent upon proving Nekrasov wrong. It was well known around that time that the “Law of Large Numbers” existed. Essentially upon repeating an independent event many times, the average outcome approaches the expected value. For e.g. if you toss a fair coin many times, you will eventually end up with a total Head to Tail ratio of 50:50. This may not happen in the first 3, 5 or even 10 tosses. But if you note the outcomes for (say) 1000+ tosses, it will eventually converge to the probability distribution — [Heads: 0.5, Tails: 0.5].

The “Law of Large Numbers” could be seen in many areas — for e.g. in crime statistics or marriage registrations of those times. It was accepted that as long as the events were independent (like in a coin toss), it could come into play. Nekrasov flipped the argument on its head and suggested that because the Law of large numbers kept propping up every where in many diverse aspects of life, therfore the the underlying acts of individuals doing these events had to be independent & voluntary — each acted on their own free will!

Maybe he put his arguments in a more ‘Mathy’ way, but essentially, this was no proof at all. Markov was furious. Leaving aside all his Math pursuits, he did a strange thing. He retired to his study with a copy of Pushkin’s Eugene Onegin & appeared to be reading it day in & out, scribbling furiously with his pencil. Was he frustrated with the likes of Nekrasov & Maths in general & hence wanting to move to literature to appreciate the fine arts instead? Not really! He was intently analyzing the vowels & consonants in the text. He showed he could predict the probability of the next letter being a vowel or a consonant based on the current letter. But why was he doing this?

What Markov showed statistically was that the letters following one another were not independent. There is a certain probability of a vowel following a vowel & a certain probability of a vowel following a consonant & vice-versa. In spite of these dependencies, the law of large numbers still applied to statistics such as the ratio of total number of vowels to consonants in a large text. He thus demolished the free will argument of Nekrasov by showing that Law of Large numbers holds even if the events involved are dependent on one another! In the process, he opened a flood gate for new ideas & thoughts!

You see, until then, the assumption of independence when building such models was a very restrictive one. Most real world situations involve a certain amount of dependency between events & hence scientists hesitated to use the Law of Large numbers for such situations. Markov was telling them they were wrong. Wow! Many real-world situations involving dependent events would now have to be re-looked at and re-modelled!

Unfortunately, that didn’t happen immediately. Markov himself was concerned only with “questions of pure analysis” saying “I refer to the question of the applicability of probability theory with indifference”. As for the rest of the world, maybe they didn’t realize the huge implications. But many years later, a seriously ill man on a hospital bed did. We will come to that story later. Let us now look at the challenges involved in Bayesian inferencing. Only then can we appreciate the solution conjured by the seriously-ill man on the hospital bed. For that solution helped in making the atom bomb. And yes, it helped in Bayesian inferencing and gave a new life to Bayesian statistics.

Bayesian Inference — Challenges

In Bayesian analysis, we make inferences on “unknown quantities of interest”. These unknown quantities could be the parameters of a model, missing data in a dataset, or even predictions. The broad technique to do so has been outlined by Bayes & covered in detail here. We start with a prior assumption about these quantities of interest and next combine this assumption with the evidence seen on the field. This evidence is our training dataset — a set of observations. We end up with what is called a posterior — our prior modified & strengthened by the evidence. We use this posterior for predictions.

Consider a random variable X with a certain probability distribution. We want to know its parameters θ. Specifically, we want to find p(θ|D) i.e. having observed the data-samples D, we want to know which θ is most likely to have generated this data. p(θ|D) is our Posterior. “|” can be read as “given”. For e.g., if the data D was generated by a Gaussian distribution, we can find out its parameters θ (the mean & variance combo) by calculating the posterior. The good Reverend Bayes has provided a neat formula to us for calculating the posterior involving priors & likelihoods & a particularly tiresome term in the denominator called the “normalising constant”. This term p(D) involves an integration across the entire possible θ space which is very difficult to solve!

But we must persist, for this kind of situation involving the posterior calculations has wide applications. Take for instance, the case of object tracking needed to develop autonomous cars. We can observe the (approx) positions of a moving object at different times. We want to simulate the velocity and the exact position given those observations. So our target is p(position, velocity | observed positions). This situation is not too different from predicting from a posterior. Here too, we simulate from a variable (true position & velocity replaces θ) given some observations. Yes, the situation can be cast into a Bayesian setup & the solution involves computing the posterior & using it for inference. So this is a problem that needs solving. What do we do?

Avoidance — The art of fighting without fighting

• In certain situations, we need not compute the normalising constant at all! We can work with the relative values. Basically, since the denominator i.e. the normalizing constant is a constant, if the numerator of p(θ1|D) > the numerator of p(θ2|D), we can say that θ1 better explains the data comapred to θ2. We can select the best θ in a possible range of θ’s.
• In general, Bayes equation can be slightly modified (see fig below) to remove the denominator and replace the “=” sign with a “∝” sign. The left hand side is the un-normalized posterior distribution/density. This guy is not a valid distribution because it will not integrate to 1 but mark my words, this formula will prove very helpful in many situations later!

• For the situation involving calculation of the MAP (used to derive our favorite ML loss functions), the normalizing constant in the denominator does not depend on θ & we can simply drop it from our computations!

• Then there is the neat conjugate prior trick which we covered here.
• In some other situations, we saw how we can side-step calculating this normalizing constant by converting an inference problem into an optimization one & approximate it using Variational Inference technique.

We see many ways to avoid calculating the normalizing constant. How about sampling techniques? Can they avoid calculating the normalizing constant and yet let us do what we want to do? For e.g., if our posterior is a Gaussian & we get many samples, then averaging these would give a value which is approximately equal to the mean of the posterior. So if we could somehow generate enough samples of the posterior, then the parameters of this “samples dataset” is a good approximate for the parameters of the posterior itself. Of course, to sample from a distribution, there is the rather mildly inconvinient pre-requisite of knowing what that distribution is in the first place i.e. to sample from p(θ|D), we need to compute its denominator (the Normalizing Constant) & get the final equation, in order to sample from it.

Still, as it turns out, this(sampling) is the way to go! But before we find out how, did you know that even if we knew the posterior, sampling can turn out to be a pretty ugly task for most distributions? I think we need to understand this first if are to truly appreciate some of the cool sampling algos invented.

Sampling is a tough business…

It is initially difficult to believe that sampling is tough. After all, it just involves randomly selecting values that follow the pattern described by the distribution. For e.g., when we throw a dice, the value we get is a sample. Throw it 100 times & we end up with a dataset of samples with us. Surely, somebody has had the good sense to write a library that allows us (Applied AI practitioners) to sample from any known distribution? Well, as it turns out sampling for even a simple Gaussian, can be difficult (though possible). Can you intuitively think of any way to sample from a Gaussian? Maybe not!

Thankfully, we can generate samples from a Uniform distribution using the standard libraries. In fact, this allows us to further generate samples from a few other simple distributions as well. The algo goes something like this:-

• Say we want to generate samples for a distribution p
• Generate a sample from a continuous Uniform distribution between 0,1
• Now, calculate the sample x of the distribution p whose CDF is the sample generated in the above step. Yes, we need to know the CDF formula for this step & in fact know how to invert it! Remember CDF is nothing but a function that maps a value x to the cumulative probability of observing a value up to x. The Inverse CDF Function — ICDF does the inverse. It takes a (cumulative) probability & returns a sample x!

So if we know the ICDF for a distribution p, then we can sample from a uniform distribution to get us a probability between 0–1 and feed this probability to the ICDF and have it return us a sample from p. Cool, we can then easily generate samples for p. Unfortunately, this ICDF is usually not easy to calculate. Try calculating it for a simple standard Gaussian N(0,1)! So, it will often be the case that we may know the PDF for a distribution but still not know how know how to (easily) sample from it.

We also have other challenges in estimating the posterior via sampling. In high-dimensional spaces, the volume of the space increases exponentially, making data points extremely sparse. Uniform sampling strategy would be extremely inefficient, as almost all samples would land in a low probability region & are useless. How do we get samples from high density zones where it matters? Indeed, these were the problems that plagued Statisticians & AI researchers for a long time. Thanks to the occasional vices of a few brilliant men (including the sick man on the hospital bed), an algorithm called MCMC was invented which was so powerful that it was classified & kept away from human eyes for a long, long time. When it was finally released, it was considered to be amongst the 10 most powerful algos of the 20th century.

MCMC allows us to simulate any posterior distribution p(θ|D) & allows us to sample from it! It is magical because it does this (hold your breath) without even knowing the form of posterior! Once we get the samples, we can estimate virtually all quantities of interest of the posterior directly from the samples. Bayesian statistics would never remain the same again! But we are getting ahead of ourselves. Before going to MCMC, let us see what other ways were invented to estimate the quantities of interest in a complex unknown distribution like a posterior or any complex distriution “p” for that matter.

<The rather lengthy preamble above may help the applied AI practitioner. I first learnt about MCMC by just using it directly in a particular situation. I didn’t have the opportunity to appreciate the genius of the algorithm, nor deeply understand the challenges of the posterior or evaluate other available options!>

Monte Carlo Techniques —Rejection & Importance Sampling

Before MCMC, let us look at a couple of simpler options first. Situation: We have a complex distribution p (usually a posterior but it could be anything) & we want to estimate the parameters of p. If we could get hold of a reasonable number of “samples” of p, we could calculate the parameters for that “samples dataset” & this will approximate the parameters of p. Unfortunately, we don’t know how to sample from p! We don’t know how to calculate the Normalizing Constant in Bayes formula but we can calculate the numerator.

The solution intuitively is to start with a simple proposal distribution q & take samples from here. We then set up an acceptance probability & “accept” certain samples such that it appears as if the final set of (accepted) samples were generated by p, our target distribution.

Rejection Sampling:

First, scale the proposal distribution q big time by multiplying it with a constant c such that q(x) ≥ p(x) for nearly all x. Now sample observations x liberally from q, but reject all samples which are above p(x) — thus resulting exactly in the distribution p eventually.

Source: https://gregorygundersen.com/blog/2019/09/01/sampling/

In the above fig, our target distributiom p is bi-modal — the one with 2 peaks in blue. Assume, this is difficult to sample from. Gundersen uses the standard Gaussian as a proposal distribution q. You may be familiar with the bell shaped Gaussian curve above. He scales this Gaussian big time & then uses Rejection sampling to eliminate all samples (brown) where q(x) > p(x). What remains are only the samples in blue! These samples look as if they are coming from p but they were actually sampled from the Gaussian q. We can now use these samples to approximate p’s parameters.

The actual algo goes something like:

1. Sample from q. Say we get a sample zi
2. Sample ui​ from Uniform(0, c * zi). Note, the range here is from 0 to c * zi
3. If ui​ < p(zi​), Accept zi​ else reject

The uniform sampling at step 2 may seem puzzling but is needed. Consider the 2D scenario below. For a given x’, we want to uniformly sample one value from all possible values on the y-axis. This ensures that at the end of the day all accepted samples are uniformly distributed over the area under p(x).

From a Math perspective, it it easier to re-write step 2 & step 3 of the algo as:

• Step 2: Generate a sample from Unform(0,1). This is ui
• Step 3: If ui​ ≤ p(zi​) / {c * q(zi)}, Accept zi​ else reject

Both algorithms are exactly the same! This one is slightly more convinient since we can talk of the ratio of p(zi​) / {c * q(zi)}.

• This ratio is bounded between 0 and 1
• The number of times N that this algo needs to be run to successfully generate a sample is a random variable with a geometric distribution. Turns out this is exactly = c. Hence c should be kept low.
• As a rule of thumb, if q is chosen such that it resembles p closely, then c can be as low as 1. Of course, this is is only theoritical. In reality, the algo works like a charm for larger values of c to an extent. For very large values of c, inefficiencies arise since we need to reject a lot of samples!
• c grows exponentially with dimensionality, so the acceptance rate will be very small in high dimensions & this solution becomes inefficient!
• Note that we only need to know p “up to a constant of proportionality”. This means we need not know the normalizing constant of p! So if the real p is p/NC1, we can still use the un-normalized p to proceed with the above algorithm & it would still work! The constant c would then just have to factor this in. In fact, note that if we use the un-normalized form of q also (in addition to p), the ratio p(zi​)/c*q(zi) just becomes (p(zi​)/NC1) / {c*(q(zi)/NC2)} which is just p(zi​) / {YetAnotherConstant * q(zi)}. However, do note that normalizing constants are different for p and q & don’t cancel out each other, whatever your favorite AI might say :)

Importance Sampling

Importance sampling is a variation that accepts all samples but assigns importance to samples. This importance is p(x)/q(x) (of course, we need to ensure that q is non-zero when p is non-zero). While this explanation is seen in many articles & is intuitive, we might want to note a few additional points:

• For starters, Importance sampling is not really sampling — in the sense that unlike Rejection sampling, the goal is not to get samples from p
• Rather, it is a direct method to calculate the quantities of interest in a complex distribution p. Say, we want to calculate the Expectation of a function f of a Random Variable X with a distribution p:

• This integral is too complex to compute. Instead we want to take several random samples from p(x), compute f(x) and then average. Sadly, sampling from p is difficult! Importance sampling helps us approximate the above integral & gives us an unbiased estimator of E[f(x)]! How?
• As always, we try to bring in a simple distribution q. We tinker with the Expectation formula slightly. Multiply & divide by q(x)!

So, we converted the original Expectation which involved p into an Expectation involving q (multiplied by a factor p(x)/q(x))! What does this achieve? Well, getting samples from q is a breeze unlike getting samples from p. Just get several random samples from q, weight each by p(x)/q(x) and you can easily compute the Expectation E[f(x)] for the distribution p! Of course, the assumption is that we are able to evaluate p(x) i.e. we are able to calculate the normalizing constant in the denominator to get a normalized PDF for p. The solution can be extended slightly to also work for an un-normalized PDF of p using Self-Normalized Importance Sampling technique.

Both Rejection Sampling & Importance Sampling are “Monte Carlo” algorithms. Any algorithm that constructs solutions based on a large number of samples from a given distribution are referred to as Monte Carlo (MC) methods. An MC simulation typically involves collecting tons of data related to a “process” to understand the deep relationship between input & output. We then create a model that represents the process in question. Then, we experiment with simulations — generate as much simulated inputs we want, feed them to the model & produce simulated outputs. We then use these to understand what we want to understand about the “process” in question. The amount of error can typically be reduced by expending more resources (usually running time and memory). In general, if we can draw enough samples from a distribution (& this is a big IF), we can use these samples to form a MC estimate of the desired quantities of interest.

Importance sampling consumes less resources than Rejection sampling. Unfortunately, it too doesn’t work that well in high dimensions. Our last resort is MCMC! Yea, wedon’t need ICDF, we don’t need the ability to directly sample from p, we dont need to compute the normalizing constant and it also works for higher dimensions!!! But before we go there, note that Importance sampling is important in its own right. It can often lead to a reduction in variance compared to plain Monte Carlo techniques. This variance reduction is achieved by the fact that it wisely chooses samples focussing the important regions which have more impact on parameters being estimated. This leads to reliable estimations with reduced variance.

MCMC

The first MC stands Monte Carlo, so it is one of those algos that are based on random sampling from a distribution to give a solution. The second MC stands for Markov Chains. So it has something to do with our purist, Markov who we read about earlier. Stan Ulam developed the original Monte Carlo simulation method for the Manhattan project — the one where they created the atom bomb. Some say he was a gambler, others maintain that it was his uncle who was a gambler and Ulam only fell in love with Cards while being hospitalised for many weeks for Encephalitis. Since there were no Insta reels to while away his time, he focussed his attention towards how to maximize his chances of winning at Solitaire which led to the Monte Carlo discoveries.

To create the atom bomb, one needed lots of Uranium to experiment. But Uranium was rare to come by. Ulam used realistic simulations based on his learnings during the time he was hospitalised & further honed at the Monte Carlo casino in Las Vegas. These simulations replaced the actual data in the experiments. These techniques were further enhanced by his colleagues like Nicholas Metropolis & MCMC was developed. MCMC is a generic algorithm & there are many implementations , like the Metropolis Hastings, Gibbs etc.

Arianna Rosenbluth wrote the first implementation of MCMC. An Olympic qualified athlete, she was also a great physicist & a terrific programmer. After the MANIACI had been built, she collaborated with Metropolis, Marshall Rosenbluth (her husband) & Augusta/Edward Teller to code the first MCMC implementation the world has seen. Arianna was around till recently — she perished in the COVID wave at the age of 93. Her colleague Mary Tsingou, one of the first programmers of MANIACI is still around at 96. Programming at that time was not the ‘import & use’ variety of yesterday nor the ‘English vibe coding’ of today. Rather, it involved translating workflows into programs on punch cards (FORTRAN had not been invented yet), while using a 5 KB total memory. Debugging IDEs didn’t exist & one had to troubleshoot imagining all the flows & variables in the mind itself. It was no mean achievement by some of the finest & earliest programmers the world has seen.

MCMC uses a couple of core ideas we have seen before (a) Since sampling from p is expensive, start with a simple distribution q & somehow create a situation where the samples coming out from it resemble p (b) Somehow bring in a “ratio of distributions” into the equation so that the normalizing constant cancels out. But there is one more key ingredient needed. Before we go there, we need an idea of Markov Chains and their connection with “probability distributions” in general. We must go back to Markov’s discovery!

Markov Chains:

“He filled a shelf with a small army of books and read and read; but none of it made sense. .. They were all subject to various cramping limitations: those of the past were outdated, and those of the present were obsessed with the past.”
― Alexander Pushkin, Eugene Onegin

Markov predicted the probability of the next letter in Pushkin’s poem being a vowel or consononant. He did this based on the simplest possible dependency situation — just by looking at the current letter. In other words, the model can “forget” all the previous history of alphabets & be memory-less. Only the last alphabet matters! We call this a Markovian property. He found that this was a reasonable assumption & he could use only the current alphabet type to predict the next alphabet type. Before we look at the massive repurcussions of this experiment, let us understand Markov Chains.

Many situations in our daily life can be modelled as a series of transitions betwen States. The transition happens one step at a time. Surprisingly, many of these situations follow the Markovian property i.e. the next state only depends on the current state, and not any historic information. Such situations can be modelled as a Markov Chain (MC). For e.g., if one wanted to model the Weather, then MC is a good choice — There are a bunch of Weather States (rainy, cloudy, sunny etc) with transitions happening between these states on a daily basis. Moreover, if it has rained today it is likely to rain tomorrow as well — a Markovian property. The last state contains info about the history.

Get Allohvk’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

The Markovian property assumption is not very restrictive. We could always enrich the state-space such that the current state includes whatever additional info & context necessary (from the past) to predict the next state i.e. a finite span of history can be encoded into the current state. For e.g., tomorrow’s weather could be made dependent on both yesterday & today’s by creating a state model in which each state is a two-day sequence combo. Of course, the number of States increase exponentially, so it is a tradeoff.

We discussed letters & weather but the State space can be anything — numbers, cricket scores, stock performances whatever. The transition probabilities can be written in a neat n × n matrix form when the chain has n possible states. The (i, j) entry is the probability of transitioning from State_i to State_j. This makes it easy to simulate any time step in the chain — just do a matrix multiplication! So, once we model a situation as an MC, we can make predictions about their long time behaviour!

Let us go back to Markov’s encounter with the fine arts. After manually counting 20K letters in Pushkin’s poem, he found there were 8,638 vowels & 11,362 consonants. This ratio of Vowels to Consonants was in line with the expected value i.e. take any other piece of literature & one ends with a somewhat similar ratio. So Pushkin’s text adhered to the Law of Large numbers! Now all he had to do was to show that the letters weren’t independent!

The probability of a randomly chosen letter being a vowel in Pushkin’s text = 8,638/20K i.e. 0.43. If adjacent letters were independent, then the probability of two vowels in succession would be 0.43 ^ 2, about 3800 instances. Markov showed (by paper, pen & simple counting) that the actual number was 1365 i.e. a third of what was expected if the letters had been independent! This means Pushkin’s text violated the principle of independence even though it adhered to the Law of Large numbers w.r.t Vowels to Consonants ratio.

Markov must have been excited after using an undisputable data-driven approach to demolish Nekrasov’s arguments on free will! But his brief foray into the fine arts was the proverbial butterfly flapping that set loose a chain of events that literally changed the world. Let us look at Markov’s original Chain which was a by-product of his experiments. He took the simplest case of 2 States (current letter is a vowel / current letter is a consonant) & checked the probability of the next letter being Vowel or Consonant. This situation has 4 transitions. He manually calculated the transition probabilities by counting the number of instances in the text for each transition.

The World’s first Markov Chain? (except that it probably wasn’t called so then)

Pushkin’s poem starts with a consonant M “My uncle, a most worthy gentleman,”. Given the above state-transition matrix, what would be the probability of the second letter being a vowel or a consonant? Logically, it is [.526, .474]. This is the probability distribution after one letter i.e. one time-step. We can use a vector-matrix multiplication to get the same result.

What would be the probability distribution after 2 letters?

The chance of being a vowel is .341 & consonant is .659. We get the same results by multiplying the initial input [0,1] with the transition-matrix twice

Above representation is very easy to program. Let us quickly calculate-

• After 1 letter, we saw that the probability distribution is: [0.526, 0.474]
• After 2 letters, we saw that the probability distribution is: [0.341, 0.659]
• After 5 letters the probability distribution is: [0.391, 0.609]
• After 100 letters the probability distribution is: [0.389, 0.611]
• After 1000 letters the probability distribution is: [0.389, 0.611]

Hey, It is not changing! But the real fun part is this. It doesn’t matter whether you start with a consonant or vowel. Let us start with a vowel now:

• After 1 letter the probability distribution is: [0.175 0.825]
• After 2 letters the probability distribution is: [0.465 0.535]
• After 5 letters the probability distribution is: [0.386 0.614]
• After 100 letters the probability distribution is: [0.389, 0.611]
• After 1000 letters the probability distribution is: [0.389, 0.611]

We see the probability distribution becoming stationary after a few time steps & it doesn’t really matter what the initial input was. This is crazy! This means the next letter is always going to be a vowel with probability 0.389 & consonant with p=0.611! Basically, if P is our transition matrix and current probabilities are π, then π_new which is = π * P is π itself i.e. π = πP.

This happens if there are 2 states or 2000 states or even if the distribution is no longer discrete but a continuous one! So here is a leap of imagination! What if we could “create” a Markov Chain whose stationary distribution is p — our posterior or the complex distribution from which we couldn’t sample! If we manage to get such a Markov Chain, then we can sample from it to our heart’s content after a certain number of timesteps (the time it takes to reach the statioary distribution) & this would be like sampling from p!

Markov Chains help us answer questions like “How long does it take for a knight making random moves on a chessboard to return to the initial square”or even “What do the hyperlinks between web pages say about their relative popularities?” Brin & Page modeled the entire web as a Markov chain, with web pages as states and links between pages as transitions. A Random surfer bot starts on a random page (doesn’t matter which), follows the links & keeps track of the %age time spent on each page. Over time, the ratio settles and the scores give a some measure of the relative importance of the pages. MCMC follows a similar Random Walk algo. We choose a good starting point & “select” a new sample from a proposal & walk there. The “proposal” & the “selection” are the key. We need to do it right & we need to do it for a while to get to our destination — the steady state of an MC with distribution p.

How to converge to a Markov Chain with stationary distribution p:

It is a testimonial to the rapid advancement of IQ in a modern-day human that an “introduction” to MCMC explains the algorithm directly. If I were in the audience, I would expect a more simple “how do we simulate a Markov Chain” & then move onto the question of how to simulate one with a stationary distribution p & reach the steady state of that MC. Let us attempt the former:

• Start with any random sample x. This will be out current state.
• Now, generate another sample y given the current state x using a proposal distribution q. This proposal can be any simple, sampleable distribution q. This sample y is our next state. It is a conditional — q(y|x).
• Repeat. We get a sequence of samples that acts like an Markov Chain.

But the second step is tricky. How exactly do we get y given an x & some proposal distribution. And how in the world is y going to be dependent on x? Here is the thought. Let x, y be the values of the parameters of the distribution p we are trying to investigate. Let our proposal come from a Gaussian q centered on the current state! For e.g., say we want to find the unknown parameters of a distribution p given some data, i.e. we want the posterior p(θ|Data). Say p is a Gaussian & we want to calculate its parameter mu(mean) with MCMC:

• Let our starting point be any mu. This is our Current State.
• Now take a proposed sample from a Gaussian distribution q (no relation to the Gaussian p) centered around this Current State mu assuming some standard deviation(proposal_width) i.e. y ~ N(mu, proposal_width^2). This sample will be our New State. Notice how it is dependent on the previous State. Note also that the proposal need not necessarily be a Gaussian.

If we keep repeating the above, we get a Markov Chain. The State keeps moving randomly based on the initial starting point & the proposal_width. Great, this is a first step. Now how can we make this MC forget its current distribution & march towards a stationary distribution p? This is literally a trillion dollar question. Note that once we get to that state (the stationary state), we will continue to stay there… our algorithm’s job is done & we can sample away to our heart’s content without disturbing the equilibrium. So how do we get there? Turns out, this is not too complicated. Instead of accepting the next sample blindly, we need to accept it based on an equation.

Pray, what is this equation? We need to derive it from the detailed balance condition which states that π is a stationary distribution for an MC if the below condition is satisfied for all states in the distribution:

where x & y are any two states & T is the transition probability. So given a sample x, we somehow want to generate a sample y such that:

We have sample x. We also generated y from a proposal q, given the sample x. All we need to do is to accept the sample if the above condition is satisfied else simply reject it! Of course, this is easier said that done. We don’t know the denominator of p, the deadly Normalizing Constant(NC). Why don’t we replace ‘p’ by the un-normalized p divided by NC.

Since all our samples are going to come from q, the T(x, y) i.e. the Transition probability from state x to state y will be the probability of q generating the sample y given the last state is x i.e. q(y|x). Also, we may accept or reject the sample. Let the Acceptance probability be Ax->y. Then:

If our proposal is based on a symmetric distribution like the Gaussian, then q(x|y) = q(y|x). Cancelling the NCs, the final equation becomes:

The last step looks confusing but isn’t. The simple rule is — If the state y better explains the distribution p than x, we simply jump!

• If the ratio of u_p(y) to u_p(x) is ≥ 1, then Ax->y = 1. Why? The un-normalized posteriors are the product of the likelihoods & the priors. If u_p(y) > u_p(x), it means mu = y explains the data better than mu = x.
• What happens if the ratio of u_p(y) to u_p(x) is < 1. Hmm, then Ay->x = 1 & so Ax->y is simply the ratio u_p(y) / u_p(x).
• Combining, we arrive at the last step above! Thanks to this cool link.

So, MCMC proposes new parameter values & decides to accept or reject them based on a ratio that compares the posterior probability of the new state y to the current state x. If this ratio >1, then it jumps to y. If not, then it is probably a good idea to stay put in x. But we do occassionally need to wander even in wrong directions. Sometimes in life, a wrong turn might actually turn out to be a shortcut! So if ratio<1, let the acceptance probability be that ratio itself. If we keep doing this, we walk to the stationary distribution p & after a period of time, the samples we get will be from an MC with a stationary distribution p. Wow, don’t you agree this is a neat algo?

BTW, how is it different from Importance sampling? There, the samples were independent of one another leading to wastage. Here, the samples are related — the current sample depends on the last one. If we happen to be in a region of high probability, then the next sample proposed is also likely to be from that neighborhood. Our algo tries to keep us in that good neighborhood.

Also, while we we assumed a Gaussian model for p, in reality we don’t even need to know what type of distribution p is. Just define any model & let loose MCMC on it. Secondly, we used a Gaussian proposal model q which is symmetric, so q(x|y) = q(y|x). If q is non-symmetric, then the Acceptance formula slightly changes to allow for this. If q is symmetric, the algo is called the Metropolis algorithm. If, not it is called the Metropolis-Hastings.

Practical Condiderations for MCMC:

• While above, we talked of one single chain, one can also simulate an ensmble of many independent chains in parallel starting from a set of arbitrarily chosen points, sufficiently distant from each other. These chains are stochastic “walkers” — their walk designed such that more time is spent in the most probable areas of the target distribution.
• Some conditions need to be satisfied for an MC to converge to a steady state (1) Irreducible — every state can be, reached from every other state; (2) Aperiodic — it doesn’t get trapped in cycles; (3) Positive recurrent — the expected return time to a particular state is finite. For such a situation, a stationary distribution exists & is independent of the initial distribution. The MCMC we saw above satisfies all of these properties.
• MCMC methods can become “stuck” & fail to fully explore the entire probability space. Choosing a good proposal is an art. If proposals are too close to the current point, then almost every move is accepted, but it takes a long time to explore the whole space. If proposals are too far away, then almost all moves are rejected, which slows convergence. Convergence is usually determined by certain standard criterias.
• Tuning a Metropolis-Hastings for a large parameter-set situation gets unwieldy & (slightly) modified implementations like Gibbs, HMC is used.

Applications of MCMC:

• You can write a simple MCMC Metropolitan algo from scratch with a few lines of code. The visuals are good. Another nice link is this.
• In this Competition, MCMC was used to locate dark matter in data of the night sky! Dark matter cannot be seen directly but its presence can be judged by observing its effect on the distortion of galaxy shapes (ellipticity) surrounding it. The task involved guessing the position & mass of the dark matter given the observed galaxy ellipticities. So, we need to infer the parameters (positions & masses) from noisy data given some observations. A perfect fit for a Bayesian approach. We need p(θ|Data). We can’t calculate it directly but we sure can use MCMC to get enough samples & then estimate θ from these samples. See this & this.
• Code to estimate Ebola disease propagation with noisy data observation
• Using Simulated annealiing (a form of MCMC) to solve the Travelling Salesman problem. See here the answer by Willem Hendriks. Also this.
• MCMC in cryptography — Decoding prisoners messages. Wow!

Why is MCMC so powerful, yet so simple? Why does it even work? In nature, when a system is in equilibrium, then for every forward elementary process (a chemical reaction or a collision A → B), there exists a reverse process (B → A) that occurs at the same average rate. This is the core of detailed balance condition based on which our Acceptance formula was derived. So perhaps there is our clue. Who better than nature to ensure that a system moves to an equilibrium state efficiently. Cerny’s quote is bang on!.

References: The Strange Math That Predicts Anything (link). Also this & this.

Other articles I have written as part of my Back to the AI Basics series:

1. Random Variables & Probability Distributions explained — Intuition, Basic Math & application in AI/ML
2. A deeper dive into Distributions — Intuition & application in models like Stable diffusion & GMM clustering
3. Of Poets, Pragmatists & the Practitioner — Perspectives on Ridge, Lasso, Dropout regularization, viewed from a Bayesian lens
4. MCMC & the magical art of Sampling without Sampling — Story, Intuition & the gentle Math behind the greatest algorithm of the 20th century
5. Secrets of the VAE — Appreciation without the apprehension
6. The Enchanting world of GNNs — MPNNs, GCN, GAT, GTN, GraphSAGE
7. Notes 1: Rise of the EBMs: A gentle introduction to Energy Based Models & why they may replace the current LLMs
8. Notes 2: Death by Notation: Unravelling the mystery behind MLE & Expectation Maximization (often hidden behind confusing notations)
9. Notes 3: The Exponential family explained: Story behind the most famous family in the ML world
10. Lastly, my personal fav: Cosine Distance vs Dot Product vs Euclidean in vector similarity search — Why my webpages will never be ranked at the top of a search engine result (well, looks like the folks at Google are having the last laugh. This is the only article of mine that has a decent search rank!)

For the yet-to-be-authored pieces, you can check directly here.