Why You Think Like a Bayesian but Were Taught Like a Frequentist

The Schoolroom Distortion

Think back to your first high school statistics class. The teacher walks to the blackboard, draws a coin, and asks a seemingly simple question: “What is the probability of flipping heads?” You and everyone else in the room answer instantly: “Fifty percent.” If pressed for an explanation, you’d probably say that if you were to flip that coin an infinite number of times, half of those flips would land on heads.

This is Frequentism, and it is the standard operating process of our educational system. It defines probability through the cold, objective lens of repeatable data. In the frequentist world, probability is understood through what would happen if we could repeat the same experiment over and over again. It is a neat, comforting framework designed for a world made of rolling dice, shuffled decks of cards, and endless time.

But there is a catch. The moment you step out of the classroom, the laboratory walls crumble. Real life rarely offers us the luxury of infinite trials. You cannot marry someone a thousand times to measure the probability of a happy marriage, nor can a company launch the exact same product campaign on Instagram a million times to test the engagement gain. In the messy arena of human existence, the frequentist definition of probability starts to feel less like a tool and more like a straitjacket.

The Reality: Built to be Bayesian

Αnd at this very point is where the grand illusion of our education becomes apparent: while we were trained to think like frequentists on paper, we seem to be naturally inclined to reason in Bayesian-like ways.

In the real world, probability isn’t about counting repetitions in the infinite future; it is about quantifying uncertainty in the present. This is the core of Bayesian statistics. Instead of demanding endless data before making a judgment, a Bayesian starts with a Prior ,  an initial belief based on past experience or intuition. Then, as new evidence arrives, they update that belief to arrive at a Posterior probability.

We don’t need a math degree to do this. The human brain often behaves like a Bayesian prediction engine. Our ancestors in the savannah didn’t have the luxury of waiting for a rustling bush to move ninety-nine times to calculate a p-value before running from a predator. They had a powerful prior: “Rustling bush equals danger.” They saw a tiny bit of new data  (maybe a flicker of yellow fur ) instantly updated their probability, and survived. It is not hard to see why a fast, Bayesian-like way of updating beliefs could have been useful for survival.

Of course, none of this means we are flawless statisticians. Decades of work by Kahneman and Tversky showed something paradoxical: when asked to reason about probabilities explicitly, with numbers on paper, humans are notoriously bad. We routinely ignore base rates and get textbook Bayesian problems wrong. But this is exactly the point. We run the Bayesian engine beautifully when we don’t think about it; we just can’t read its dashboard. Our intuition is the posterior; our conscious math is the bug.

The Supermarket and the Waiting Game

To see this evolutionary machinery in action, we don’t need to fight off tigers. We just need to look at how we navigate ordinary, adult life.

Imagine walking into a premium supermarket while traveling in a foreign country. You spot a high-end chocolate bar on a shelf, but there is no price tag. A strict frequentist approach would have a harder time here: you have zero observations for this exact item in this exact store. Theoretically the price could be two euros or two hundred, and both are equally unknown.

But you don’t panic, because your brain instantly deploys a prior. You know roughly what chocolate costs in general. You register the elegant lighting, the wooden shelving, the fact that everything else here is somehow imported. Before you have even touched the wrapper, you are carrying a probability distribution in your head, centred somewhere around four euros with a tail that would not be shocked by seven.

Notice what just happened. Nobody handed you data on this product. You built a prior out of context like the store, the country, the packaging. This is not sloppy thinking; it is the most information-efficient move available to you.

Then comes the cashier and says “That will be 14 euros.” Your prediction error spikes as the number sits far out in the tail of what you were expecting. And by the time you walk out the door, you have quietly revised something bigger than the price of one chocolate bar: your whole sense of what “premium” costs in this country. The next unlabelled shelf you meet, you will guess higher or maybe you will run away.

If pricing chocolate feels too trivial, consider the higher stakes of a first date. The date goes well, the conversation flows, you laugh at the same jokes, you float home. You believe there is roughly an 80% chance they want a second date. Before falling asleep, you send a text: “I had a wonderful time tonight.”

Now the evidence starts arriving. A reply in two minutes nudges your belief up towards 95% and you sleep like a baby. But two hours pass. Then five. The next morning your screen is still blank and this silence is deeply improbable under your original hypothesis. Without ever opening a statistics textbook, your brain grinds through a brutal round of updating. Eighty percent becomes fifty and fifty becomes thirty. Your prior stays exactly where it was, of course, because a prior is what you believed before. What is quietly dying overnight is your posterior.

A frequentist would frame the question differently. Instead of asking “what is the probability that this person likes me?”, they would ask what would happen to the response rate across repeated, comparable situations. But you don’t have a thousand attempts, of course. You have one.

The Mathematical Straitjacket

If the Bayesian approach is so natural, so deeply embedded in our cognitive wiring, why does the educational system still force-feed us frequentism? Why did most of the 20th century treat Bayes as a curiosity rather than a tool? The answer is a mix of fierce philosophical warfare and one dirty mathematical secret.

In the early 1900s, the architects of modern statistics (Ronald Fisher, and later Jerzy Neyman and Egon Pearson) were on a mission to make science objective. When they looked at Bayes’ theorem, they recoiled at one specific term: the prior. The idea of a scientist injecting personal belief into a mathematical equation felt like heresy. Fisher was deeply skeptical of it. Where do you even get a prior from, and even more importantly, how do you defend it to a referee or in a conference?

It is a fair question, and it deserves a fair answer. Notice that your prior about the chocolate bar didn’t come from nowhere. It came from every price you had ever seen, filtered through the lighting of that particular store. A prior isn’t a wish but compressed experience, written down honestly instead of smuggled in through the back door. The frequentists make assumptions too. They just doesn’t have to declare them loudly.

Fisher wasn’t persuaded. To rescue science from subjectivity, that generation built a different framework: p-values, t-tests, confidence intervals. The ambition was to make statistical conclusions depend on the observed data and a clearly specified sampling procedure, rather than on a scientist’s prior beliefs about the unknown parameters.

Bayes never quite disappeared, to be fair. Harold Jeffreys was arguing publicly with Fisher throughout the 1930s, and Alan Turing was using Bayesian reasoning at Bletchley Park to break Enigma, a work that stayed classified for decades. But these were exceptions, and the reason Bayes stayed at the margins had less to do with philosophy than most people assume.

Even a bayesian statistician who would have won the philosophical argument would immediately hit a wall, and the wall was made of arithmetic. Bayes’ theorem contains one term that, for realistic models, was often impossible to compute. That single term is what kept the entire approach out of practical reach, and it is worth looking at it directly.”

Anatomy of a Nightmare

Let’s actually look at the monster Fisher refused to fight. Bayes’ theorem itself is deceptively simple:

Here, $θ$ (theta) is everything you don’t know: all the parameters of your model, bundled together. D is your data. The left side, $P(θ|D)$is what you want: the posterior; your updated belief after seeing the data. The numerator is friendly. $P(D|θ)$is the likelihood (how well a specific guess of $θ$explains the data), and $P(θ)$ is your prior. For any single guess of $θ$, both are easy to compute.

The trouble lives entirely in the denominator. That integral has a name, the marginal likelihood, and it is normally written in the compact form P(D). A pretty logical qestion would have been: “why is it an integral at all?” or “why isn’t it simply a number you can look up?”

Here is the intuition: P(D) asks: “How probable is my data, full stop?” Not “how probable is my data if the price effect variable is 0.4,” but how probable it is overall, without committing to any particular version of reality.

The catch is that your data doesn’t have a probability in a vacuum. It only has a probability given some assumption about the world. Imagine looking at the exact same sales data under two different versions of reality. In the first, price has a strong effect on sales: whenever the price drops, sales jump. Under that version of the world, a dataset showing exactly that pattern would be completely unsurprising, so $P(D|θ)$ would be high. Now imagine a second version where price has almost no effect at all. The data hasn’t changed, but suddenly those same jumps in sales after every price cut become much harder to explain, so $P(D|θ)$ becomes low. Change θ again (perhaps promotions are doing most of the work, or seasonality is stronger than you thought) and the probability of observing that exact same dataset changes again. The data is fixed; what changes is the world you are asking the data to make sense under. Each candidate value of θ describes one such world, and each world assigns your data a different probability.

So how do you get one single number out of infinitely many stories? You do the only fair thing: you take a weighted average over all of them. For every possible θ, you ask “how likely is my data in this scenario?”. The answer is $P(D|θ)$and then you weigh it by “how plausible was this scenario to begin with?” which is $P(θ)$; at the end you sum it all up.

If θ could only take a handful of discrete values, this would be an ordinary sum:

But many parameters like “the effect of price on sales or a possible promotion” don’t come in a handful of options. They live on a continuum: the effect could be -0.50, or -0.501, or -0.5017, and so on, infinitely divisible. And when you sum over a continuum instead of a list, the sum becomes an integral. That’s all the integral is here: a sum over infinitely many scenarios.

In statistics, this move has a name: marginalization. You are “integrating out” the parameters, averaging over every version of reality you can imagine, to find how expected your data is across all of them. It is the mathematical equivalent of asking: “Across all the possible worlds I can imagine, how likely would I be to see this data?” Some worlds count more than others, because my prior says they were more plausible to begin with.$P(D)$combines all of those possibilities into one number.

So far, so elegant. As a piece of mathematical philosophy, marginalization is beautiful: one clean expression that honestly accounts for every scenario you can imagine. The problem is not writing this integral down. The problem is computing it. For a coin flip this integral is trivial but in real life we want to create real models that are far from being coin flips. Suppose ourmodel has just three unknowns: a slope, an intercept, and a noise level. Now your integral is three-dimensional:

A realistic model, say, a Marketing Mix Model (MMM) with fifteen channels, each with its own effect size, plus carry-over and saturation parameters (which are non-linear effects), can easily involve thirty or fifty dimensions. And multi-dimensional integrals do not just get “a bit harder” as dimensions grow. They explode. If you tried to approximate a 30-dimensional integral on a grid using just 10 points per dimension, you would need 10³⁰ evaluations. That is more calculations than there are stars in the observable universe, for a crude approximation of one denominator.

Worse, for most realistic combinations of priors and likelihoods, no closed-form solution exists at all. It is not that the integral is tedious; it is that, in general, no formula for its answer can be written down. Fisher’s generation had pencils, paper, and slide rules. Against a 30-dimensional integral with no analytical solution, they never stood a chance.

Bayes wasn’t held back because the theorem was wrong. It was held back because, for many realistic problems, the computation was brutal.

The Computational Prison Break

For two hundred years, that integral sat there like a locked door. Then, in the second half of the 20th century, two things happened at once: computers became fast, and a family of algorithms called Markov Chain Monte Carlo (MCMC) revealed something remarkable. You don’t need to solve the integral. You don’t even need to compute it. You just need to make it irrelevant.

Remember what we actually want: the posterior, $P(θ|D)$. And remember what’s blocking us: the denominator, $P(D)$. But now ask a slightly different question. Instead of “what is the probability of this specific θ?”, ask: “how much more probable is θ₁ compared to θ₂?” Write that comparison as a ratio, and watch what happens:

The $P(D)$ appears in both the top and the bottom. It is the same constant, the same impossible 30-dimensional integral, in both places. So it cancels completely. What survives is:

Let’s take a look at what’s left: only likelihoods and priors which are the two quantities we said from the beginning that are easy to compute for any single guess. The monster integral didn’t get slain. It just got divided by itself.

Looking closer at the above equation we cannot help noticing that a ratio only tells you which of two points is better. How do you turn that into a full picture of the posterior? This is the part where the “walking” begins.

Imagine the space of all possible parameter values as a landscape: a vast, dark, 30-dimensional terrain where altitude represents posterior probability. Peaks are parameter combinations that explain your data well; valleys are combinations that don’t. You cannot see this landscape, and you can never compute the total volume under it, since that would be the integral. But standing at any single point, you can always measure your local altitude relative to a neighbouring point, using the above ratio-formula.

So the algorithm becomes an explorer with a simple rule. From wherever you stand, propose a step in some direction and compare altitudes. If the new point is higher, go there. If it is lower, sometimes go there anyway, with a probability equal to that ratio (if the proposed point is half as probable, you accept the move half the time). This occasional downhill wandering is not sloppiness. It is what keeps the explorer from getting trapped on the first small hill it happens to find. Finally, repeat, for thousands of time.

And here is the beautiful part, the mathematical guarantee that makes all of this work: after enough steps, the amount of time the explorer spends in each region becomes proportional to that region’s probability. The walker naturally lingers on high peaks and rarely visits the lowlands.If you keep a diary of everywhere it has been, that diary becomes a portrait of the posterior. Not solved but traced (remember that word as it will show up again when we finally write some code). Like a tourist without a map who, after a week of wandering a city, has drawn the map with their own footsteps: dense scribbles in the lively squares, faint lines in the empty alleys.

Why Random Walking Is Not Enough

At this point, it looks like we have won. The impossible integral is gone, the explorer can wander through the posterior using nothing more than likelihoods and priors, and after enough steps its footprints reconstruct the distribution we could never calculate directly. For small problems, that really is almost the end of the story.

But there is a catch. The word enough is doing an enormous amount of work in that sentence. In two or three dimensions, our simple explorer can move around reasonably well. Give it thirty parameters, however, and “enough steps” can become a painfully large number. The algorithm might still be correct but at the same time is becomes painfully inefficient.

The reason is one of the most counter-intuitive facts in high-dimensional geometry.

Your first instinct might be that the sampler should simply climb towards the peak of the posterior and stay there. But that would be optimization, not sampling. The goal of MCMC is not to find the single most probable point; it is to explore the whole posterior in the right proportions.

This distinction becomes surprisingly important in high dimensions. The very top of the posterior may have the highest probability density, but it occupies a tiny amount of space. Move a little farther away and each individual point becomes less probable, but there are vastly more points available. In enough dimensions, that extra space more than compensates for the lower density. As a result, most of the posterior probability is not concentrated at the peak, but in a region surrounding it.

This region is called the typical set. It is where a correctly working sampler should spend most of its time. The mode tells you where the posterior is highest; the typical set tells you where most of its probability actually lives. And as the number of dimensions grows, efficiently moving through that region becomes much harder.

And this is where our random walker runs into trouble. It needs to explore this enormous typical set, but it has no idea which direction to go. At every step, it simply proposes a nearby point and uses the probability ratio we saw earlier to decide whether to move there. This simple strategy has a name: the Metropolis algorithm.

In a few dimensions, blindly trying directions can work reasonably well. In thirty dimensions, it becomes painfully inefficient. Make the steps large enough to cross the posterior quickly, and many proposals will jump into regions of much lower probability and get rejected. Make the steps small enough to be accepted frequently, and the walker barely moves. You are left with an unpleasant choice: reject constantly or crawl.

The important point is that Metropolis has not suddenly become wrong. Given enough time, it still targets the correct posterior. The problem is that in high dimensions, enough time can become an awfully long time.

Giving the Explorer Momentum

The fix, and the reason modern Bayesian modelling actually works in practice, is to stop letting the explorer guess.

Notice what the random walker is throwing away. At any point in the landscape, you don’t only know your altitude. You can also compute the slope, because the likelihood and the prior are ordinary functions you can differentiate. It’s like the walker having a compass in its pocket and never looking at it.

Hamiltonian Monte Carlo uses that compass. Instead of proposing a blind step, it treats the posterior as a physical surface, flips it upside down so that high-probability regions become valleys rather than peaks, and rolls a frictionless ball across it. You give the ball a random shove, and then simply let it move the way an object moves on a curved surface: accelerating downhill, decelerating uphill, curving as the terrain curves. After following that trajectory for a while, you stop, and wherever the ball has ended up becomes your proposed new point.

The consequences are dramatic. A random step in thirty dimensions is almost always a bad step. A trajectory that uses the geometry of the surface can move through high-probability regions much more efficiently, so it can travel a long distance and still propose somewhere the sampler is happy to accept. Instead of a thousand rejected nudges, you get one long, informed move.

The momentum matters too. Without it, our ball would simply slide downhill towards the bottom of the valley — the single most probable point. But a ball moving on a frictionless surface does not settle there. It carries momentum, passes through the bottom, climbs back uphill, and keeps moving. That is exactly what we want from a sampler: not to find the single best point and stop, but to keep exploring the region where the posterior probability lives.

Under the hood, this is Hamiltonian mechanics in the literal sense: define a potential energy from the log posterior, add a kinetic energy term for the momentum, and evolve the system through Hamilton’s equations. The core mathematical machinery survived the trip from 19th-century physics into a Python library, which is either a coincidence or a hint about how deep the connection between probability and geometry runs.

The minus sign is doing something important here. A parameter combination with a large$P(D|θ) P(θ)$, has a low potential energy. In other words, the minus sign turns the peaks of our probability landscape into valleys of potential energy. We have taken a statistical problem and deliberately rewritten it to look like a problem from physics.

We then invent a momentum for our parameters and give it kinetic energy. The parameters do not literally have mass or velocity; these are artificial variables introduced purely to make the sampler move. Add the potential and kinetic energies together, and you have a Hamiltonian system. From there, the algorithm follows Hamilton’s equations to generate a long, informed trajectory through the posterior.

That single expression is also where the last trace of our monster disappears. P(D) is a constant, so leaving it out only shifts the whole landscape up or down by a fixed amount. The shape stays identical, and the ball rolls exactly the same way.

What NUTS Adds

One awkward question remains: how long do you let the ball roll? Stop too early and you are back to small, random-walk-sized steps. Roll too long and the trajectory curves around and comes back towards where it started, burning computation to go nowhere.

The No-U-Turn Sampler, which is the default engine inside PyMC, removes that decision from your hands. It extends the trajectory while checking whether the path has begun to double back on itself, and uses that U-turn as a signal that further travel is no longer useful. The length adapts to the local shape of the posterior at every single iteration.

This is what is happening when you call pm.sample() and get a progress bar instead of a research project. It also gives you something unexpectedly useful: when the sampler cannot follow the trajectory reliably, it reports a divergence. That is not a nuisance warning to be silenced. It is the algorithm telling you that your posterior has a shape it cannot navigate, which is usually a sign that the model itself needs rethinking.

The integral that made Bayesian inference computationally hopeless for so many realistic models never disappeared. For posterior sampling, it simply stopped being necessary. We found a way to trace the shape of the mountain without ever measuring the total volume beneath it.

The Model, in Practice

But, enough philosophy!

All of this may still sound like an elegant story about algorithms, so let’s make it concrete. The example comes from Marketing Mix Modeling: untangling how much promotions, price changes, seasonality, and advertising each contribute to sales. Marketing data is short, noisy, and full of channels that move together, which makes Bayesian methods less a philosophical preference than a practical necessity.

Here is what a small but honest Marketing Mix Model looks like in PyMC, with five ingredients: a baseline, a promo effect, a price effect, seasonality, and a media channel with diminishing returns.

import pymc as pmimport pytensor.tensor as ptimport numpy as np# week_of_year, promo_flag, price, media_spend and observed_sales# are NumPy arrays with one entry per week of historywith pm.Model() as mmm:    # 1. Baseline: sales we would get with no marketing at all    intercept = pm.Normal("intercept", mu=100, sigma=20)    # 2. Promo: we believe promotions help, so we truncate below zero    beta_promo = pm.TruncatedNormal("beta_promo", mu=0.4, sigma=0.2, lower=0)    # 3. Price: we believe higher prices hurt, so we truncate above zero    beta_price = pm.TruncatedNormal("beta_price", mu=-0.5, sigma=0.2, upper=0)    # 4. Seasonality: a yearly wave, letting the data learn its amplitude    amplitude = pm.HalfNormal("amplitude", sigma=10)    seasonality = amplitude * pt.cos(2 * np.pi * week_of_year / 52)    # 5. Media with diminishing returns: the first 100k of ad spend    #    works harder than the next 100k (saturation via a log curve)    beta_media = pm.HalfNormal("beta_media", sigma=0.3)    media_effect = beta_media * np.log1p(media_spend)    # Put it all together    mu = (intercept          + beta_promo * promo_flag          + beta_price * price          + seasonality          + media_effect)    sigma = pm.HalfNormal("sigma", sigma=10)    sales = pm.Normal("sales", mu=mu, sigma=sigma, observed=observed_sales)    trace = pm.sample(2000, tune=1000)

And there is that word again: trace.

Remember our explorer’s diary? That is essentially what pm.sample() is giving us. It does not solve the posterior into a neat formula. It sends the sampler through the parameter space and records the values it visits. Those recorded draws are our trace: thousands of plausible combinations of the intercept, promo effect, price effect, seasonality, media effect, and noise, visited in the proportions dictated by the posterior.

That is why the output is so much richer than a single “best” coefficient. We did not ask PyMC to find the best answer. We asked it to trace the distribution of plausible answers.

The lower=0 on promo and the upper=0 on price are doing real philosophical work: you are telling the model, mathematically, “I refuse to consider a universe where promotions hurt sales or higher prices help them.” Fisher would have called that unscientific bias. A marketer would call it Tuesday.

The np.log1p on media spend encodes a different kind of prior knowledge, a structural one: the belief that advertising saturates. Doubling your ad budget does not double its effect, and the logarithm bakes that diminishing-returns shape directly into the model, before a single data point is seen.

And the numbers themselves are not invented either. A prior like mu=0.4, sigma=0.2 on the promo effect comes from somewhere specific: last year’s version of the same model, a geo experiment, or the elicited judgement of people who have been running these promotions for a decade. The sigma is where you declare how much you trust that source. Set it tight and you are saying the history is reliable and the data will have to fight hard to move you. Set it loose and you are letting two years of noisy weekly data overrule it. That negotiation is the actual craft of Bayesian modelling: deciding how much weight your existing knowledge should carry before the new data arrives.

This is the answer to Fisher’s objection, incidentally. He asked where a prior comes from, as though the alternative were to have none. But leaving prior knowledge out of the model does not make that knowledge disappear. Nobody in the room really believes that a price effect of +50 is just as reasonable a possibility as one of −0.5; a Bayesian model simply gives you a formal way to say so. The choice was never between assumptions and no assumptions. It was between which assumptions you make explicit and which ones you leave outside the model.

And then pm.samplequietly does what would have been computationally unthinkable for most of Bayesian history. It walks its thousands of MCMC steps around this six-dimensional space, never once computing P(D), and hands you back not point estimates but full distributions: how big each effect plausibly is, and how sure you are entitled to be about it.

That last part is the real gift. A frequentist regression can quantify uncertainty too, of course; it might give you a coefficient surrounded by a very wide confidence interval. But that interval is not a probability distribution over the coefficient itself. The Bayesian posterior is. If most of its mass lies between 0.2 and 0.6, you can say something much closer to what you actually mean: given the model, the data, and the prior, the promo effect is probably somewhere in that range. An honest “I’m not sure” is built directly into the answer.

And that changes the conversation you can have with the business. A p-value can tell you whether the data are difficult to reconcile with a particular null hypothesis, but that is rarely the question sitting in a CMO’s head. They want to know: what is the probability this channel returns more than it costs? With posterior samples, that question becomes almost embarrassingly simple. Count the samples where ROI exceeds one and divide. If 87% of them do, you can say: “Given what we know and the model we have built, there is an 87% probability that this channel has an ROI above one.” That is a language a business can actually make decisions with.

References

The Metropolis Algorithm (1953)

• Title: “Equation of State Calculations by Fast Computing Machines” (Published in: The Journal of Chemical Physics)

The “Base Rate & Lawyer-Engineer” Study (1973)

• Title: “On the Psychology of Prediction” (Published in Psychological Review)

The “Availability Heuristic” Study (1973)

• Title: “Availability: A Heuristic for Judging Frequency and Probability” (Published in Cognitive Psychology)

The Overview of Heuristics and Biases (1974)

• Title: “Judgment under Uncertainty: Heuristics and Biases” (Published in Science)

The “Framing & Asian Disease” Study (1981)

• Title: “The Framing of Decisions and the Psychology of Choice” (Published in Science)

Hamiltonian Monte Carlo (1987)

• Title: “Hybrid Monte Carlo” (Published in: Physics Letters B)

Hamiltonian dynamics in MCMC (2011)

• Title: “MCMC using Hamiltonian dynamics“ (Handbook of Markov Chain Monte Carlo)