Math 04: The Bayesian Model, Likelihood & Inference

The previous notebooks built the parts — carryover, saturation, the baseline. This page wires them into a single generative probability model, writes down its likelihood and priors exactly as model/base.py::_build_model does, then runs the full inference pipeline on the Aurora data: prior predictive → NUTS → convergence → posterior intervals → posterior-predictive check → component decomposition. Every quoted number is measured from the baked, seeded notebook run (random_seed=0, 2 chains × 400 tune + 400 draws).

The joint model: priors × likelihood

The framework first standardizes the KPI, \(\tilde y_t = (y_t - \bar y)/s_y\) — for Aurora, \(\bar y = 839.0\) and \(s_y = 77.0\) (thousand dollars/week, \(T = 104\) weeks). On that scale it builds the additive linear predictor from math 00:

with the exponential saturation \(\operatorname{sat}(u) = 1 - e^{-\lambda u}\) (not Hill — see math 02) and the default adstock (since June 2026) a parametric in-graph geometric kernel, \(\operatorname{adstock}_c(x)_t = \sum_{k} w_k\, x_{t-k}\) with \(w_k \propto \alpha_c^{\,k}\) and a learned per-channel decay rate \(\alpha_c \sim \mathrm{Beta}(1,3)\) — a real decay rate, unlike the legacy path's \(\mathrm{Beta}(2,2)\) blend weight \(m_c\), which is demoted to an opt-out via use_parametric_adstock=False (see math 01). The likelihood is Gaussian on the standardized KPI:

A Bayesian model is this likelihood times a prior \(p(\theta)\) over every unknown — here \(\theta\) collects roughly twenty parameters (\(\beta_c, \lambda_c, \alpha_c, \gamma_j\), seasonal amplitudes, intercept, slope, \(\sigma\)). The object of inference is the posterior, by Bayes' theorem for the MMM:

The evidence \(p(\tilde y) = \int p(\tilde y\mid\theta)\,p(\theta)\,d\theta\) is an intractable integral over all \(\sim\!20\) dimensions, so we never compute it. Instead we sample from the unnormalized posterior — and a histogram of those draws is our posterior estimate. Every credible interval, ROAS band, and decomposition slice below comes from those draws.

The priors — what _build_model actually writes

The framework's priors are weakly informative: loose enough to let the data speak, tight enough to keep the sampler out of absurd regions and to encode mild domain knowledge (media effects are non-negative; carryover is bounded).

Two choices deserve a second look. The \(\mathrm{Gamma}(\mu{=}1.5, \sigma{=}1.0)\) on \(\beta_c\) is positive-only — it bakes in "more media does not reduce sales" (an experiment-calibrated prior would replace it: that is math 05). The \(\mathrm{Beta}(1,3)\) on the decay rate is right-skewed with mean \(0.25\): a-priori the model leans toward short memory and must be argued into long carryover by the data — a deliberate, conservative commitment. The notebook verifies the documented moments by assertion: the Gamma has mean exactly 1.5, sd exactly 1.0, and zero mass below 0; the Exponential has mean \(1/0.5 = 2\); the Beta(1,3) has mean 0.25 with most of its mass below 0.5.

Computed live in your browser — the same math as the notebook. The three media priors as density curves. Note the Gamma's hard wall at zero (positivity), the Exponential's mean at 2 (dashed), and the Beta(1,3)'s right-skewed mass with mean 0.25 (short memory unless the data argues otherwise): each encodes one modeling commitment.

Prior predictive: simulate sales before seeing data

Before fitting, ask: do these priors, run forward through the generative model, produce KPI series that look like plausible coffee sales? That is the prior predictive distribution

estimated by drawing \(\theta^{(s)} \sim p(\theta)\) and then \(\tilde y^{(s)} \sim p(\cdot \mid \theta^{(s)})\) — no data involved. A sane prior predictive brackets the observed range without being absurdly wide. One subtlety to expect: because \(\beta_c\) is positive-only and saturated media only ever adds, the envelope leans high — the prior's honest default belief is "media probably helps a bit." And under the parametric-adstock default (June 2026) the prior is noticeably more concentrated than the legacy blend's was, so expect bulk-not-blanket coverage of the observed weeks.

Illustrative recreation (seeded) — prior paths simulated in your browser from the same prior family (Gamma β, Exponential λ, Normal seasonality, HalfNormal σ) pushed through a simplified Aurora-shaped generative model; hit the button for a fresh batch. The measured notebook run (mmm.sample_prior_predictive(samples=300)): prior 1–99% envelope [703, 1602] vs observed range [690, 1035] (mean 839); 63% of observed weeks inside the band; prior-predictive median 1100 — same order as the data, leaning high, exactly as the positivity prior predicts.

NUTS in plain math: why gradients beat random walks

With \(\sim\!20\) correlated parameters, a random-walk Metropolis sampler is hopeless: to avoid rejection it must take tiny steps, and tiny random steps explore a 20-dimensional posterior at a glacial, diffusive crawl. (For the gentle build-up from grid Bayes to Metropolis to HMC, see workshop 02.)

Hamiltonian Monte Carlo fixes this by treating the negative log-posterior as a potential energy surface, \(U(\theta) = -\log\big[p(\tilde y \mid \theta)\, p(\theta)\big]\), and the current parameter vector as a particle. Each iteration draws a random momentum \(p\), then rolls the particle along the surface with leapfrog integration — discrete steps alternating gradient kicks and position drifts:

Because the proposal follows the gradient of the log-posterior, it glides along high-probability ridges instead of dithering — distant, high-acceptance proposals even in high dimensions. NUTS (the No-U-Turn Sampler) removes the one hard tuning knob, the trajectory length: it keeps integrating until the path starts to curve back on itself — a "U-turn" — then stops. The price of admission is that the whole model must be differentiable: the in-graph PyTensor adstock and saturation transforms of math 01–02 exist precisely to supply \(\nabla U\).

The notebook then fits for real: mmm.fit(draws=400, tune=400, chains=2, cores=1, random_seed=0) (cores=1 is required on macOS). Sampling 800 + 800 draws took 20 seconds.

Convergence: did the sampler work?

MCMC only yields the right posterior if the chains mixed. Two standard checks. \(\hat R\) (Gelman–Rubin) compares between-chain variance \(B\) to within-chain variance \(W\):

If the chains converged to the same distribution, \(B\) and \(W\) agree and \(\hat R \to 1\). Rule of thumb: \(\hat R < 1.01\) is great, \(> 1.1\) is a red flag. ESS corrects for autocorrelation: \(\text{ESS} \approx N \big/ \big(1 + 2\sum_{k\ge 1}\rho_k\big)\) — \(N\) correlated draws are worth fewer than \(N\) independent ones; a few hundred effective draws usually stabilizes means and intervals.

Measured — \(\hat R\) for the 14 key parameters from the baked mmm.summary() table. Every bar hugs 1: all 14 are below the 1.1 red line (max 1.010 on this subset), and ESS ranges 196–797. The sampler worked. Hold that thought — §Posterior shows what clean sampling does not guarantee.

What does \(\hat R\) actually see? Try it on two toy chains. Well-mixed chains wander over the same region, so between-chain variance matches within-chain variance and \(\hat R \approx 1\). Stuck chains each look healthy alone — the classic trap — but disagree about where the posterior is, and \(\hat R\) blows up.

Computed live in your browser — the same math as the notebook's \(\hat R\) definition. Two simulated 400-draw chains; the annotation shows \(\hat R\) computed from the formula above on exactly these draws. The trace-plot shape to memorize: a "fuzzy caterpillar" with both chains overlapping = good; two flat ribbons at different levels = chains that never found each other (each looks fine alone — only comparing chains exposes it).

Posterior intervals: what did we learn about each channel?

The posterior is a distribution, not a point. The natural summary is a forest plot: posterior mean plus a credible interval (here ArviZ's 94% HDI). The \(\beta_c\) intervals show both how big each media effect is and how sure the model is; the \(\lambda_c\) (saturation-rate) intervals are wide — that is the model honestly reporting that observational data weakly identifies saturation speed (the \(\beta\)–\(\lambda\) equifinality of math 02 and math 05); and the adstock decay rates \(\alpha_c\) come out wider than their short-memory \(\mathrm{Beta}(1,3)\) prior — prior–data tension, not learning (quantified in §Did the data teach us?).

Measured — posterior mean and 94% HDI per channel from the baked summary table. \(\beta\): Search 1.62 [0.04, 3.47] is the largest, TV 1.06 [0.03, 2.51]. \(\lambda\): every HDI is wide, spanning much of the prior — nothing is pinned. \(\alpha\): means 0.38–0.55, all above the prior's 0.25, with wide bands — the data pushing toward longer carryover than the Beta(1,3) prior wants.

Decision-makers want ROAS, not \(\beta\). The framework's MMMAnalyzer.compute_channel_roi() propagates the full posterior through a zero-out counterfactual, giving every channel a contribution HDI — divided by spend, a ROAS band. And because Aurora is synthetic, the true ROAS of every channel is known, so we can grade the model:

Measured — model ROAS (94% band from the contribution HDI) vs Aurora's known truth. The model is wrong everywhere: TV reads 0.33 vs a true 2.14, Display 0.67 vs 2.11 (both act through an unmodeled awareness mediator the single-equation model cannot see), while Search reads 1.05 vs a true 0.66 (its spend chases latent demand). The sampler converged perfectly on the wrong answer — clean sampling is not correct attribution. Fixing this takes mediation structure (extensions) or experiment calibration (math 05), not more draws.

The mathematics of that fix — how a lift test enters the model and pulls the posterior off the equifinality ridge — is worked through in math 05 · calibration, and its effect is measured across the full adversarial-worlds sequence in stress 05 · the gauntlet.

Posterior predictive check: does the fitted model reproduce the data?

The posterior predictive replays the generative model forward using the fitted posterior:

A well-fit model tracks the observed series and keeps most weeks inside its band — the band's width is the model's honest forecast uncertainty. mmm.predict(return_original_scale=True, hdi_prob=0.9) returns exactly this on the sales scale.

Illustrative recreation (seeded) of the posterior-predictive figure — an Aurora-shaped series with a 90% band of the same character. The measured notebook values: corr(predicted, observed) = 0.962, in-sample R² = 0.925, and 92% of weeks fall inside the 90% band — essentially nominal coverage. The fit is genuinely good. That is precisely what makes the ROAS bias above silent: every check on this page is green. See workshop 04 for how to read these checks as an analyst.

Decomposition: where do the sales come from?

The additive structure of \(\mu_t\) means fitted sales split exactly into interpretable pieces. Un-standardizing (\(\hat\mu_t = \tilde\mu_t\, s_y + \bar y\), with \(\bar y\) folded into the intercept) gives, in original sales units,

and mmm.compute_component_decomposition() returns each piece (dec.intercept, dec.trend, dec.seasonality, dec.controls_total, dec.media_by_channel) already on the sales scale. Because the model emits a Deterministic for every term in \(\mu_t\), the decomposition is additive and complete: baseline + \(\sum_c\) media reconstructs the fitted line.

Illustrative recreation (seeded), scaled to the measured component sizes — stacked baseline + per-channel media vs the predict() mean (dashed). Measured levels ($000s/week): baseline mean 747.4 (std 61.4), media means TV 19.4, Search 33.1, Social 21.1, Display 18.4 — the organic baseline (747) dwarfs total media (92). Seasonality and controls are mean-zero by construction but carry real amplitude (std 27.0 and 24.8): read them as swing around the level, not contributions to the average. Additivity, measured: corr(baseline + ∑ media, predict()) = 0.9999.

Did the data actually teach us each parameter?

One last honesty check the forest plots cannot give. A tight posterior interval tells you what you believe after the data — not how much of that belief came from the data versus the prior. A confident-looking posterior can be a faithful echo of an informative prior. The framework's compute_parameter_learning() compares prior draws to the posterior with three sample-based diagnostics:

\(c \to 1\): the data pinned it. \(c \approx 0\): prior-dominated. \(c < 0\): the posterior is wider than the prior — a flag, but an ambiguous one on its own. Read it with the shift: if \(|z| \geq 1\) the posterior moved at least one prior-sd without narrowing — the evidence dominated the location (verdict relocated), and the extra width reflects the likelihood's flatness in the newly-favored region, or heavy tails left by chains that mixed into it slowly. The diagnostic separates those two with an IQR-based contraction_robust (tail-insensitive) and post_ess_bulk (low ESS = suspect the sampler, not the posterior).

Measured — contraction per parameter from the baked compute_parameter_learning(prior_samples=2000) table, sorted so the least-learned parameters land on top. The headline under the parametric-adstock default is prior–data tension, not prior-domination: every adstock_alpha decay rate has negative contraction (c = −0.85 to −0.66, overlap 0.49–0.77) — the short-memory \(\mathrm{Beta}(1,3)\) prior commits to mean 0.25, the data argue for longer carryover, and the posterior comes out wider than the prior. For Search the likelihood outright wins the argument over location: shift z = +1.50 — the posterior moved 1.5 prior-sd without narrowing — verdict relocated (the evidence dominated where the parameter sits; the width is the likelihood's flatness, not "the data taught nothing"). The other three alphas stay verdict weak: the likelihood fights the prior without winning — observational data barely identifies carryover, the equifinality of math 01 as a measurement. Meanwhile \(\sigma\), the intercept, the controls, the saturation rates, and the Fourier amplitudes are clearly learned (c = 0.67–0.995); and beta_Search has c = −0.177 with shift z = +0.13 — its posterior, too, is wider than its prior. That is the demand-confounded channel (stress 03) flagging prior-data tension that a forest plot alone would hide.