After you read this essay, you will be able to use Bayesian reasoning on your own problem — and read the advanced books without fear.
You have seen the Bayes theorem in school. You may have opened a professional toolbox such as PyMC and bounced off it. You suspect that you are missing something, but you cannot name it. This essay names it.
School presents the Bayes update rule like this:
\[p(A|B) = \frac{p(B|A)\,p(A)}{p(B)}\]
The equation is correct. It is easy to prove. And yet, when you face a real problem, it does not tell you what to do. Even professional data scientists get stuck here. The equation fails us in practice for two reasons.
Reason 1: the notation hides where each term comes from. All four terms are written with the same letter \(p\), as if they were the same kind of thing. They are not:
In paper-and-pencil mathematics, symbols are shortened to save hand movements. Computer science learned long ago that such savings are counter-productive. We will un-shorten them.
Reason 2 — the important one: the equation says nothing about the model. School trains us to expect that the given equation is all we need: the teacher said \(F = ma\), and every exercise follows from it. Bayes does not work like that. The Bayes rule tells you nothing about how the world works. Before the rule can tell you anything, you must create a model of the world yourself — and construct its likelihood function. Only then can Bayes tell you how good each version of your model is. No amount of staring at the update equation produces the likelihood function; it is not in there.
This is the missing link. Books discuss priors at length — informative, uninformative, how to defend yours — but the update is useless until you understand where the likelihood function comes from.
And the answer is: you invent it.
To keep the concepts tied to reality, we start with a plant. A plant is word, used to describe any object or process with inputs and outputs — a coin, a sensor, a market, a patient. The plant produces data \(D\), which we can observe.
We describe the plant with a model \(M\).
There can be more than one good model for the same plant; this topic alone deserves a book. However, let’s just take this for granted now, and press on.
Let’s say that the model (that you have just invented) takes a hidden (latent) parameter \(\phi\), so that, loosely, \(D = M(\phi)\). We cannot measure \(\phi\) directly. We can only observe \(D\) and hold beliefs about \(\phi\).
One more grounding point: in Bayesian reasoning, nothing on your desk is random. The likelihood values, the prior, the posterior — all are deterministic computations. The only place where randomness lives is inside the plant, and your model of the plant has to describe how this randomness produces data \(D\) as a function of \(\phi\).
It is hence very important to keep the distinction between these three things.
The three are measured on the same 0-to-1 scale, which is probably why history gave them the same letter \(p\). To keep them apart, this essay gives each its own unit:
| Unit | Concept | Example reading |
|---|---|---|
| \([R]\) | belief | 0.7 R: “I am 70% sure this value of \(\phi\) is the true one.” |
| \([\Omega]\) | probability | 0.8 Ω: “sampling forever, the plant gives this outcome 80% of the time.” |
| \([S]\) | sample fraction | 0.9 S: “in this data set, this outcome occurred 90% of the time.” |
Give each term its own letter and let the letter state the term’s job:
(We write the script letter \(\ell\), not plain \(l\), because plain \(l\) is too easy to confuse with \(1\), \(I\), and \(|\). Where \(\ell\) cannot be rendered, write \(L\) — the capital is also standard for likelihood in statistics.)
With these letters, the update rule reads:
\[b(\phi|D) \leftarrow b(\phi) \times \frac{\ell(D|\phi)}{n()}\]
Every symbol now has one meaning, and the equation itself shows the flow: posterior = prior scaled by how well this candidate \(\phi\) explains the data.
Classic texts call this term the “marginal probability” or the “evidence”. They still write it as \(p(D)\), and thereby cause maximal confusion. In fact, it needs a sum over all \(\{\phi\}\); its inputs include your beliefs, and it has none of the properties of the other \(p\)’s.
The honest, full signature is:
\[n() = n(D,\ \ell(),\ \{\phi\},\ \mu(\{\phi\})) = \sum_{\alpha \in \{\phi\}} \ell(D|\alpha)\, \mu(\alpha)\]
where \(\alpha \in \{\phi\}\) means “for each candidate \(\phi\) from the set \(\{\phi\}\), take one, temporarily name it \(\alpha\) and use that one only”
Note that the prior enters through its measure \(\mu() = b()\,/\,1\mathrm{R}\), not as raw belief \(b\) — the reason becomes visible in the unit check below.
The sum runs over a fresh symbol \(\alpha\), so it does not collide with the single \(\phi\) in the update rule. (In the classic notation, the \(\phi\) inside the sum silently means “all \(\phi\)s” while the \(\phi\) outside means “this one \(\phi\)” — a notation conflict that breaks ordinary substitution rules.)
Also, note that this function takes a function \(\ell()\) as a parameter (!). and this is where many books cause confusion again, by plugging in their custom “model of the world” too early.
\(n()\) is not a property of the world. It is a normalising constant; a mixing rule; mixing weights come computed from your model and your prior, whose only job is to make the posterior sum to 1. That is why this essay calls it the normaliser.
If the “probability” of data appearing to depend on your beliefs ever surprised you — this is why: \(n()\), previously denoted as \(p(D)\), was never a probability of anything.
The full, substitution-safe update rule is:
\[b(\phi|D) \leftarrow b(\phi) \times \frac{\ell(D|\phi)}{\sum_{\alpha \in \{\phi\}} \ell(D|\alpha)\, \mu(\alpha)}\]
We gave belief, probability, and sample fraction their own units. Units obey two rules, the same rules as in physics:
Both symbols of the prior already obey these rules. \(\mu() = b()\,/\,1\mathrm{R}\) is itself a rule-1 operation: \([R]/[R]\) cancels, in the open, and out comes a plain number. Dividing by \(1\,\mathrm{R}\) means “divide by your total belief budget” — so \(\mu\) is the share of belief on a candidate, and shares always sum to 1.
Now push the units through the whole update and watch every step balance.
The normaliser is an expectation. Read the formula as what it is: a weighted average of the likelihood, taken under the measure \(\mu\):
\[n() = \sum_{\alpha \in \{\phi\}} \ell(D|\alpha)\,\mu(\alpha) = \mathbb{E}_{\mu}\!\left[\,\ell\,\right]\]
In words: the probability you expected the data to have, before you saw it. Expectations have carried units correctly for three hundred years, always the same way: the averaged quantity keeps its unit, the weights are dimensionless. Compute the average weight of a population, \(\bar{w} = \sum_i w_i f_i\): the \(w_i\) are in kilograms, the population fractions \(f_i\) are plain numbers, and the average of kilograms is kilograms. Same here: an average of \([\Omega]\)-quantities is an \([\Omega]\)-quantity. Each term of the sum carries \([\Omega \cdot 1] = [\Omega]\), all terms match, rule 2 permits the sum, and \(n()\) carries \([\Omega]\).
The update balances. The Bayes factor \(\ell(D|\phi)\,/\,n()\) is \([\Omega]/[\Omega]\) — a pure number: “how much better than expected did this candidate predict the data”. The posterior is \(b(\phi) \times \text{(pure number)}\), which carries \([R]\). Belief in, belief out. No step breaks; no symbol carries two meanings.
(If you had written raw \(b\) inside the normaliser’s sum instead of \(\mu\), the units would refuse: each term would carry \([\Omega \cdot R]\), and the posterior would come out with no unit at all. That refusal is the unit system working — it rejects the mis-parse that treats an expectation’s weights as if they were the quantity itself. The \(\mu\) symbol makes the wrong version unwritable.)
The unit bookkeeping now derives three facts that we earlier only asserted:
Rule 2 also rejects a whole family of classic mistakes before they produce a wrong number:
The toolbox enforces all of this (bayesian_reasoning/units.py,
built on sympy.physics.units). Beliefs are tagged
Q(values, R), likelihoods come back tagged with
OMEGA, and a forbidden operation raises a
UnitError that states the rule and hints at the correct
next step. The conversion from belief to measure is one explicit,
visible call: prior.measure(), implementing literally \(\mu() = b()\,/\,1\mathrm{R}\).
Better asked: what should I believe about the bias of this coin?
The plant is a coin. Our model \(M\): each toss gives 1 with probability \(\phi\) and 0 with probability \(1-\phi\). A fair coin has \(\phi = 0.5\); if there was a \(\phi = 0.1\) it would mean a coin that is heavily biased towards zeros. The model ignores the coin landing on its edge — models are allowed to approximate.
From the model description, directly:
\[\ell(D{\equiv}1\,|\,\phi) = \phi \qquad \ell(D{\equiv}0\,|\,\phi) = 1-\phi\]
For several independent tosses, the likelihoods multiply. This function did not come from the Bayes rule and it did not come from the world. We built it from our model. This step — not the update — is where the real work of Bayesian reasoning happens.
In code (bayesian_reasoning/core.py):
def coin_likelihood(observations, phi_grid):
obs = numpy.asarray(observations).reshape(-1, 1) # observations in rows
phi = numpy.asarray(phi_grid).reshape(1, -1) # candidates in columns
per_toss = numpy.where(obs == 1, phi, 1.0 - phi)
return Q(per_toss.prod(axis=0), OMEGA) # probabilities carry ΩWe consider three candidate values, \(\{\phi\} = [0.2, 0.5, 0.8]\), and start with an uninformative prior, \(b(\{\phi\}) = [\tfrac13, \tfrac13, \tfrac13]\).
We toss twice and get \(D = [0, 0]\). Then:
prior = Q([1/3, 1/3, 1/3], R) # beliefs carry R
likelihood = coin_likelihood([0, 0], phi_grid) # [0.64, 0.25, 0.04] Ω
posterior = belief_update(prior, likelihood) # [0.69, 0.27, 0.04] R
Two zeros in a row shift belief strongly towards the low-bias candidate — but \(\phi = 0.5\) keeps a healthy 27%. Two tosses are weak evidence, and the posterior says so, quantitatively.
Run the same data through four different priors:
A prior of “I am 80% sure the coin is fair” survives two zeros almost intact. A prior of “biased towards ones” gets pulled hard towards the low candidates. No prior is punished for existing; each is revised by exactly the strength of the evidence.
Feed tosses one at a time, using each posterior as the next prior:
By thirty tosses of a \(\phi = 0.2\) coin, the belief in \(\phi = 0.2\) is close to 1, and the wrong candidates are extinguished. Note the dips: single contrary observations move belief the wrong way, and further data recovers it. That is not a flaw; that is what honest reasoning under uncertainty looks like.
Run it yourself:
pixi run python examples/coin_toss.py
pixi run python examples/make_figures.pyDoes this scale beyond coins? Yes. Take the proposition “taller people are heavier” and ask: what should we believe about the proportionality coefficient?
Create the model. Input: height \(h\). Output: weight \(w\). Two hidden parameters: the coefficient \(\phi_{prop}\) and an uncertainty width \(\phi_\sigma\). The model:
\[w = h \cdot \phi_{prop} + N(0, \phi_\sigma)\]
Construct the likelihood function. We want the likelihood of an observed weight to fall as it gets further from the model’s noise-free prediction. With \(\Delta = D_w - h \cdot \phi_{prop}\), a function with the right properties is the familiar Gaussian shape:
\[\ell(D_w\,|\,\phi_{prop}, \phi_\sigma) = \exp\!\left(-\frac{\Delta^2}{2\phi_\sigma^2}\right)\]
Large errors are heavily penalised; a large \(\phi_\sigma\) forgives distance. Every value lands in \([0, 1]\). That is all the update rule requires — the values across candidates need not sum to anything in particular, because the normaliser handles the scaling.
A constructed likelihood can fail quietly — run it and look. Run the update with this likelihood over a 2D candidate grid, and it recovers \(\phi_{prop}\) nicely — but the belief in \(\phi_\sigma\) runs to the largest candidate on the grid, every time. Why: a bigger \(\phi_\sigma\) forgives every observation more, and this function never charges for the forgiveness. As built, \(\phi_\sigma\) is unlearnable. The charge is the missing prefactor of the full Gaussian density:
\[\ell(D_w\,|\,\phi_{prop}, \phi_\sigma) = \frac{1}{\phi_\sigma\sqrt{2\pi}} \exp\!\left(-\frac{\Delta^2}{2\phi_\sigma^2}\right)\]
A wide \(\phi_\sigma\) now pays rent: it spreads its probability thinly, so a data point near the line rewards it less than a tight \(\phi_\sigma\) does. With the prefactor in place, both parameters become learnable:
From 25 simulated people (true \(\phi_{prop} = 45\), \(\phi_\sigma = 5\)), the highest-belief
candidate lands at \(\phi_{prop} =
44.5\), \(\phi_\sigma = 5.0\).
Run it: pixi run python examples/linear_regression.py. Note
that the toolbox needed zero changes for this — a new
model brings its own likelihood function, and belief_update
does not care what the candidates mean.
This is a likelihood function, not the likelihood function. If your plant’s noise is not Gaussian, build a different one. The freedom to choose the model is not a bug in Bayesian reasoning; it is the whole point. But as the \(\phi_\sigma\) trap shows: after you construct a likelihood, test that it can actually learn each parameter you care about.
Good — that is the correct prior. Contact me on LinkedIn and tell me why I am wrong.