LUNAROPS · OPERATIONAL UPLINK 100% UPTIME 1,247d POSTS 893 JEFF.MOON@LUNAROPS.DEV UTC --:--:--

Bayesian Statistics for Engineers

17 min read Mathematics for Engineers
bayesian-statisticsprobabilitystatisticsmcmcab-testinginference
Contents

The single most useful thing an engineer can learn about statistics is that “Bayesian” and “frequentist” are not rival religions arguing about the same question — they are precise answers to different questions, and most of the heat in the debate comes from people confusing which question they actually need answered. A frequentist computes the probability of observing your data given a hypothesis: P(data | hypothesis). A Bayesian computes the probability of a hypothesis given your data: P(hypothesis | data). Those are not the same number, they are not even the same kind of number, and the famous failures of statistical reasoning in engineering — misread p-values, peeking at A/B tests until they go significant, confidence intervals interpreted as if they were credible intervals — almost all come from quietly substituting one for the other. The argument of this post is narrow and practical: learn the Bayesian machinery not because frequentism is wrong (it is not, and treating “frequentist” as an insult marks you as someone who has not done the reading), but because Bayesian methods directly produce the quantity engineers usually want — a probability distribution over the thing they are trying to decide — and because in a handful of common engineering situations they are simply the better tool.

What follows is the working subset: Bayes’ theorem as belief updating, priors and what they cost you, conjugate distributions and why they were a big deal before computers, a fully worked A/B test with real conversion counts, MCMC at the level you need to actually run PyMC or Stan without lying to yourself about convergence, and an honest accounting of where the whole approach is more trouble than it is worth.

Bayes’ Theorem Is a Rule for Updating Beliefs

Bayes’ theorem is a one-line consequence of the definition of conditional probability, and writing it down is the easy part. The hard part is taking seriously what each term means.

                    P(D | H) * P(H)
   P(H | D)  =  -----------------------
                        P(D)

   posterior  =  likelihood * prior / evidence

Read it as a verb. You start with a prior P(H): your belief about the hypothesis before seeing data. Data D arrives. The likelihood P(D | H) tells you how well each candidate hypothesis explains that data. You multiply them, normalize, and out comes the posterior P(H | D): your updated belief. The denominator P(D), the evidence or marginal likelihood, is just the integral of the numerator over every possible hypothesis — the constant that makes the posterior integrate to one:

   P(D) = integral over theta of  P(D | theta) * P(theta) d-theta

That integral is the whole story of why Bayesian computation is hard. For most interesting models it has no closed form, and approximating it well is what MCMC exists to do. But the structural insight is cleaner than the arithmetic: the posterior is proportional to likelihood times prior. You can often ignore the evidence term entirely while you reason about the shape of the answer, then restore normalization at the end.

The mental model is sequential and cumulative. Today’s posterior is tomorrow’s prior. Process the data one batch at a time or all at once — for a fixed model the answer is identical, because multiplication is associative. This is not a metaphor; it is an algebraic property, and it is exactly why Bayesian sequential testing behaves so differently from fixed-horizon frequentist testing, which we will come back to.

   prior belief            evidence            updated belief
   ------------    x    --------------   =    ----------------
    P(theta)            P(data|theta)            P(theta|data)
   (what you knew)     (what data says)        (what you now know)
                              |
                              v
                   normalize by P(data)

Priors: Informative, Weakly-Informative, and Flat

The prior is the part newcomers fear and experts respect. The fear is that you are “putting your thumb on the scale.” The respect is that you are always putting your thumb on the scale — a frequentist analysis encodes assumptions too (the model, the stopping rule, the choice of test), it just hides them. Bayesian analysis forces the assumption into the open where it can be argued about.

Priors come in three practical flavors:

  • Informative priors encode real knowledge. If you have shipped twelve previous versions of a feature and conversion always lands between 8% and 12%, a prior concentrated there is not cheating; it is using information you would be negligent to discard. This is where Bayes shines on small data — the prior does real work when the data is too thin to speak for itself.
  • Weakly-informative priors rule out the absurd without committing to much. A prior that says “the effect is probably within a few orders of magnitude of zero” stabilizes estimation, prevents the sampler from wandering into nonsense, and barely moves the answer once you have a moderate amount of data. This is the modern default, recommended by the Stan team for almost everything.
  • Flat (uninformative) priors try to “let the data speak.” They are seductive and frequently a trap: flat in one parameterization is not flat in another, and a uniform prior over an unbounded parameter is improper (does not integrate to one), which can silently produce an improper posterior. “Uninformative” is often informative in ways you did not intend.

The honest danger is real: with little data and a strong prior, the prior drives the result. The discipline is to (1) report the prior explicitly, (2) run a sensitivity analysis by re-fitting under a few reasonable alternative priors, and (3) confirm that with enough data the prior washes out — which, for any well-specified model, it provably does. If your conclusions flip when you nudge a defensible prior, you do not have a conclusion; you have a prior wearing a conclusion’s clothes.

Conjugate Distributions: Closed-Form Posteriors

Before cheap computation, the only Bayesian problems you could solve by hand were the ones where the math closed up neatly. A prior is conjugate to a likelihood if the posterior comes out in the same family as the prior. You update by adjusting parameters instead of computing integrals — the posterior is exact, instant, and analytic.

The three you will actually meet:

Likelihood Conjugate prior Posterior Update rule (intuition)
Binomial (successes/trials) Beta(a, b) Beta(a + s, b + f) add successes to a, failures to b
Poisson (event counts) Gamma(k, theta) Gamma(k + sum x, theta / (1 + n*theta)) accumulate counts and exposure
Normal (known variance) Normal(mu0, tau0^2) Normal(weighted mean, smaller variance) precision-weighted average of prior and data

The Beta-Binomial pair is the workhorse for anything that is a rate or a probability — conversion rates, click-through, defect rates, the probability a coin is fair. The Beta(a, b) distribution lives on [0, 1] and is wonderfully interpretable: a and b act like pseudo-counts of prior successes and failures. Beta(1, 1) is uniform — one imaginary success, one imaginary failure, total ignorance. Observe s real successes and f real failures, and the posterior is just Beta(1 + s, 1 + f). No integral. That is the entire update.

Conjugacy still matters for three reasons. Pedagogically it lets you see Bayes updating without a sampler in the way. Computationally it gives exact answers with no convergence to worry about, which is why conjugate sub-blocks appear inside larger Gibbs samplers. And practically, an enormous fraction of real engineering questions — is this rate higher than that rate? — are exactly Beta-Binomial problems you can solve in three lines of code. The information-theoretic reading of the prior as encoded knowledge connects to the broader idea of entropy as quantified surprise: a sharp prior is a low-entropy belief, a flat one is high-entropy.

A Worked Bayesian A/B Test

Here is the canonical case, with real numbers. You run an experiment. Variant A converts 90 of 1000 visitors; variant B converts 110 of 1000. Did B actually win, and by how much, and how sure are you?

Put a uniform Beta(1, 1) prior on each conversion rate (total ignorance, so the data dominates). Apply the conjugate update:

   rate_A ~ Beta(1 + 90,  1 + 910)  = Beta(91,  911)
   rate_B ~ Beta(1 + 110, 1 + 890)  = Beta(111, 891)

The posterior mean of a Beta(a, b) is a / (a + b):

   E[rate_A] = 91  / 1002 = 0.0908   (9.08%)
   E[rate_B] = 111 / 1002 = 0.1108   (11.08%)

So far this looks like the obvious point estimates. The Bayesian payoff is that we have entire distributions, so we can ask the question a product manager actually cares about: what is the probability that B is genuinely better than A? That is P(rate_B > rate_A), and there is no clean closed form for the difference of two Betas — but it is trivial by Monte Carlo. Draw samples from each posterior and count.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import numpy as np
from scipy.stats import beta

rng = np.random.default_rng(0)

# Posterior parameters after the conjugate update
a_A, b_A = 91, 911     # Beta(1+90, 1+910)
a_B, b_B = 111, 891    # Beta(1+110, 1+890)

# Draw from each posterior
N = 200_000
samples_A = beta.rvs(a_A, b_A, size=N, random_state=rng)
samples_B = beta.rvs(a_B, b_B, size=N, random_state=rng)

# Probability B beats A
p_B_better = np.mean(samples_B > samples_A)

# Posterior distribution of the lift (relative improvement)
lift = (samples_B - samples_A) / samples_A

print(f"P(B > A)            = {p_B_better:.4f}")
print(f"E[absolute lift]    = {np.mean(samples_B - samples_A):.4f}")
print(f"E[relative lift]    = {np.mean(lift):.3f}")
print(f"95% credible lift   = "
      f"[{np.percentile(lift, 2.5):.3f}, {np.percentile(lift, 97.5):.3f}]")

Running this gives roughly:

   P(B > A)            = 0.9242
   E[absolute lift]    = 0.0200
   E[relative lift]    = 0.232
   95% credible lift   = [-0.030, 0.560]

Read those outputs carefully, because they are exactly what people think a p-value gives them and never does. P(B > A) ≈ 0.92 means: given the data and the prior, there is about a 92% probability that B’s true conversion rate exceeds A’s. That is a probability about the world, the thing the stakeholder asked for. The 95% credible interval on relative lift means there is a 95% probability the true lift lies in that range — and notice it includes zero at the low end, which honestly says “probably better, but we cannot rule out no improvement.” A decision-maker can act on these numbers directly: if shipping B is cheap and reversible, 92% is plenty; if it is expensive and permanent, maybe gather more data until the credible interval clears zero.

Contrast with the frequentist two-proportion test. A z-test on 90/1000 versus 110/1000 yields a p-value around 0.13 — not significant at the conventional 0.05 threshold. The p-value answers: “if A and B had identical true rates, how often would we see a gap at least this large?” About 13% of the time. That is a statement about hypothetical repeated experiments under a null hypothesis, not a probability that B is better. Both analyses are correct; they answer different questions. The Bayesian “92% chance B wins” and the frequentist “not significant at 0.05” are not in contradiction — they are different sentences.

Why Peeking Is Fine for Bayesians and Fatal for Frequentists

Here is where the practical difference bites hardest. In a fixed-horizon frequentist test, the p-value’s guarantee — at most a 5% false-positive rate — holds only if you decide the sample size in advance and look once. If you peek every day and stop the moment p < 0.05, you are running many correlated tests, and your real Type I error inflates toward 30%, 50%, higher. This is the single most common way A/B testing platforms lie to their users. The frequentist fixes (sequential testing, alpha-spending, group-sequential designs) exist and work, but they require you to plan your peeking in advance.

The Bayesian posterior carries no such curse. The posterior P(rate_B > rate_A) is a valid statement about your current belief no matter when you compute it, because it conditions only on the data in hand and never references a stopping rule or a sampling distribution of hypothetical experiments. You can update after every visitor and stop whenever the posterior probability or the expected loss crosses a threshold you set from business logic. There is no multiple-comparisons penalty for looking, because you are not comparing against a null sampling distribution at all. This optional-stopping property is, for online experimentation, the most underrated practical advantage of the entire Bayesian framework. (It is not literally free — naive stopping rules can still bias your estimate of effect size, and you should use a decision rule based on expected loss rather than chasing a probability threshold — but the catastrophic Type I inflation simply does not occur.)

MCMC at a Working Level

Conjugacy is a luxury. The moment your model has multiple interacting parameters, hierarchical structure, or a non-conjugate likelihood, the evidence integral P(D) becomes intractable and there is no closed-form posterior. Markov Chain Monte Carlo is the workaround: instead of computing the posterior, you build a random walk whose stationary distribution is the posterior, run it for a long time, and treat the samples it visits as draws from the posterior. You never compute the normalizing constant — and that is the magic, because MCMC only needs the unnormalized posterior (likelihood times prior), which you can always evaluate.

The foundational algorithm is Metropolis-Hastings, and its logic fits in a paragraph. From your current parameter value theta, propose a nearby value theta’. Compute the ratio of unnormalized posterior densities r = p(theta’) / p(theta) (the unknown normalizer cancels in the ratio — this is the whole trick). If r >= 1, the proposal is more probable, so accept it. If r < 1, accept it with probability r; otherwise stay put. Repeat millions of times. The chain spends time in each region of parameter space in proportion to its posterior probability, so a histogram of the visited values approximates the posterior. Simple, correct, and often painfully slow, because in high dimensions a random-walk proposal either takes tiny steps (and crawls) or large ones (and gets rejected).

Modern tools — Stan, PyMC, NumPyro — use Hamiltonian Monte Carlo (HMC) and its self-tuning variant NUTS (the No-U-Turn Sampler). The idea borrows from physics: treat the negative log-posterior as a potential energy landscape, give the parameter a random momentum, and simulate it rolling along the surface. Because the trajectory follows the gradient of the posterior, HMC proposes distant, high-acceptance moves and explores correlated, high-dimensional posteriors vastly more efficiently than random-walk Metropolis. NUTS removes the need to hand-tune trajectory length. You rarely write any of this; you specify a model and the sampler handles it:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import pymc as pm

# Same A/B test, now as a sampled model (overkill here, but shows the pattern)
with pm.Model() as model:
    rate_A = pm.Beta("rate_A", alpha=1, beta=1)
    rate_B = pm.Beta("rate_B", alpha=1, beta=1)

    pm.Binomial("obs_A", n=1000, p=rate_A, observed=90)
    pm.Binomial("obs_B", n=1000, p=rate_B, observed=110)

    # Deterministic quantity we care about
    pm.Deterministic("diff", rate_B - rate_A)

    idata = pm.sample(2000, tune=1000, chains=4, target_accept=0.9)

print(pm.summary(idata, var_names=["rate_A", "rate_B", "diff"]))
print("P(B > A) =", (idata.posterior["diff"] > 0).mean().item())

The summary table reports the posterior mean and credible interval — and, critically, the diagnostics. You do not trust MCMC output until you check that the sampler converged, and that check is non-negotiable:

  • R-hat (Gelman-Rubin): compares variance within chains to variance between chains. It should be very close to 1.00. Anything above ~1.01 means your chains have not mixed — they disagree about where the posterior is — and your results are garbage.
  • Effective sample size (ESS): MCMC samples are autocorrelated, so 8000 draws might carry the information of only 800 independent ones. Low ESS means wide Monte Carlo error on your estimates; you want it in at least the hundreds, ideally thousands, for the quantities you report.
  • Divergences: HMC reports when its numerical integrator breaks down, usually in regions of sharp posterior curvature (the classic symptom of an unidentified or badly parameterized hierarchical model). Even a handful of divergences can mean the sampler is systematically missing part of the posterior. Do not ignore them.
  • Trace plots: the eyeball test. Healthy chains look like flat, overlapping “fuzzy caterpillars.” Trends, sticking, or chains exploring different regions are all alarms.

The honest cost of MCMC: it is slow (seconds to hours), it can fail silently if you skip the diagnostics, and getting a complex hierarchical model to sample cleanly can require real expertise in reparameterization (the non-centered parameterization trick, for instance, exists specifically to kill divergences in hierarchical models). This is the tax you pay for abandoning conjugacy. For a simple A/B test it is pure overkill — use the three-line conjugate version. For a model that pools information across forty experiments with partial pooling, there is no alternative, and the tax is worth it.

Where Bayesian Actually Beats Classical

Stripped of ideology, here is where the Bayesian approach earns its keep in engineering work:

  • Small samples and rare events. When you have ten data points or are estimating a failure rate from two failures in fifty thousand trials, the data alone is too thin to support a stable estimate. A defensible prior regularizes the answer and, just as importantly, gives you honest uncertainty instead of a point estimate that pretends to a precision it does not have.
  • Hierarchical / partial-pooling models. This is the killer app. Estimating conversion for 200 landing pages, some with millions of visitors and some with a dozen? A hierarchical model shares strength across pages: data-rich pages get their own estimate, data-poor pages get pulled toward the population mean exactly as far as the data warrants. This “partial pooling” is the principled middle ground between treating every page separately (overfitting the sparse ones) and lumping them all together (ignoring real differences). It is awkward to do well in the frequentist world and natural in the Bayesian one.
  • Online and sequential decisions. Optional stopping without penalty, plus multi-armed bandit algorithms like Thompson sampling, which is just “act in proportion to your posterior probability of being best.” Bayesian posteriors are built for streaming, sequential decision-making.
  • Intervals that mean what people think. A 95% credible interval has a 95% probability of containing the parameter, full stop. A 95% confidence interval does not — it means 95% of such intervals constructed over hypothetical repeated experiments would contain the true value, which is not a statement about the interval you are holding. Engineers and stakeholders consistently misread confidence intervals as credible intervals; the Bayesian version is the one they already believe they have.
  • Propagating uncertainty into decisions. Because the output is a full distribution, you can push it straight through a cost function and compute an expected loss or expected value, turning a statistical result directly into a decision. That pipeline — posterior to expected utility to action — is the entire point in an engineering setting.

And the flip side, stated just as plainly, because pretending otherwise is what gives Bayesians a bad name. Frequentist methods are fine or better when: you have a large sample (the prior washes out, the methods agree, and a closed-form t-test is faster and standard); you genuinely have no defensible prior and do not want to defend an arbitrary one; you are in a regulatory or adversarial context (clinical trials, safety certification) where pre-registered, prior-free procedures with guaranteed long-run error rates are the established, auditable standard; or you simply want a fast, well-understood, off-the-shelf answer and the philosophical nicety does not change the decision. There is nothing embarrassing about reaching for a t-test.

Verdict

Learn Bayes, but learn it as a tool and not a tribe. The framework’s real gift to engineers is that its output — a probability distribution over the unknown quantity — is the thing you actually need to make a decision, and it produces that without the interpretive traps that make p-values and confidence intervals so easy to misuse. For small-data estimation, hierarchical pooling, online sequential testing, and any problem where uncertainty must flow into a downstream choice, it is frequently the better instrument, and the conjugate cases cost you three lines of code. But the prior is a real responsibility, MCMC is a real expense with diagnostics you cannot skip, and on large samples or in regulated, prior-free contexts the frequentist toolkit is faster, standard, and entirely correct. The mature position is bilingual: know that the two paradigms answer different questions, pick the one that matches the question in front of you, and stop treating “frequentist” as an insult. The competent engineer owns both.

Sources

Comments