Hello, I'm Maneshwar, and I'm building LiveReview — a blast-radius aware AI code review built for your business-critical systems. Star us to help devs discover the project, give it a try, and share your feedback to help improve the product. There is an algorithm that got invented in 1953 on a machine with less memory than the favicon of this page. It is used to forecast the weather. It is used to fit models of black hole mergers. It sits under the hood of every serious Bayesian statistics library you have ever imported. And most working developers have never heard its name. It is called Markov chain Monte Carlo, MCMC to its friends, and the reason it feels like a black box is that people usually explain it backwards. They start with detailed balance and ergodicity and stationary distributions, and by the time they get to the part that is actually clever you have closed the tab. So let's do it forwards. MCMC is two ideas bolted together, and both of them are simple enough to explain at a bar. Idea one: you can measure things by throwing stuff at them Suppose I ask you for the value of pi and take away your calculator. You could derive it. People have. It is unpleasant. Or you could draw a square of side 2R, draw a circle of radius R inside it, and start throwing darts at the square while blindfolded. A dart that lands uniformly in the square has some probability of landing inside the circle, and that probability is just the ratio of the areas. Circle is pi * R^2. Square is 4 * R^2. Ratio is pi / 4. So throw a lot of darts, count how many landed inside, multiply by four, and you have pi. That is the entire Monte Carlo method. It is named after the casino, because the people who invented it at Los Alamos were doing neutron diffusion calculations and needed a codename, and one of them had an uncle who kept borrowing money to gamble in Monte Carlo. In five lines: import random hits = sum(1 for _ in range(10_000_000) if random.random()**2 + random.random()**2 <= 1.0) print(4 * hits / 10_000_000) # 3.1417... Nobody in that snippet solved an integral. They counted. The catch, and this is the catch that the whole rest of the article exists to fix, is that random.random() gave us independent samples for free. We knew the shape we were sampling from. It was a square. Sampling uniformly from a square is easy. In every problem you actually care about, you do not know how to sample from the shape. That is the problem. Idea two: a process that only remembers where it is right now A Markov chain is a sequence of states where the next state depends only on the current one. Not on how you got here. Not on the previous forty steps. Just on here. This is called the Markov property, and the shorthand for it is that the chain is memoryless. The textbook example is weather, so let's use the textbook example. Three states: rainy, cloudy, sunny. Fixed probabilities for hopping between them. From rainy there is a 60% chance you go to cloudy and a 40% chance you stay put. From cloudy there is a 50% chance you go to sunny. And so on. Now run it. Rainy, cloudy, sunny, sunny, cloudy. Step five consulted step four and nothing else. It has no idea that the run started in the rain. Here is the property that makes this useful instead of merely cute. If the chain can reach every state and does not get permanently trapped anywhere, then the fraction of time it spends in each state converges to a fixed set of numbers. Run it a thousand steps and you might be rainy 30% of the time. Run it a million and it is still 30%. It has settled. That set of numbers is the stationary distribution, and it does not depend on where you started. Sit with that for a second, because it is the whole trick. A Markov chain, left alone, generates samples from a distribution. Not a distribution you chose. Just whatever distribution falls out of the transition rules you happened to write down. So what if you ran that backwards? What if you had a distribution you wanted, and you designed the transition rules so that its stationary distribution was exactly that thing? Then you would have a machine that spits out samples from a distribution you were never able to sample from directly. That is MCMC. That is the entire idea. Everything else is engineering. The distribution nobody can compute Time to be concrete about what "a distribution you cannot sample from" actually means, because otherwise this is all very abstract. Say you are doing Bayesian inference. You have data, you have a model with some parameters, and you want to know which parameter values are consistent with what you observed. A non-Bayesian method hands you one best fit number and a standard error. Bayesian inference hands you a whole distribution over parameters, called the posterior, which tells you both what is likely and how confident you are allowed to be about it. You get it from Bayes' theorem, which is four symbols and one enormous problem. The numerator is fine. The likelihood is "how well do these parameters explain my data", which you can evaluate. The prior is "what did I believe before I saw the data", which you wrote down yourself. The denominator is where it falls apart. It is called the evidence, and it is an integral over every possible combination of every parameter. It exists purely to make the whole thing sum to one. With one parameter, grid it at 100 points, evaluate 100 times, done before your coffee lands. With five parameters, that is ten billion evaluations. With ten parameters, which is a small model by any modern standard, you are at 10^20 and the sun has opinions about your timeline. This is the curse of dimensionality, and it is not a performance problem you can engineer around. Grid methods die here. Every time. So you are stuck holding a distribution that you want to draw samples from, and you cannot even evaluate it, because evaluating it requires a normalising constant you cannot compute. Which sounds terminal. The cancellation that saves everything Here is the move. MCMC never asks "what is the posterior probability at this point". It only ever asks "is this new point better or worse than the one I am standing on, and by how much". That is a ratio. And in that ratio, the evidence appears on the top and on the bottom. So it cancels. You never compute the impossible thing. You just arrange never to need it. The consequence is that you only need the posterior up to a constant, which is just likelihood times prior, which you can always evaluate. That is the load-bearing insight of the entire field. Metropolis-Hastings, the whole thing The oldest MCMC algorithm is Metropolis et al., 1953, later generalised by Hastings in 1970. It is short enough to hold in your head. You are standing at some parameter value. You want to take a step. Propose. Draw a candidate near where you are, usually from a Gaussian centred on your current position. Score. Compute R = p(proposed) / p(current), using the unnormalised posterior, because that is all you have and all you need. Decide. If R >= 1, the new spot is better. Move there. Always. If R < 1, the new spot is worse. Move there anyway, with probability R. That last line is the one people skip past, and it is the one that makes the algorithm work. If you only ever accepted uphill moves, you would have written a hill climber. It would sprint to the nearest peak, sit on it, and report that peak as the answer with total confidence, having never seen the rest of the distribution. The occasional deliberately bad move is what lets the chain roll down one hill and find another. It is what turns a greedy optimiser into a sampler. And the acceptance rule is not arbitrary. It is constructed so that the chain's stationary distribution is exactly the posterior you handed it. The chain visits high probability regions often and low probability regions rarely, in precisely the right proportion. So it looks like a random walk, and it kind of is, but it is a random walk with a rigged floor. Here is the whole algorithm, and I mean the whole algorithm: import numpy as np def metropolis(log_post, start, n_steps, step_size): x, lp = start, log_post(start) chain = [] for _ in range(n_steps): candidate = x + np.random.normal(0, step_size, size=np.shape(x)) lp_new = log_post(candidate) # log space, so the ratio is a subtraction and nothing overflows if np.log(np.random.rand()) < lp_new - lp: x, lp = candidate, lp_new chain.append(x) # note: appended even when we rejected return np.array(chain) Two things in there that trip people up. We work in log space. Posteriors underflow to zero fast in float64, so the ratio becomes a subtraction of log densities. np.log(rand()) < lp_new - lp is the same rule, numerically survivable. We append x even on rejection. Staying put is a real outcome. If you only recorded accepted moves you would systematically under-count the sharp peaks, which are exactly the regions where most proposals get rejected. As a flowchart: flowchart TD A[Pick a starting value] --> B[Propose a new value<br/>from a Gaussian around it] B --> C[Compute R = p_new / p_current] C --> D{R >= 1?} D -- yes, uphill --> E[Accept the move] D -- no, downhill --> F{Coin flip lands<br/>under R?} F -- yes --> E F -- no --> G[Reject, stay put<br/>and record the old value again] E --> H[Record the sample] G --> H H --> I{Still in burn-in?} I -- yes --> J[Throw this sample away] I -- no --> K[Keep it] J --> B K --> L{Enough samples?} L -- no --> B L -- yes --> M[Average the samples<br/>to get anything you want] classDef decision fill:#f4d35e,stroke:#b8991f,color:#1a1a1a classDef start fill:#e9ecef,stroke:#6c757d,color:#1a1a1a classDef good fill:#5ee6c8,stroke:#1f9c86,color:#1a1a1a classDef bad fill:#ff9a5c,stroke:#c65f22,color:#1a1a1a classDef work fill:#6ea8ff,stroke:#2f5fbf,color:#1a1a1a class D,F,I,L decision class A,M start class E,K,H good class G,J bad class B,C work The step size is the one knob, and it will bite you step_size looks like a tuning detail. It is not. It is the difference between a chain that works and a chain that lies to you. Make the proposals too narrow and almost everything gets accepted, because you barely moved. The chain shuffles across the distribution at a glacial pace and you will need millions of samples to see anything. Make them too wide and almost everything lands somewhere terrible and gets rejected. The chain stands still for hundreds of steps at a time, and your "10,000 samples" are actually about forty distinct values repeated. Both failures look like success from the outside. You get your array of 10,000 numbers either way. The folk rule for the simple random walk version is to aim for an acceptance rate around 20 to 25%, which comes from some genuinely lovely asymptotic work by Roberts, Gelman and Gilks. Print the acceptance rate. Always print the acceptance rate. This is also why nobody writes the loop above in production. Modern samplers like Stan and PyMC use Hamiltonian Monte Carlo and NUTS, which use gradients of the posterior to propose smart, distant moves instead of blind local wobbles, and tune themselves. The idea is identical. The proposal is just far less stupid. Burn-in, or: throwing away the work you paid for You have to start the chain somewhere, and your somewhere is probably wrong. If your initial guess lands in a region of terrible posterior probability, the chain will spend a while wandering out of the wilderness before it finds the part of parameter space that matters. Those early samples are real samples, they cost real compute, and they are garbage. They describe your bad guess, not the posterior. But remember the Markov property. The chain has no memory of where it started. Once it reaches the high probability region, it stays there, and nothing about its future behaviour is contaminated by the trek it took to get there. So the fix is embarrassingly blunt. Delete the first few hundred or few thousand samples. That is burn-in. It is not a hack, it is a direct consequence of memorylessness. The way you check is the trace plot on the left of that diagram. Plot parameter value against iteration. A healthy chain looks like a fuzzy horizontal caterpillar. A chain still climbing in, or drifting, or sitting flat for long stretches, is telling you something and it is not good news. In practice people run several chains from different starting points and check they all converge to the same place, which is what the R-hat statistic measures. What you do with a pile of samples Once you have samples from the posterior, the hard part is over and everything downstream is embarrassingly easy. Want the mean of a parameter? Average the samples. Want a 95% credible interval? Sort them and take the middle 95%. Want the probability that a parameter is bigger than 10? Count how many are, divide by how many you have. chain = metropolis(log_post, start=0.0, n_steps=50_000, step_size=0.8) samples = chain[5_000:] # drop burn-in samples.mean() # point estimate np.percentile(samples, [2.5, 97.5]) # 95% credible interval (samples > 10).mean() # P(theta > 10), directly Every question you can ask about a distribution is an expectation, and every expectation is approximated by an average over samples. That is the Monte Carlo half, quietly doing its job at the end. Which is where the name finally makes sense. The Markov chain gives you the samples. The Monte Carlo gives you the answers. Neither half works alone. Where the two halves meet flowchart TD P[A posterior you cannot integrate] --> Q{Can you draw<br/>independent samples?} Q -- yes --> MC[Plain Monte Carlo<br/>throw darts, average them] Q -- no --> R{Can you at least<br/>evaluate it up to a constant?} R -- no --> STUCK[You are genuinely stuck] R -- yes --> CHAIN[Build a Markov chain<br/>whose stationary distribution<br/>is that posterior] CHAIN --> WALK[Walk it for a long time] WALK --> S[Correlated samples that still<br/>average to the right answer] MC --> ANS[Means, intervals, probabilities] S --> ANS classDef decision fill:#f4d35e,stroke:#b8991f,color:#1a1a1a classDef start fill:#e9ecef,stroke:#6c757d,color:#1a1a1a classDef chip fill:#5ee6c8,stroke:#1f9c86,color:#1a1a1a classDef accel fill:#9d8cff,stroke:#5b4bcc,color:#1a1a1a classDef bad fill:#ff9a5c,stroke:#c65f22,color:#1a1a1a class Q,R decision class P start class MC,S,ANS chip class CHAIN,WALK accel class STUCK bad One honest caveat, because I do not want to oversell this. MCMC samples are correlated. Consecutive steps are near each other by construction, so 10,000 MCMC samples carry less information than 10,000 independent ones. The quantity that matters is the effective sample size, and it can be brutally smaller than the number you generated. Every decent library reports it. Look at it. MCMC also does not automatically work. It converges eventually, and "eventually" is doing real work in that sentence. A posterior with two well separated modes and a deep valley between them can keep a random walk chain trapped on one mode for longer than you are willing to wait, and the chain will look perfectly healthy the whole time. Convergence diagnostics are not paranoia, they are the job. So why does this matter to you Because the pattern generalises well past statistics. You have an object you cannot enumerate. You cannot compute its total. But you can compare two candidates cheaply, and you can take a random step. That is enough. That is the whole precondition. Simulated annealing is this. So is a large chunk of statistical physics. So is the PageRank random surfer, which is a Markov chain whose stationary distribution is the ranking. So is every diffusion model generating an image right now, iteratively walking noise toward a distribution it learned rather than one you wrote down. The 1953 paper was about hard spheres in a box. SIAM later put the Metropolis algorithm on its list of the top ten algorithms of the twentieth century, next to the FFT and QR decomposition. So next time someone asks how you do inference over a distribution you cannot compute, you have an answer. You do not compute it. You build a chain that walks through it for you, and you let the walking do the arithmetic. Your team's attention is limited, and the deluge of AI-generated code is making it harder to keep production code safe without slowing you down. I'm building LiveReview, a blast-radius aware AI code review built for your business-critical systems. Instead of presenting every diff with equal emphasis, LiveReview scores each change by blast radius — how far its impact reaches through your call graph — so you can focus attention where it actually matters. Spend code review effort where business risk is highest — not spread evenly across every diff. Try LiveReview on your codebase: