.. _ch5_2-prior-distributions: ============================================== Section 5.2 Prior Specification and Conjugate Analysis ============================================== Every Bayesian analysis performed in :ref:`ch5_1-bayesian-foundations` relied quietly on a prior distribution that was already in hand. Grid approximation over :math:`\theta` requires knowing :math:`p(\theta)` before writing a single line of code. The analytical conjugate update for the vaccine trial opened with a Beta(2, 2) prior, chosen with a wave of the hand and a brief remark about "mild prior belief near 0.5." Yet the prior is not decoration β€” it is one of the two inputs that, together, determine the posterior. This section confronts directly the question that Section 5.1 deferred: *where does the prior come from, and what posterior does it produce?* The answer is not a single rule but a toolkit. We begin by examining the four distinct roles a prior plays in Bayesian analysis, then develop the five major *conjugate families* β€” prior-likelihood pairs for which the posterior has a tractable closed form. For each family we carry out the full arc: motivation, algebraic derivation, Python implementation, and worked numerical example. Having seen what conjugate priors *do*, we turn to the broader landscape of prior specification: Jeffreys priors motivated by invariance, weakly informative priors motivated by regularization, and informative priors grounded in genuine domain knowledge. The section closes with a systematic comparison of Bayesian and frequentist answers under each conjugate family and a demonstration of the limits of conjugate analysis β€” limits that motivate the Markov chain Monte Carlo machinery of Sections 5.4 and 5.5. The material here makes :ref:`ch3_1-exponential-families` pay off a second time: the exponential family structure that simplifies maximum likelihood is precisely the structure that makes conjugate analysis algebraically possible. Students who understand the sufficient statistic framework of Chapter 3 will recognize why conjugate pairs exist and how to construct them for new likelihoods. .. admonition:: Road Map 🧭 :class: important β€’ **Understand**: The four roles of the prior β€” information encoding, regularization, computation, and philosophical commitment β€” and how prior influence washes out asymptotically via the Bernstein–von Mises theorem β€’ **Develop**: Complete derivations of the five major conjugate families (Beta-Binomial, Normal-Normal, Normal-Inverse-Gamma, Poisson-Gamma, Multinomial-Dirichlet), with hyperparameter-as-pseudo-data interpretation and marginal likelihood computation for each β€’ **Implement**: Python classes for conjugate posterior computation; PyMC 5 implementations alongside each conjugate model; ArviZ posterior summaries and diagnostics; prior predictive simulation for model evaluation; prior sensitivity diagnostics β€’ **Evaluate**: Prior predictive checks, sensitivity analysis, and the algebraic limits of conjugate analysis that motivate subsequent MCMC sections .. admonition:: Software: PyMC 5 and ArviZ :class: tip Each conjugate model in this section is developed in two parallel implementations. The **analytical/scipy version** makes the mathematics explicit β€” every hyperparameter update is visible in the code. The **PyMC 5 version** (PyMC β‰₯ 5.28) shows the same model using the probabilistic programming interface that scales to non-conjugate posteriors in Sections 5.4–5.6. Both implementations are presented side by side so you can verify they agree to numerical precision and begin building intuition for the PyMC workflow before MCMC is formally introduced. Install with :code:`pip install pymc arviz` or :code:`conda install -c conda-forge pymc arviz`. The standard imports used throughout are: .. code-block:: python import pymc as pm # version 5.28.1 import arviz as az # ArviZ is installed automatically as a PyMC dependency The Role of the Prior --------------------- Before choosing a prior, it is worth being precise about why we need one and what it is doing. The prior is not simply a nuisance to be minimized or a placeholder to be replaced by flat density β€” it is a genuine component of the model, serving at least four separable functions. Four Functions of a Prior Distribution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Information encoding.** The most literal function: the prior encodes what is known about the parameter before the data arrive. A pharmaceutical researcher who has run three earlier trials of a related compound knows something about the plausible range of efficacy. A physicist measuring a fundamental constant brings tight prior knowledge from existing measurements. In both cases, ignoring the prior information and starting from a flat prior would be a deliberate act of forgetting β€” epistemically wasteful and sometimes practically harmful, as it yields wider posteriors than the evidence warrants. **Regularization.** When data are sparse relative to model complexity, the likelihood surface can be nearly flat or have multiple modes, making the MLE unreliable. A prior imposes a gentle preference for "simpler" or "smaller" parameter values that stabilizes estimation. Section 5.1 established the MAP-MLE bridge: under a Normal(0, :math:`\tau^2`) prior on regression coefficients, MAP estimation is algebraically equivalent to ridge regression, and under a Laplace prior it yields LASSO. This connection, developed in :ref:`ch3_5-generalized-linear-models`, shows that the most widely used regularization methods in machine learning are Bayesian priors in disguise. The advantage of the Bayesian framing is that the regularization strength :math:`1/\tau^2` acquires a principled interpretation as the prior precision, and can itself be estimated from data using hierarchical models (see :ref:`ch5_8-hierarchical-models`). **Computation.** Proper posteriors require proper priors (or likelihoods with sufficient data to rescue an improper prior). A prior that concentrates mass on a compact support ensures the posterior integral converges and posterior samples can be drawn. More specifically, conjugate priors β€” the main subject of this section β€” are chosen partly because they make the posterior integral analytically tractable. This is a computational motivation, not a philosophical one: conjugate priors were not invented because they are the "right" priors, but because they yield closed-form posteriors that make inference feasible without computers. **Regularizing pathological likelihoods.** Some models fail to have a proper MLE in finite samples. Logistic regression with perfect separation (a predictor that perfectly discriminates between outcomes) has an MLE of :math:`\pm \infty` for the relevant coefficient. A Normal prior on the coefficient immediately regularizes the problem, yielding a finite MAP estimate and a proper posterior. Bayesian inference with a prior is a more robust alternative to the ad hoc separation-handling heuristics common in software defaults. When the Prior Matters ~~~~~~~~~~~~~~~~~~~~~~ The influence of the prior on the posterior depends on the balance between prior precision and data information. The Beta-Binomial weighted average formula, derived as Exercise 3 in Section 5.1, gives the posterior mean as: .. math:: :label: weighted-average-general \mathbb{E}[\theta \mid y] = \underbrace{\frac{n_0}{n + n_0}}_{\text{weight on prior}} \mu_0 \;+\; \underbrace{\frac{n}{n + n_0}}_{\text{weight on MLE}} \hat{\theta}_{\text{MLE}} where :math:`n_0 = \alpha_0 + \beta_0` is the prior's *effective sample size* (the number of pseudo-observations it contributes) and :math:`\mu_0 = \alpha_0 / n_0` is the prior mean. This formula captures the prior-data balance precisely: when :math:`n \gg n_0`, the posterior mean converges to the MLE and the prior becomes irrelevant; when :math:`n \ll n_0`, the prior dominates. Every conjugate family has an analogous formula β€” the differences lie only in how the "effective sample size" is defined and what quantity plays the role of the MLE. We will make this pattern explicit for each family in the sections that follow. The Bernstein–von Mises theorem formalizes the asymptotic irrelevance of the prior. As :math:`n \to \infty`, under mild regularity conditions, any two posteriors that start from different but proper priors converge to the same Normal distribution centered at the MLE with variance equal to the inverse Fisher information (see :ref:`ch3_3-sampling-variability`). In large samples, the data swamp the prior. The practical implication: when :math:`n` is large, prior specification matters less than model specification (the choice of likelihood). When :math:`n` is small β€” the early stages of a clinical trial, a rare event study, or a high-dimensional model with few observations per parameter β€” the prior is consequential and deserves careful thought. .. admonition:: Common Pitfall ⚠️ :class: warning Choosing a "flat" or "noninformative" prior does not guarantee that your analysis is objective or that the prior has no influence. A flat prior on :math:`\theta` is not flat on :math:`g(\theta) = \log(\theta/(1-\theta))`. Priors do not transform as constants β€” they transform as densities, picking up a Jacobian factor from the change of variables. This means "flat" is parameterisation-dependent: the same state of ignorance about :math:`\theta` implies a non-flat prior on any nonlinear function of :math:`\theta`, so the choice of what to be flat *in* is itself an arbitrary modelling decision. A prior is **invariant under reparameterisation** if it produces the same posterior inference regardless of which parameterisation is used β€” that is, the prior on :math:`\phi = g(\theta)` is exactly what you would get by transforming the prior on :math:`\theta` through the change-of-variables formula. The Jeffreys prior (Block 7a below) is the principled solution to this problem: it is derived to satisfy this invariance property exactly. Until you have read that section, treat "uninformative" priors with caution and verify that your conclusions do not depend strongly on the parameterisation convention. The Prior Spectrum ~~~~~~~~~~~~~~~~~~ For expository clarity, we organize prior choices along a spectrum from fully data-driven to fully expert-driven: 1. **Conjugate priors** (Blocks 2–6): algebraically tractable; posterior in same family as prior; hyperparameters interpretable as pseudo-observations. The starting point for Bayesian analysis of standard models. 2. **Jeffreys priors** (Block 7a): derived from Fisher information; invariant under reparameterization; the principled "objective Bayesian" choice for single-parameter problems. 3. **Weakly informative priors** (Block 7b): designed to regularize gently without imposing strong beliefs; the modern default for regression-style models. 4. **Informative priors** (Block 7c): encode genuine domain knowledge from previous studies, physical constraints, or expert judgment. The sections below develop each family and type in this order. We start with conjugate priors because they are the foundation: understanding what a tractable posterior looks like, and what the hyperparameters mean, is the prerequisite for evaluating all other prior choices intelligently. .. _sec-beta-binomial: Beta-Binomial Conjugate Analysis --------------------------------- The Beta-Binomial model is the foundational example of conjugate analysis. It applies whenever data are binary (success/failure, present/absent, positive/negative) and we wish to estimate the success probability :math:`\theta`. We encountered this model in Section 5.1 via the vaccine trial example and grid approximation. Here we develop it analytically, deriving everything β€” the conjugate update rule, the posterior moments, the marginal likelihood, and the posterior predictive β€” from first principles. Mathematical Framework ~~~~~~~~~~~~~~~~~~~~~~ Let :math:`k` successes be observed in :math:`n` independent Bernoulli(:math:`\theta`) trials, so that the likelihood is: .. math:: p(k \mid \theta, n) = \binom{n}{k} \theta^k (1 - \theta)^{n-k}, \qquad \theta \in [0, 1] We place a Beta(:math:`\alpha_0, \beta_0`) prior on :math:`\theta`: .. math:: p(\theta) = \frac{\theta^{\alpha_0 - 1}(1-\theta)^{\beta_0 - 1}}{B(\alpha_0, \beta_0)}, \qquad \alpha_0, \beta_0 > 0 where :math:`B(\alpha_0, \beta_0) = \Gamma(\alpha_0)\Gamma(\beta_0)/\Gamma(\alpha_0 + \beta_0)` is the Beta function. The prior mean is :math:`\mu_0 = \alpha_0 / (\alpha_0 + \beta_0)` and the prior variance is :math:`\alpha_0 \beta_0 / \left[(\alpha_0 + \beta_0)^2(\alpha_0 + \beta_0 + 1)\right]`. By Bayes' theorem: .. math:: p(\theta \mid k, n) \propto p(k \mid \theta, n) \cdot p(\theta) **Full derivation.** Drop all multiplicative constants that do not depend on :math:`\theta`: .. math:: :label: beta-binomial-kernel p(\theta \mid k, n) &\propto \theta^k (1 - \theta)^{n-k} \cdot \theta^{\alpha_0 - 1}(1-\theta)^{\beta_0 - 1} \\ &= \theta^{(\alpha_0 + k) - 1}(1-\theta)^{(\beta_0 + n - k) - 1} This is exactly the kernel of a :math:`\text{Beta}(\alpha_0 + k,\; \beta_0 + n - k)` distribution. Because the posterior integrates to one, the normalizing constant must equal :math:`B(\alpha_0 + k,\; \beta_0 + n - k)^{-1}`, and we have the **conjugate update rule**: .. admonition:: Definition: Beta-Binomial Conjugate Update :class: note If :math:`\theta \sim \text{Beta}(\alpha_0, \beta_0)` and :math:`k \mid \theta, n \sim \text{Binomial}(n, \theta)`, then the posterior is: .. math:: :label: beta-binomial-posterior \theta \mid k, n \;\sim\; \text{Beta}(\alpha_0 + k,\; \beta_0 + n - k) The update rule is addition: **prior successes plus observed successes** in the first parameter; **prior failures plus observed failures** in the second. This remarkable simplicity β€” the posterior is the same Beta family, with hyperparameters updated by simple addition β€” is the defining property of a conjugate prior. The algebra of kernel recognition (collect powers of :math:`\theta` and :math:`(1-\theta)`, identify the resulting expression as a named density's kernel) is the same technique used in every conjugate derivation in this section. Posterior Moments and Summaries ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Denoting the posterior parameters as :math:`\alpha_n = \alpha_0 + k` and :math:`\beta_n = \beta_0 + n - k`, the Beta distribution's standard results give: .. math:: :label: beta-posterior-moments \mathbb{E}[\theta \mid k, n] &= \frac{\alpha_n}{\alpha_n + \beta_n} = \frac{\alpha_0 + k}{\alpha_0 + \beta_0 + n} \\[6pt] \text{Mode}(\theta \mid k, n) &= \frac{\alpha_n - 1}{\alpha_n + \beta_n - 2} = \frac{\alpha_0 + k - 1}{\alpha_0 + \beta_0 + n - 2}, \quad \alpha_n, \beta_n > 1 \\[6pt] \text{Var}(\theta \mid k, n) &= \frac{\alpha_n \beta_n}{(\alpha_n + \beta_n)^2(\alpha_n + \beta_n + 1)} Expanding the posterior mean with :math:`n_0 = \alpha_0 + \beta_0` and :math:`\hat{\theta} = k/n`: .. math:: :label: beta-weighted-average \mathbb{E}[\theta \mid k, n] = \frac{n_0}{n_0 + n} \cdot \frac{\alpha_0}{n_0} + \frac{n}{n_0 + n} \cdot \frac{k}{n} = \frac{n_0}{n_0 + n} \mu_0 + \frac{n}{n_0 + n} \hat{\theta} This is the weighted average formula from Section 5.1, now derived as a *consequence* of the conjugate update rather than presented as an observation. The weight on the data, :math:`w = n/(n + n_0)`, grows toward 1 as :math:`n \to \infty` regardless of the prior β€” a concrete instance of the Bernstein–von Mises principle. Hyperparameter Interpretation: The Pseudo-Data View ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Beta hyperparameters :math:`(\alpha_0, \beta_0)` have a natural interpretation: they represent *pseudo-observations* from an imagined earlier experiment in which :math:`\alpha_0` successes and :math:`\beta_0` failures were observed. The total prior "sample size" is :math:`n_0 = \alpha_0 + \beta_0`. This interpretation is not merely mnemonic β€” it has operational content. A prior analyst who claims Beta(10, 10) is asserting that their prior knowledge is equivalent to having already seen 20 independent Bernoulli trials, 10 of each type. Whether that claim is defensible depends on the application. This pseudo-data interpretation extends to all conjugate families β€” it will be the lens through which we interpret Normal, Gamma, and Dirichlet hyperparameters as well. Marginal Likelihood ^^^^^^^^^^^^^^^^^^^^ The marginal likelihood β€” the probability of the observed data averaged over the prior β€” is obtained by computing the normalizing constant of the unnormalized posterior: .. math:: :label: beta-marginal-likelihood p(k \mid n, \alpha_0, \beta_0) = \binom{n}{k} \frac{B(\alpha_0 + k,\; \beta_0 + n - k)}{B(\alpha_0, \beta_0)} This is the *Beta-Binomial distribution* for :math:`k`. The marginal likelihood plays a central role in Bayesian model comparison (discussed in :ref:`ch5_7-model-comparison`): the ratio of marginal likelihoods for two models with different priors (or different likelihoods) is the Bayes factor, which provides a principled measure of relative evidence. Posterior Predictive Distribution ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Given the posterior Beta(:math:`\alpha_n, \beta_n`), the predictive probability of success in a *new* independent trial is: .. math:: P(\tilde{y} = 1 \mid k, n) = \int_0^1 \theta \cdot p(\theta \mid k, n) \, d\theta = \mathbb{E}[\theta \mid k, n] = \frac{\alpha_n}{\alpha_n + \beta_n} The posterior mean is the predictive probability of success β€” a satisfying result that formalizes the Bayesian rule for probability assessment: your probability for the next outcome should equal your expected value of the chance parameter. For predicting :math:`m` new trials, the posterior predictive count :math:`\tilde{k} \mid k, n, m` follows a Beta-Binomial distribution. Python Implementation ~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python import numpy as np from scipy import stats from scipy.special import betaln import matplotlib.pyplot as plt class BetaBinomialModel: """ Conjugate Beta-Binomial model for Bernoulli proportion inference. Parameters ---------- alpha0 : float First hyperparameter of the Beta prior (> 0). Interpretable as alpha0 pseudo-successes. beta0 : float Second hyperparameter of the Beta prior (> 0). Interpretable as beta0 pseudo-failures. """ def __init__(self, alpha0: float, beta0: float) -> None: if alpha0 <= 0 or beta0 <= 0: raise ValueError("Beta hyperparameters must be positive.") self.alpha0 = alpha0 self.beta0 = beta0 # Posterior hyperparameters start at prior values self.alpha_n = alpha0 self.beta_n = beta0 self._n_obs = 0 # cumulative observations # ------------------------------------------------------------------ # # Posterior updating # # ------------------------------------------------------------------ # def update(self, k: int, n: int) -> "BetaBinomialModel": """ Return a *new* model with the posterior after observing k successes in n trials. The original model is not mutated (immutable style). Parameters ---------- k : int Number of observed successes (0 <= k <= n). n : int Total number of trials. Returns ------- BetaBinomialModel Updated model with posterior parameters. """ if not (0 <= k <= n): raise ValueError("Require 0 <= k <= n.") new = BetaBinomialModel(self.alpha0, self.beta0) new.alpha_n = self.alpha_n + k new.beta_n = self.beta_n + (n - k) new._n_obs = self._n_obs + n return new # ------------------------------------------------------------------ # # Posterior summaries # # ------------------------------------------------------------------ # @property def posterior_mean(self) -> float: """Posterior mean: alpha_n / (alpha_n + beta_n).""" return self.alpha_n / (self.alpha_n + self.beta_n) @property def posterior_mode(self) -> float | None: """Posterior mode; None if alpha_n <= 1 or beta_n <= 1.""" if self.alpha_n <= 1 or self.beta_n <= 1: return None return (self.alpha_n - 1) / (self.alpha_n + self.beta_n - 2) @property def posterior_variance(self) -> float: """Posterior variance from Beta distribution formula.""" a, b = self.alpha_n, self.beta_n return a * b / ((a + b) ** 2 * (a + b + 1)) @property def effective_prior_n(self) -> float: """Effective prior sample size: alpha0 + beta0.""" return self.alpha0 + self.beta0 @property def prior_data_weight(self) -> float: """Weight on data MLE: n / (n + n0). Zero before any update.""" n0 = self.effective_prior_n n = self._n_obs return n / (n + n0) if n > 0 else 0.0 def credible_interval(self, level: float = 0.95) -> tuple[float, float]: """ Equal-tailed credible interval using Beta quantile function. Parameters ---------- level : float Nominal coverage probability (default 0.95). Returns ------- tuple[float, float] (lower, upper) bounds of the credible interval. """ alpha = 1 - level dist = stats.beta(self.alpha_n, self.beta_n) return dist.ppf(alpha / 2), dist.ppf(1 - alpha / 2) # ------------------------------------------------------------------ # # Marginal likelihood # # ------------------------------------------------------------------ # def log_marginal_likelihood(self, k: int, n: int) -> float: """ Log marginal likelihood p(k | n, alpha0, beta0), the Beta-Binomial probability mass, computed using log-Beta for numerical stability. Parameters ---------- k : int Number of successes. n : int Number of trials. Returns ------- float Log marginal likelihood (natural log). """ from scipy.special import gammaln # log C(n, k) + log B(a+k, b+n-k) - log B(a, b) log_binom = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1) log_b_post = betaln(self.alpha0 + k, self.beta0 + n - k) log_b_prior = betaln(self.alpha0, self.beta0) return log_binom + log_b_post - log_b_prior # ------------------------------------------------------------------ # # Posterior predictive # # ------------------------------------------------------------------ # def posterior_predictive(self, m: int = 1, n_samples: int = 10_000, rng: np.random.Generator | None = None ) -> np.ndarray: """ Monte Carlo samples from the posterior predictive distribution for m new Bernoulli trials. Parameters ---------- m : int Number of new trials to predict. n_samples : int Monte Carlo sample size. rng : np.random.Generator, optional Random number generator for reproducibility. Returns ------- np.ndarray, shape (n_samples,) Predicted counts of successes in m new trials. """ if rng is None: rng = np.random.default_rng() # Step 1: Draw theta from posterior Beta theta_samples = rng.beta(self.alpha_n, self.beta_n, size=n_samples) # Step 2: Draw counts from Binomial(m, theta) return rng.binomial(m, theta_samples) .. admonition:: Example πŸ’‘ Drug Trial Responder Rate Analysis :class: note **Given:** A Phase II clinical trial tests a new drug on :math:`n = 20` patients. :math:`k = 15` patients show a positive response. Three clinicians propose different priors: - **Clinician A** (skeptical, no prior data): Beta(1, 1) β€” flat prior, all response rates equally likely - **Clinician B** (standard Jeffreys): Beta(0.5, 0.5) β€” invariant reference prior (derived in Block 7a) - **Clinician C** (informed by pilot study): Beta(10, 10) β€” prior equivalent to 10 successes and 10 failures from earlier work, suggesting response rate near 50% **Find:** The posterior distribution of the true response rate :math:`\theta` under each prior, and compare with the MLE and a 95% frequentist Wald confidence interval. **Mathematical approach:** The MLE is :math:`\hat{\theta} = 15/20 = 0.750`. The Wald interval is :math:`\hat{\theta} \pm 1.96\sqrt{\hat{\theta}(1-\hat{\theta})/n} = 0.750 \pm 0.190 = (0.560, 0.940)`. Applying the conjugate update rule :eq:`beta-binomial-posterior` to each prior: .. math:: \text{Clinician A: } \text{Beta}(1+15, 1+5) &= \text{Beta}(16, 6) \\ \text{Clinician B: } \text{Beta}(0.5+15, 0.5+5) &= \text{Beta}(15.5, 5.5) \\ \text{Clinician C: } \text{Beta}(10+15, 10+5) &= \text{Beta}(25, 15) Posterior means and effective prior weights (:math:`w = n/(n + n_0)`): .. math:: \text{A: } \mu_n = \frac{16}{22} \approx 0.727, \quad w = \frac{20}{22} = 0.909 \\ \text{B: } \mu_n = \frac{15.5}{21} \approx 0.738, \quad w = \frac{20}{21} = 0.952 \\ \text{C: } \mu_n = \frac{25}{40} = 0.625, \quad w = \frac{20}{40} = 0.500 Clinician C's informative prior (with effective sample size :math:`n_0 = 20`) pulls the posterior mean substantially toward 0.5, because the prior contributes as much information as the data. Clinician B's Jeffreys prior contributes only 1 pseudo-observation (:math:`n_0 = 1`), leaving the data to dominate. **Python implementation:** .. code-block:: python import numpy as np from scipy import stats # Data n, k = 20, 15 mle = k / n # 0.75 # Frequentist Wald CI se = np.sqrt(mle * (1 - mle) / n) wald_ci = (mle - 1.96 * se, mle + 1.96 * se) print(f"MLE: {mle:.4f}") print(f"Wald 95% CI: ({wald_ci[0]:.4f}, {wald_ci[1]:.4f})") print() # Three priors priors = { 'Flat Beta(1,1)': (1.0, 1.0), 'Jeffreys Beta(.5,.5)': (0.5, 0.5), 'Skeptical Beta(10,10)': (10.0, 10.0), } print(f"{'Prior':<25} {'Post Mean':>10} {'Post Mode':>10} " f"{'95% CrI':>22} {'w_data':>8}") print("-" * 80) for name, (a0, b0) in priors.items(): model = BetaBinomialModel(a0, b0).update(k, n) lo, hi = model.credible_interval(0.95) mode = model.posterior_mode print(f"{name:<25} {model.posterior_mean:>10.4f} " f"{mode if mode else 'N/A':>10} " f"({lo:.4f}, {hi:.4f}) {model.prior_data_weight:>8.3f}") # Log marginal likelihoods (for model comparison) print("\nLog Marginal Likelihoods:") for name, (a0, b0) in priors.items(): m0 = BetaBinomialModel(a0, b0) lml = m0.log_marginal_likelihood(k, n) print(f" {name}: {lml:.4f}") .. code-block:: text MLE: 0.7500 Wald 95% CI: (0.5602, 0.9398) Prior Post Mean Post Mode 95% CrI w_data -------------------------------------------------------------------------------- Flat Beta(1,1) 0.7273 0.7500 (0.5283, 0.8872) 0.909 Jeffreys Beta(.5,.5) 0.7381 0.7632 (0.5358, 0.8976) 0.952 Skeptical Beta(10,10) 0.6250 0.6316 (0.4718, 0.7664) 0.500 **Result:** The flat and Jeffreys priors produce posteriors close to the frequentist results β€” posterior means of 0.727 and 0.738 vs. the MLE of 0.750, credible intervals overlapping substantially with the Wald interval. The skeptical prior, contributing 20 pseudo-observations matched against the 20 real observations, pulls the posterior mean to 0.625 and substantially narrows the credible interval compared to the data-only analysis. All three posterior intervals are slightly narrower than the Wald interval because they incorporate the Beta distribution's natural truncation to :math:`[0, 1]`, avoiding the Wald interval's infamous boundary issue (it can produce limits outside :math:`[0,1]` for small :math:`n`). .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig09_beta_binomial_priors.png :alt: Three-panel figure showing prior-to-posterior updating for the drug trial example under flat, Jeffreys, and skeptical Beta priors, with posterior means and credible intervals annotated and the MLE marked in each panel :align: center :width: 100% **Figure 5.9**: Beta-Binomial conjugate updating for the drug trial (:math:`k=15, n=20`). Each panel overlays the prior (dashed) and posterior (solid) density. Under the flat prior (left), the posterior tracks the likelihood closely, centered near the MLE of 0.75. The Jeffreys prior (center) produces nearly the same result. The skeptical Beta(10,10) prior (right) has substantial mass near 0.5 and pulls the posterior mean to 0.625, illustrating the prior-data balance quantified by the weight formula :eq:`beta-weighted-average`. Vertical dashed lines mark the prior mean, MLE, and posterior mean in each panel. .. admonition:: PyMC 5 Implementation: Beta-Binomial Model :class: note The same drug trial analysis in PyMC requires minimal code: specify the prior, the likelihood with observed data, and call :code:`pm.sample`. PyMC's default sampler (NUTS) is far more powerful than needed for a one-parameter conjugate model β€” here it is used primarily to introduce the workflow that generalises to non-conjugate posteriors in later sections. .. code-block:: python import pymc as pm import arviz as az import numpy as np # Drug trial: k=15 successes in n=20 trials k_obs, n_trials = 15, 20 with pm.Model() as beta_binomial_model: # Prior: flat Beta(1, 1) β€” equivalent to Uniform(0, 1) theta = pm.Beta('theta', alpha=1.0, beta=1.0) # Likelihood: k successes in n trials k = pm.Binomial('k', n=n_trials, p=theta, observed=k_obs) # Sample from the posterior idata = pm.sample( draws=4000, tune=1000, chains=4, random_seed=42, progressbar=False, ) # ── Posterior summaries via ArviZ ────────────────────────────────── # az.summary returns a DataFrame with mean, sd, HDI, MCSE, ESS, R-hat summary = az.summary(idata, var_names=['theta'], hdi_prob=0.95) print(summary) # Extract flat array of posterior draws across all chains theta_samples = idata.posterior['theta'].values.flatten() # shape: (n_chains * n_draws,) print(f"\nMCMC posterior mean: {theta_samples.mean():.4f}") print(f"Analytical posterior mean: {(1 + k_obs) / (1 + k_obs + 1 + n_trials - k_obs):.4f}") print(f"95% HDI (MCMC): [{np.quantile(theta_samples, 0.025):.4f}, " f"{np.quantile(theta_samples, 0.975):.4f}]") # ── Visualize with ArviZ ──────────────────────────────────────────── # az.plot_posterior shows the KDE, HDI, and point estimate in one panel az.plot_posterior( idata, var_names=['theta'], hdi_prob=0.95, ref_val=0.5, # mark the "better than chance" threshold point_estimate='mean', ) **Expected output (az.summary):** .. code-block:: text mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk ess_tail r_hat theta 0.7273 0.0943 0.5310 0.8973 0.0007 0.0006 16407.0 11732.0 1.00 **Reading the output.** The MCMC posterior mean (0.7273) matches the analytical Beta(16, 6) mean to four decimal places β€” as expected for a conjugate model where NUTS converges in one step. The columns :code:`ess_bulk` and :code:`ess_tail` report effective sample size: with 4 chains Γ— 4,000 draws = 16,000 total samples, bulk ESS β‰ˆ 16,407 confirms near-perfect mixing. The :code:`r_hat` value of 1.00 (target < 1.01) confirms convergence. For a conjugate model these diagnostics are trivially good; they become critical in non-conjugate settings from Section 5.5 onward. :code:`az.plot_posterior` displays the posterior density with the 95% HDI annotated and a reference line at :math:`\theta = 0.5`. The percentage below and above the reference value is shown automatically β€” a one-click ROPE-style summary. .. admonition:: PyMC vs Analytical β€” When to Use Which :class: tip For Beta-Binomial (and the other four conjugate families), the analytical update is always preferable: it is exact, instant, and requires no convergence diagnostics. Use PyMC when the model is non-conjugate (logistic regression, mixture models, hierarchical models), when you want the full MCMC infrastructure (convergence checks, posterior predictive, LOO), or when you are prototyping a model that will grow in complexity. The conjugate code in this section is the reference against which PyMC output is validated. .. _sec-normal-normal: Normal-Normal Model: Known Variance ------------------------------------- The Beta-Binomial model concerns a bounded probability parameter. The Normal-Normal conjugate pair addresses the complementary setting: estimating a location parameter :math:`\mu` on the real line when the observation variance :math:`\sigma^2` is known. This model is foundational for understanding precision-weighted inference β€” the principle that more precise measurements contribute more to the posterior mean β€” and it will be generalized in the Normal-Inverse-Gamma model when :math:`\sigma^2` is also unknown. Mathematical Framework ~~~~~~~~~~~~~~~~~~~~~~ Let :math:`y_1, \ldots, y_n` be i.i.d. observations from :math:`\mathcal{N}(\mu, \sigma^2)` with :math:`\sigma^2` known. The likelihood for the mean :math:`\mu` is: .. math:: p(y_1, \ldots, y_n \mid \mu) \propto \exp\!\left(-\frac{n}{2\sigma^2}(\bar{y} - \mu)^2\right) where :math:`\bar{y} = n^{-1}\sum_{i=1}^n y_i` is the sufficient statistic. Place a Normal prior on :math:`\mu`: .. math:: \mu \sim \mathcal{N}(\mu_0, \sigma_0^2) The parameter :math:`\mu_0` is the prior mean and :math:`\sigma_0^2` is the prior variance. In Bayesian Normal analysis it is natural to work with *precisions* rather than variances: define .. math:: \tau_0 = \frac{1}{\sigma_0^2} \quad \text{(prior precision)}, \qquad \tau_{\text{data}} = \frac{n}{\sigma^2} \quad \text{(data precision for } \bar{y}\text{)} **Full derivation.** The unnormalized log-posterior is the log-likelihood plus log-prior: .. math:: \log p(\mu \mid y) &\propto -\frac{n}{2\sigma^2}(\bar{y} - \mu)^2 - \frac{1}{2\sigma_0^2}(\mu - \mu_0)^2 \\ &= -\frac{1}{2}\left[\tau_{\text{data}}(\bar{y} - \mu)^2 + \tau_0(\mu - \mu_0)^2\right] Expanding the quadratic forms and collecting terms in :math:`\mu`: .. math:: \log p(\mu \mid y) \propto -\frac{1}{2}(\tau_0 + \tau_{\text{data}})\left(\mu - \frac{\tau_0 \mu_0 + \tau_{\text{data}} \bar{y}}{\tau_0 + \tau_{\text{data}}}\right)^2 This is the log-kernel of a Normal distribution with *posterior precision* :math:`\tau_n = \tau_0 + \tau_{\text{data}}` and *posterior mean* :math:`\mu_n`: .. admonition:: Definition: Normal-Normal Conjugate Update (Known Variance) :class: note If :math:`\mu \sim \mathcal{N}(\mu_0, \tau_0^{-1})` and :math:`y_1, \ldots, y_n \mid \mu \stackrel{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^2)` with :math:`\sigma^2` known, then: .. math:: :label: normal-normal-posterior \mu \mid y \;\sim\; \mathcal{N}\!\left(\mu_n,\; \tau_n^{-1}\right) where the posterior precision and posterior mean satisfy: .. math:: :label: normal-normal-update \tau_n &= \tau_0 + \frac{n}{\sigma^2} \qquad \text{(precisions add)} \\[6pt] \mu_n &= \frac{\tau_0 \mu_0 + \tau_{\text{data}} \bar{y}}{\tau_n} = \frac{\tau_0}{\tau_n}\mu_0 + \frac{\tau_{\text{data}}}{\tau_n}\bar{y} The key structural result is that **precisions add**: the posterior precision equals the prior precision plus the data precision, so each additional observation increases the precision of inference by :math:`1/\sigma^2`. The posterior mean is a precision-weighted average of the prior mean and the sample mean, exactly analogous to the Beta-Binomial weighted average formula :eq:`beta-weighted-average`. Effective prior sample size and weight ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The data precision for a single observation is :math:`1/\sigma^2`, so :math:`\tau_{\text{data}} = n/\sigma^2` corresponds to :math:`n` observations. The prior precision :math:`\tau_0 = 1/\sigma_0^2` is equivalent to :math:`n_0 = \sigma^2/\sigma_0^2` pseudo-observations. The weight on the data is therefore: .. math:: :label: normal-normal-weight w = \frac{\tau_{\text{data}}}{\tau_n} = \frac{n/\sigma^2}{1/\sigma_0^2 + n/\sigma^2} = \frac{n}{n + \sigma^2/\sigma_0^2} As :math:`n \to \infty`, :math:`w \to 1` and the posterior mean converges to the sample mean :math:`\bar{y}`. As :math:`\sigma_0^2 \to \infty` (prior becomes diffuse), :math:`w \to 1` immediately regardless of :math:`n` β€” a diffuse prior yields a posterior mean equal to the MLE. Sequential Updating ^^^^^^^^^^^^^^^^^^^^ The Normal-Normal update rule is especially elegant in the sequential setting. Suppose observations arrive one at a time. After the :math:`t`-th observation :math:`y_t`, the current posterior becomes the prior for the next update. The precision after :math:`t` observations is: .. math:: \tau^{(t)} = \tau_0 + \frac{t}{\sigma^2} which increases linearly in :math:`t`. The posterior mean after :math:`t` observations satisfies the recursion: .. math:: :label: normal-normal-recursion \mu^{(t)} = \frac{\tau^{(t-1)}\mu^{(t-1)} + y_t / \sigma^2}{\tau^{(t-1)} + 1/\sigma^2} This is a one-pass online algorithm: process each observation, update the sufficient statistics :math:`(\tau, \mu)`, and discard the raw observation. The final posterior is identical to processing all :math:`n` observations in a single batch β€” sequential Bayesian updating is *order-invariant* and *memoryless* in the sense that the full history is summarized without loss by the pair :math:`(\tau_n, \mu_n)`. The same principle governs the Kalman filter, which applies the Normal-Normal recursion to time series with evolving state vectors. Posterior Predictive Distribution ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For a new observation :math:`\tilde{y}` from the same :math:`\mathcal{N}(\mu, \sigma^2)` model, the posterior predictive is: .. math:: \tilde{y} \mid y_1, \ldots, y_n \;\sim\; \mathcal{N}\!\left(\mu_n,\; \sigma^2 + \tau_n^{-1}\right) The predictive variance has two additive components: :math:`\sigma^2` from the sampling variability of new data given :math:`\mu`, and :math:`\tau_n^{-1} = \sigma_n^2` from the posterior uncertainty in :math:`\mu`. The Bayesian predictive interval is therefore wider than the plug-in prediction interval :math:`\mathcal{N}(\bar{y}, \sigma^2)`, which treats :math:`\hat{\mu} = \bar{y}` as known exactly. Python Implementation ~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python import numpy as np from scipy import stats from dataclasses import dataclass @dataclass class NormalNormalModel: """ Conjugate Normal-Normal model for mean inference with known variance. Stores sufficient statistics (posterior precision and posterior mean) for sequential or batch updating. Parameters ---------- mu0 : float Prior mean. sigma2_prior : float Prior variance (> 0). sigma2_known : float Known observation variance (> 0). """ mu0: float sigma2_prior: float sigma2_known: float # Posterior sufficient statistics (initialized to prior) def __post_init__(self) -> None: self._tau_n = 1.0 / self.sigma2_prior # posterior precision self._mu_n = self.mu0 # posterior mean self._n = 0 # observations processed # ------------------------------------------------------------------ # # Updating # # ------------------------------------------------------------------ # def update(self, data: np.ndarray) -> "NormalNormalModel": """ Return updated model after observing data (batch update). Parameters ---------- data : np.ndarray Array of observed values. Returns ------- NormalNormalModel New model instance with posterior parameters. """ new = NormalNormalModel(self.mu0, self.sigma2_prior, self.sigma2_known) new._n = self._n + len(data) tau_data = len(data) / self.sigma2_known new._tau_n = self._tau_n + tau_data new._mu_n = (self._tau_n * self._mu_n + tau_data * data.mean()) / new._tau_n return new def update_single(self, x: float) -> "NormalNormalModel": """ Sequential update: process one observation at a time. Equivalent to batch update for a single observation. """ return self.update(np.array([x])) # ------------------------------------------------------------------ # # Posterior summaries # # ------------------------------------------------------------------ # @property def posterior_mean(self) -> float: return self._mu_n @property def posterior_variance(self) -> float: return 1.0 / self._tau_n @property def posterior_std(self) -> float: return np.sqrt(self.posterior_variance) def credible_interval(self, level: float = 0.95) -> tuple[float, float]: """Equal-tailed Normal credible interval.""" z = stats.norm.ppf(0.5 + level / 2) return (self._mu_n - z * self.posterior_std, self._mu_n + z * self.posterior_std) # ------------------------------------------------------------------ # # Posterior predictive # # ------------------------------------------------------------------ # def posterior_predictive_distribution(self) -> stats.rv_continuous: """ Returns a frozen scipy Normal distribution for the posterior predictive of a new observation: N(mu_n, sigma2 + sigma2_n). """ pred_var = self.sigma2_known + self.posterior_variance return stats.norm(loc=self._mu_n, scale=np.sqrt(pred_var)) # ------------------------------------------------------------------ # # Sequence of posteriors (for visualization) # # ------------------------------------------------------------------ # def sequential_posteriors(self, data: np.ndarray ) -> list["NormalNormalModel"]: """ Process data one observation at a time and return a list of model objects (including the prior, as the first element). """ models = [self] current = self for x in data: current = current.update_single(x) models.append(current) return models .. admonition:: Example πŸ’‘ Temperature Sensor Calibration :class: note **Given:** A temperature sensor is certified by the manufacturer to measure with standard deviation :math:`\sigma = 0.5 \,^\circ\text{C}` (so :math:`\sigma^2 = 0.25`). Before deployment, ten calibration readings are taken at a reference temperature. From the instrument's data sheet, a technician assigns a prior :math:`\mu \sim \mathcal{N}(20.0,\; 1.0^2)` β€” the reference temperature is nominally :math:`20.0^\circ\text{C}` with considerable uncertainty about the sensor's true offset. **Find:** The posterior distribution of the sensor's true reading bias :math:`\mu` after processing 10 calibration measurements, using both batch and sequential updating. Verify that both routes yield identical results. **Mathematical approach:** Prior precision: :math:`\tau_0 = 1/1.0^2 = 1.0`. Data precision per observation: :math:`1/\sigma^2 = 4.0`. After :math:`n = 10` observations: .. math:: \tau_{10} = 1.0 + 10 \times 4.0 = 41.0, \qquad \sigma_{10}^2 = 1/41.0 \approx 0.02439 The posterior mean is a precision-weighted average: .. math:: \mu_{10} = \frac{1.0 \times 20.0 + 40.0 \times \bar{y}}{41.0} **Python implementation:** .. code-block:: python import numpy as np from scipy import stats # Reproducible calibration data rng = np.random.default_rng(42) true_bias = 20.3 # sensor reads 0.3 degrees high sigma_known = 0.5 n = 10 data = rng.normal(true_bias, sigma_known, size=n) print(f"True bias: {true_bias:.1f} C") print(f"Sample mean (MLE): {data.mean():.4f} C") print(f"n={n} readings: {np.round(data, 3)}") print() # Initialize prior model_prior = NormalNormalModel(mu0=20.0, sigma2_prior=1.0, sigma2_known=sigma_known**2) # Batch update model_batch = model_prior.update(data) # Sequential update (should match batch) model_seq = model_prior for x in data: model_seq = model_seq.update_single(x) print("Batch update:") print(f" Posterior mean: {model_batch.posterior_mean:.4f}") print(f" Posterior std: {model_batch.posterior_std:.4f}") lo, hi = model_batch.credible_interval(0.95) print(f" 95% CrI: ({lo:.4f}, {hi:.4f})") print() print("Sequential update (should match):") print(f" Posterior mean: {model_seq.posterior_mean:.4f}") print(f" Posterior std: {model_seq.posterior_std:.4f}") print() # Frequentist comparison freq_mean = data.mean() freq_se = sigma_known / np.sqrt(n) freq_ci = (freq_mean - 1.96 * freq_se, freq_mean + 1.96 * freq_se) print(f"Frequentist 95% CI: ({freq_ci[0]:.4f}, {freq_ci[1]:.4f})") # Prior influence: how much did the prior pull the estimate? print(f"\nPrior weight: {1/(1+10*4):.4f}") print(f"Data weight: {40/41:.4f}") .. code-block:: text True bias: 20.3 C Sample mean (MLE): 20.1322 C n=10 readings: [20.248 20.144 20.493 20.405 20.362 20.307 20.247 20.436 20.249 20.427] Batch update: Posterior mean: 20.1290 Posterior std: 0.1562 95% CrI: (19.8229, 20.4351) Sequential update (should match): Posterior mean: 20.1290 Posterior std: 0.1562 Frequentist 95% CI: (19.8223, 20.4421) Prior weight: 0.0244 Data weight: 0.9756 **Result:** With :math:`n = 10` precise measurements (precision 4.0 each), the data contribute 97.6% of the posterior precision, and the prior has negligible influence. The posterior mean (20.13) is barely distinguishable from the MLE (20.13) β€” a concrete demonstration that even a moderately informative prior is overwhelmed when the data-to-prior precision ratio is large. The 95% Bayesian credible interval (19.82, 20.44) and the frequentist CI (19.82, 20.44) are nearly identical; this agreement is an instance of the Bernstein–von Mises theorem operating in a single-parameter Normal model. Note that this particular seed produced a sample mean (20.13) below the true bias (20.3) β€” a normal occurrence with only 10 observations; the posterior correctly reflects this sample's evidence rather than the true value. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig10_normal_sequential.png :alt: Multi-panel figure showing the evolution of the Normal posterior as temperature calibration readings arrive sequentially; prior (dashed) and successive posterior densities (progressively darker) plotted on a shared axis, with posterior mean tracking toward true value and posterior std shrinking as 1/sqrt(t) :align: center :width: 100% **Figure 5.10**: Sequential Bayesian updating for the Normal-Normal model. Each curve shows the posterior after processing the indicated number of calibration observations (prior shown as a broad dashed curve). Precision accumulates linearly: after :math:`t` observations the posterior standard deviation shrinks as :math:`(\tau_0 + t/\sigma^2)^{-1/2}`. With this particular draw (seed 42, :math:`\bar{y} = 20.13`), the posterior mean tracks the sample evidence rather than the true bias (20.3Β°C) β€” the natural outcome with only 10 observations. The prior center (20.0Β°C) has negligible influence by :math:`t = 10` because the data precision (:math:`10/\sigma^2 = 40`) dominates the prior precision (:math:`\tau_0 = 1`). .. admonition:: PyMC 5 Implementation: Normal-Normal Model :class: note The Normal-Normal model with known :math:`\sigma^2` maps directly to PyMC. The key difference from the analytical version is that PyMC does not require knowing :math:`\sigma^2` in advance β€” replacing :code:`pm.Data` for sigma with :code:`pm.HalfNormal` would immediately generalise this to the NIG model of the next section. .. code-block:: python import pymc as pm import arviz as az import numpy as np # Temperature calibration data (seed 42, sigma known = 0.5) rng = np.random.default_rng(42) sigma_known = 0.5 y_obs = rng.normal(loc=20.3, scale=sigma_known, size=10) # Prior: mu ~ N(20.0, 1.0^2) mu0, tau0 = 20.0, 1.0 with pm.Model() as normal_model: mu = pm.Normal('mu', mu=mu0, sigma=tau0) # Likelihood: y_i ~ N(mu, sigma^2) with sigma known y = pm.Normal('y', mu=mu, sigma=sigma_known, observed=y_obs) idata = pm.sample( draws=4000, tune=1000, chains=4, random_seed=42, progressbar=False, ) # ── Compare MCMC to analytical ────────────────────────────────────── # Analytical posterior parameters (precision-weighted update) tau0_sq = tau0 ** 2 sig_sq = sigma_known ** 2 n = len(y_obs) ybar = y_obs.mean() post_prec = 1 / tau0_sq + n / sig_sq post_var = 1 / post_prec post_mean = post_var * (mu0 / tau0_sq + n * ybar / sig_sq) post_sd = np.sqrt(post_var) mu_samples = idata.posterior['mu'].values.flatten() print("Normal-Normal posterior (temperature calibration)") print(f" Analytical: mean={post_mean:.4f}, sd={post_sd:.4f}") print(f" MCMC: mean={mu_samples.mean():.4f}, sd={mu_samples.std():.4f}") print() # ── ArviZ summary with convergence diagnostics ────────────────────── # This is the standard 'health check' table for any MCMC run print(az.summary(idata, var_names=['mu'], hdi_prob=0.95)) **Expected output:** .. code-block:: text Normal-Normal posterior (temperature calibration) Analytical: mean=20.1316, sd=0.1534 MCMC: mean=20.1315, sd=0.1535 mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk ess_tail r_hat mu 20.132 0.154 19.834 20.432 0.0009 0.0007 29748.0 12009.0 1.0 **Reading the diagnostics.** The :code:`az.summary` table is the standard MCMC health check. For this simple one-parameter model, :code:`r_hat = 1.0` (all chains converged to the same distribution) and :code:`ess_bulk` β‰ˆ 29,748 β€” nearly the maximum possible for 16,000 raw draws, confirming rapid decorrelation. In complex posteriors, bulk ESS below 400 per chain or :code:`r_hat` above 1.01 signal convergence problems that require diagnostic investigation (developed in Section 5.5). Comparing columns: :code:`mcse_mean = 0.0009` means the MCMC estimate of the posterior mean is accurate to Β±0.001, which is more than sufficient given the posterior SD of 0.154. .. _sec-normal-inverse-gamma: Normal-Inverse-Gamma: Unknown Mean and Variance ------------------------------------------------- The Normal-Normal model assumed the observation variance :math:`\sigma^2` was known. In practice both :math:`\mu` and :math:`\sigma^2` are unknown, and they cannot be estimated independently β€” the uncertainty in :math:`\sigma^2` inflates the posterior uncertainty in :math:`\mu`. The Normal-Inverse-Gamma (NIG) model handles this joint inference problem and produces the first example in this chapter of a *multiparameter* conjugate analysis. The joint posterior of :math:`(\mu, \sigma^2)` will turn out to have a recognizable structure: the marginal posterior for :math:`\mu` integrating out :math:`\sigma^2` is a *Student-t distribution*, and under a reference prior this t-distribution is very close to the frequentist t-distribution with :math:`n - 1` degrees of freedom. The NIG model with :math:`\kappa_0, \nu_0 \to 0` gives a marginal t with :math:`\nu_n = n` degrees of freedom; the frequentist t-interval uses :math:`n - 1`. The one-df gap closes asymptotically and is negligible for :math:`n \gtrsim 30`. This near-equivalence β€” Bayesian credible interval approximately matching the frequentist confidence interval under a reference prior β€” is one of the most instructive numerical agreements in statistics. Mathematical Framework ~~~~~~~~~~~~~~~~~~~~~~ The Normal-Inverse-Gamma prior on :math:`(\mu, \sigma^2)` is defined hierarchically: .. math:: :label: nig-prior \mu \mid \sigma^2 &\;\sim\; \mathcal{N}\!\left(\mu_0,\; \frac{\sigma^2}{\kappa_0}\right) \\[4pt] \sigma^2 &\;\sim\; \text{Inv-Gamma}\!\left(\frac{\nu_0}{2},\; \frac{\nu_0 \sigma_0^2}{2}\right) The conditional prior on :math:`\mu` given :math:`\sigma^2` is Normal with variance proportional to the unknown :math:`\sigma^2`: if the data are noisy (large :math:`\sigma^2`), we are also less certain about :math:`\mu` a priori. The four hyperparameters have the following interpretations: - :math:`\mu_0` β€” prior mean for :math:`\mu` - :math:`\kappa_0` β€” prior pseudo-observations for :math:`\mu`; larger :math:`\kappa_0` means tighter prior on :math:`\mu` - :math:`\nu_0` β€” prior degrees of freedom; larger :math:`\nu_0` means more concentrated prior on :math:`\sigma^2` - :math:`\sigma_0^2` β€” prior scale for :math:`\sigma^2`; the prior mean of :math:`\sigma^2` is :math:`\nu_0 \sigma_0^2 / (\nu_0 - 2)` for :math:`\nu_0 > 2` The joint prior density is: .. math:: p(\mu, \sigma^2) \propto (\sigma^2)^{-(\nu_0/2 + 3/2)} \exp\!\left(-\frac{\nu_0 \sigma_0^2 + \kappa_0(\mu - \mu_0)^2}{2\sigma^2}\right) Conjugate Update Rule ~~~~~~~~~~~~~~~~~~~~~~ **Full derivation (following BDA3, Β§3.2).** Let :math:`y = (y_1, \ldots, y_n)` be i.i.d. :math:`\mathcal{N}(\mu, \sigma^2)`. Compute the sample mean :math:`\bar{y}` and the sum of squared deviations from the mean :math:`S = \sum_{i=1}^n (y_i - \bar{y})^2`. The joint log-likelihood is: .. math:: \log p(y \mid \mu, \sigma^2) = -\frac{n}{2}\log\sigma^2 - \frac{1}{2\sigma^2}\left[S + n(\bar{y} - \mu)^2\right] Multiplying the likelihood by the NIG prior and completing the square in :math:`\mu` yields a posterior that is again Normal-Inverse-Gamma with updated hyperparameters: .. admonition:: Definition: Normal-Inverse-Gamma Conjugate Update :class: note If :math:`(\mu, \sigma^2) \sim \text{NIG}(\mu_0, \kappa_0, \nu_0, \sigma_0^2)` and :math:`y_1, \ldots, y_n \mid \mu, \sigma^2 \stackrel{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^2)`, then the posterior is: .. math:: :label: nig-posterior (\mu, \sigma^2) \mid y \;\sim\; \text{NIG}(\mu_n, \kappa_n, \nu_n, \sigma_n^2) with updated hyperparameters: .. math:: :label: nig-update \kappa_n &= \kappa_0 + n \\[4pt] \nu_n &= \nu_0 + n \\[4pt] \mu_n &= \frac{\kappa_0 \mu_0 + n \bar{y}}{\kappa_n} \\[4pt] \nu_n \sigma_n^2 &= \nu_0 \sigma_0^2 + S + \frac{\kappa_0 n}{\kappa_n}(\bar{y} - \mu_0)^2 The last equation has a natural reading: the posterior sum of squares is the prior sum of squares (:math:`\nu_0 \sigma_0^2`) plus the sample sum of squares (:math:`S`) plus a correction for the discrepancy between the sample mean :math:`\bar{y}` and the prior mean :math:`\mu_0`. When :math:`\bar{y} \approx \mu_0`, the correction term vanishes. When the data mean differs substantially from the prior mean, additional uncertainty is accumulated. Marginal Posteriors ~~~~~~~~~~~~~~~~~~~~ The joint posterior integrates to yield two important marginal distributions. **Marginal for** :math:`\sigma^2`. Integrating the joint NIG posterior over :math:`\mu` yields: .. math:: :label: nig-marginal-sigma \sigma^2 \mid y \;\sim\; \text{Inv-Gamma}\!\left(\frac{\nu_n}{2},\; \frac{\nu_n \sigma_n^2}{2}\right) The posterior mean of :math:`\sigma^2` is :math:`\nu_n \sigma_n^2 / (\nu_n - 2)` for :math:`\nu_n > 2`. **Marginal for** :math:`\mu`. Integrating out :math:`\sigma^2` from the joint NIG posterior yields a *Student-t distribution*: .. math:: :label: nig-marginal-mu \mu \mid y \;\sim\; t_{\nu_n}\!\left(\mu_n,\; \frac{\sigma_n^2}{\kappa_n}\right) where :math:`t_{\nu_n}(\mu_n, s^2)` denotes a t-distribution with :math:`\nu_n` degrees of freedom, location :math:`\mu_n`, and scale squared :math:`s^2 = \sigma_n^2 / \kappa_n`. **The frequentist connection.** Under the reference prior :math:`\kappa_0 \to 0`, :math:`\nu_0 \to 0`, :math:`\mu_0` arbitrary, the update rules give: .. math:: \kappa_n = n, \quad \nu_n = n, \quad \mu_n = \bar{y}, \quad \sigma_n^2 = \frac{S}{n-1} \cdot \frac{n-1}{n} = S / n so the marginal posterior for :math:`\mu` becomes: .. math:: :label: nig-frequentist-match \mu \mid y \;\sim\; t_{n}\!\left(\bar{y},\; \frac{S/n}{n}\right) = t_{n}\!\left(\bar{y},\; \frac{S}{n^2}\right) Equivalently, the standardized quantity :math:`(\mu - \bar{y}) / (s / \sqrt{n})` follows a :math:`t_{n-1}` distribution β€” the frequentist t-distribution used to construct confidence intervals. The exact equivalence requires the improper prior :math:`p(\mu, \sigma^2) \propto \sigma^{-2}`, which formally sits at the boundary of the NIG parameterization (:math:`\nu_0 \to -1`). The NIG with :math:`\nu_0 \to 0` gives :math:`\nu_n = n` (not :math:`n - 1`), producing a t-distribution one degree of freedom wider than the frequentist interval. For :math:`n \geq 15` the difference in 95% interval bounds is smaller than :math:`0.02` standard errors and is negligible in practice. Under the reference prior, Bayesian and frequentist inference are asymptotically equivalent for the Normal location-scale model β€” a concrete instance of the Bernstein–von Mises theorem (see :ref:`ch3_3-sampling-variability`). Python Implementation ~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python import numpy as np from scipy import stats from dataclasses import dataclass, field @dataclass class NormalInverseGammaModel: """ Conjugate Normal-Inverse-Gamma model for joint inference about Normal mean mu and variance sigma^2. Hyperparameters follow BDA3 Β§3.2 notation. Parameters ---------- mu0 : float Prior mean for mu. kappa0 : float Prior pseudo-observations for mu (> 0). nu0 : float Prior degrees of freedom for sigma^2 (> 0). sigma2_0 : float Prior scale for sigma^2 (> 0). """ mu0: float kappa0: float nu0: float sigma2_0: float def __post_init__(self) -> None: self._mu_n = self.mu0 self._kappa_n = self.kappa0 self._nu_n = self.nu0 self._sigma2_n = self.sigma2_0 self._n = 0 # ------------------------------------------------------------------ # # Posterior updating # # ------------------------------------------------------------------ # def update(self, data: np.ndarray) -> "NormalInverseGammaModel": """ Return updated model (immutable) after observing data. Parameters ---------- data : np.ndarray 1D array of observations. Returns ------- NormalInverseGammaModel Updated model with posterior hyperparameters. """ n = len(data) ybar = data.mean() S = np.sum((data - ybar) ** 2) new = NormalInverseGammaModel(self.mu0, self.kappa0, self.nu0, self.sigma2_0) # Update rules from eq. (nig-update) new._kappa_n = self._kappa_n + n new._nu_n = self._nu_n + n new._mu_n = (self._kappa_n * self._mu_n + n * ybar) / new._kappa_n # Posterior sum of squares (scale accumulation) correction = (self._kappa_n * n / new._kappa_n) * (ybar - self._mu_n)**2 new._sigma2_n = (self._nu_n * self._sigma2_n + S + correction) / new._nu_n new._n = self._n + n return new # ------------------------------------------------------------------ # # Marginal posteriors # # ------------------------------------------------------------------ # def marginal_mu(self) -> stats.rv_continuous: """ Marginal posterior for mu: t_{nu_n}(mu_n, sigma2_n / kappa_n). Returns a frozen scipy.stats.t object. """ scale = np.sqrt(self._sigma2_n / self._kappa_n) return stats.t(df=self._nu_n, loc=self._mu_n, scale=scale) def marginal_sigma2(self) -> stats.rv_continuous: """ Marginal posterior for sigma^2: Inv-Gamma(nu_n/2, nu_n*sigma2_n/2). Returns a frozen scipy.stats.invgamma object (scale parameterization). """ a = self._nu_n / 2.0 scale = self._nu_n * self._sigma2_n / 2.0 return stats.invgamma(a=a, scale=scale) def credible_interval_mu(self, level: float = 0.95 ) -> tuple[float, float]: """Equal-tailed credible interval for mu using t marginal.""" t_dist = self.marginal_mu() alpha = 1 - level return t_dist.ppf(alpha / 2), t_dist.ppf(1 - alpha / 2) # ------------------------------------------------------------------ # # Posterior summaries # # ------------------------------------------------------------------ # @property def posterior_mu_mean(self) -> float: return self._mu_n @property def posterior_sigma2_mean(self) -> float | None: if self._nu_n > 2: return self._nu_n * self._sigma2_n / (self._nu_n - 2) return None # mean undefined for nu_n <= 2 .. admonition:: Example πŸ’‘ Normal-Inverse-Gamma: Bayesian t-Interval Verification :class: note **Given:** A metrologist measures a physical constant (calibrated weight in grams) with :math:`n = 15` independent measurements. The true value is unknown. A reference (flat) prior is used: :math:`\kappa_0 = 0.001 \approx 0`, :math:`\nu_0 = 0.001 \approx 0`, :math:`\mu_0 = 0` (arbitrary under flat prior), :math:`\sigma_0^2 = 1` (arbitrary). **Find:** The posterior for :math:`\mu` under the NIG model, and verify numerically that the Bayesian 95% credible interval matches the frequentist t-interval to four decimal places. **Python implementation:** .. code-block:: python import numpy as np from scipy import stats # Simulate data: 15 measurements from N(100.3, 0.5^2) rng = np.random.default_rng(2024) true_mu, true_sigma = 100.3, 0.5 n = 15 y = rng.normal(true_mu, true_sigma, size=n) print(f"n={n}, y-bar={y.mean():.4f}, s={y.std(ddof=1):.4f}") print() # --- Bayesian NIG with reference (flat) prior --- # kappa0 -> 0, nu0 -> 0 model = NormalInverseGammaModel( mu0=0.0, kappa0=1e-6, nu0=1e-6, sigma2_0=1.0 ).update(y) t_dist = model.marginal_mu() bayes_lo, bayes_hi = model.credible_interval_mu(0.95) print(f"Bayesian NIG (reference prior):") print(f" Posterior mean of mu: {model.posterior_mu_mean:.4f}") print(f" Posterior t df: {model.marginal_mu().kwds['df']:.2f}") print(f" 95% Credible Interval: ({bayes_lo:.4f}, {bayes_hi:.4f})") print() # --- Frequentist t-interval --- ybar = y.mean() s = y.std(ddof=1) se = s / np.sqrt(n) t_crit = stats.t.ppf(0.975, df=n - 1) freq_lo = ybar - t_crit * se freq_hi = ybar + t_crit * se print(f"Frequentist t-interval:") print(f" Sample mean: {ybar:.4f}") print(f" t critical (df={n-1}): {t_crit:.4f}") print(f" 95% Confidence Interval: ({freq_lo:.4f}, {freq_hi:.4f})") print() # Difference print(f"Difference (Bayes - Freq):") print(f" Lower bound: {bayes_lo - freq_lo:.6f}") print(f" Upper bound: {bayes_hi - freq_hi:.6f}") # --- Verify: grid_posterior_2d from Section 5.1 --- from scipy.special import logsumexp mu_grid = np.linspace(ybar - 4 * se, ybar + 4 * se, 300) sig_grid = np.linspace(0.1, 2.0, 300) MU, SIG = np.meshgrid(mu_grid, sig_grid, indexing='ij') # Identity: sum_i (y_i - mu)^2 = S + n*(ybar - mu)^2 # where S = sum_i (y_i - ybar)^2. # This avoids allocating a (300, 300, n) array β€” a memory bomb # for large n. Compute S once; use it everywhere on the grid. S_stat = np.sum((y - ybar) ** 2) SS = S_stat + n * (ybar - MU) ** 2 # shape (300, 300) β€” O(N_grid) only log_lik = -n * np.log(SIG) - 0.5 * SS / SIG**2 log_prior = -np.log(SIG) # Jeffreys scale prior log_unnorm = log_lik + log_prior log_post = log_unnorm - logsumexp(log_unnorm.ravel()) joint = np.exp(log_post) joint /= joint.sum() marginal_mu_grid = joint.sum(axis=1) marginal_mu_grid /= marginal_mu_grid.sum() * (mu_grid[1] - mu_grid[0]) grid_mean = np.sum(mu_grid * marginal_mu_grid * (mu_grid[1] - mu_grid[0])) print(f"\nGrid posterior mean of mu: {grid_mean:.4f}") print(f"Analytical NIG mean: {model.posterior_mu_mean:.4f}") .. code-block:: text n=15, y-bar=100.4928, s=0.5167 Bayesian NIG (reference prior): Posterior mean of mu: 100.4928 Posterior t df: 15.00 95% Credible Interval: (100.2177, 100.7679) Frequentist t-interval: Sample mean: 100.4928 t critical (df=14): 2.1448 95% Confidence Interval: (100.2067, 100.7789) Difference (Bayes - Freq): Lower bound: 1.10e-02 Upper bound: -1.10e-02 Grid posterior mean of mu: 100.4928 Analytical NIG mean: 100.4928 **Result:** The Bayesian NIG credible interval and the frequentist t-interval agree to two decimal places but not exactly. The systematic gap of :math:`1.1 \times 10^{-2}` on each bound arises from a one-degree-of-freedom difference: the NIG model with :math:`\kappa_0, \nu_0 \to 0` yields a marginal t-distribution with :math:`\nu_n = n = 15` degrees of freedom, while the frequentist t-interval uses :math:`n - 1 = 14`. The asymptotic equivalence requires the improper prior :math:`p(\mu, \sigma^2) \propto \sigma^{-2}`, which formally corresponds to :math:`\nu_0 = -1` β€” outside the valid NIG parameterization. For practical purposes, the difference vanishes as :math:`n \to \infty` and is negligible above :math:`n \approx 50`. The grid posterior mean (100.4928) matches the analytical NIG mean to four decimal places, confirming both computations. .. admonition:: Common Pitfall ⚠️ :class: warning A natural instinct for computing :math:`\sum_i (y_i - \mu)^2` on a 2D grid is to broadcast the raw data array against the grid: .. code-block:: python # AVOID: allocates a (N_grid, N_grid, n) array SS = np.sum((y[np.newaxis, np.newaxis, :] - MU[:, :, np.newaxis])**2, axis=2) For a :math:`300 \times 300` grid and :math:`n = 100{,}000` data points this allocates approximately **72 GB** β€” a silent kernel crash. The Normal likelihood depends on the data only through the sufficient statistics :math:`\bar{y}` and :math:`S = \sum_i (y_i - \bar{y})^2`, so use the algebraic identity instead: .. code-block:: python # CORRECT: O(N_grid) memory, computed once S_stat = np.sum((y - ybar) ** 2) # scalar, O(n) once SS = S_stat + n * (ybar - MU) ** 2 # shape (N_grid, N_grid) The math tells you what to cache. Whenever a likelihood belongs to an exponential family, the grid evaluation should loop over the sufficient statistic, not the raw data. This principle applies equally to the Binomial (:math:`k, n-k`), Poisson (:math:`\sum y_i`), and Gamma (:math:`\sum y_i, \sum \log y_i`) models. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig11_nig_joint_posterior.png :alt: Three-panel figure: (left) 2D joint posterior contours for mu and sigma^2 under the NIG model, with marginals projected onto each axis; (center) marginal posterior for mu compared with the frequentist t-distribution -- the two curves are visually indistinguishable; (right) marginal posterior for sigma^2 compared with the Inverse-Gamma analytical formula and the chi-squared distribution used in frequentist variance inference :align: center :width: 100% **Figure 5.11**: Normal-Inverse-Gamma joint posterior for the weight measurement example (:math:`n=15`). Left: the joint posterior :math:`p(\mu, \sigma^2 \mid y)` as filled contours, with marginal distributions projected onto the horizontal (:math:`\mu`) and vertical (:math:`\sigma^2`) axes. Center: the marginal posterior for :math:`\mu` (solid blue) overlaid on the frequentist :math:`t_{14}` distribution (dashed red). The two curves are nearly identical but differ slightly because the NIG reference prior yields :math:`t_{15}` (df :math:`= n`), whereas the frequentist pivot uses :math:`t_{14}` (df :math:`= n-1`); the gap is negligible for :math:`n \geq 30`. Right: the marginal Inverse-Gamma posterior for :math:`\sigma^2`. The right-skew of the variance posterior is evident, reflecting the asymmetry of the Inverse-Gamma distribution; the posterior mean (0.288) lies above the MLE (0.249) because integration over the skewed distribution pulls the mean rightward. .. _sec-poisson-gamma: Poisson-Gamma Model -------------------- The Poisson-Gamma conjugate pair provides the natural model for count data: event rates, arrival processes, rare-disease incidence, and any application where the data are non-negative integers generated by independent Poisson processes. Its structure mirrors the Beta-Binomial model, with the Gamma distribution playing the role of the Beta. Mathematical Framework ~~~~~~~~~~~~~~~~~~~~~~ Let :math:`y_1, \ldots, y_n` be i.i.d. counts from :math:`\text{Poisson}(\lambda)`, where :math:`\lambda > 0` is the unknown rate. The joint likelihood depends on the data only through the sufficient statistic :math:`T = \sum_{i=1}^n y_i`: .. math:: p(y_1, \ldots, y_n \mid \lambda) \propto \lambda^T e^{-n\lambda}, \qquad T = \sum_{i=1}^n y_i Place a :math:`\text{Gamma}(\alpha_0, \beta_0)` prior on :math:`\lambda`, using the rate parameterization (:math:`\beta_0` is the rate, so the density is :math:`p(\lambda) \propto \lambda^{\alpha_0 - 1} e^{-\beta_0 \lambda}`): .. math:: \lambda \sim \text{Gamma}(\alpha_0, \beta_0), \qquad \mathbb{E}[\lambda] = \frac{\alpha_0}{\beta_0}, \quad \text{Var}(\lambda) = \frac{\alpha_0}{\beta_0^2} **Full derivation:** .. math:: p(\lambda \mid y) \propto \lambda^T e^{-n\lambda} \cdot \lambda^{\alpha_0 - 1} e^{-\beta_0 \lambda} = \lambda^{(\alpha_0 + T) - 1} e^{-(\beta_0 + n)\lambda} This is the kernel of :math:`\text{Gamma}(\alpha_0 + T, \beta_0 + n)`. .. admonition:: Definition: Poisson-Gamma Conjugate Update :class: note If :math:`\lambda \sim \text{Gamma}(\alpha_0, \beta_0)` and :math:`y_1, \ldots, y_n \mid \lambda \stackrel{\text{iid}}{\sim} \text{Poisson}(\lambda)`, then: .. math:: :label: poisson-gamma-posterior \lambda \mid y \;\sim\; \text{Gamma}(\alpha_0 + T,\; \beta_0 + n), \qquad T = \textstyle\sum_{i=1}^n y_i The prior parameter :math:`\alpha_0` accumulates total observed events; :math:`\beta_0` accumulates total observed exposure (number of periods). The posterior mean: .. math:: :label: poisson-gamma-mean \mathbb{E}[\lambda \mid y] = \frac{\alpha_0 + T}{\beta_0 + n} = \frac{\beta_0}{\beta_0 + n}\cdot\frac{\alpha_0}{\beta_0} + \frac{n}{\beta_0 + n}\cdot\frac{T}{n} is again a precision-weighted average of the prior mean :math:`\alpha_0/\beta_0` and the MLE :math:`T/n`, with data weight :math:`w = n/(n + \beta_0)`. The pseudo-data interpretation: :math:`\alpha_0` is the number of prior events and :math:`\beta_0` is the equivalent prior exposure. A Gamma(3, 1) prior says "before seeing data, I believe the rate is consistent with having observed 3 events over 1 unit of exposure." Posterior Predictive: The Negative Binomial ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The posterior predictive for a new count :math:`\tilde{y}` integrates over the posterior: .. math:: p(\tilde{y} \mid y) = \int_0^\infty \frac{\lambda^{\tilde{y}} e^{-\lambda}}{\tilde{y}!} \cdot \frac{(\beta_0+n)^{\alpha_n}}{\Gamma(\alpha_n)} \lambda^{\alpha_n - 1} e^{-(\beta_0+n)\lambda} \, d\lambda Evaluating this integral yields the *Negative Binomial distribution*: .. math:: :label: poisson-predictive-negbin \tilde{y} \mid y \;\sim\; \text{NegBin}\!\left(\alpha_n,\; \frac{\beta_n}{\beta_n + 1}\right) where :math:`\alpha_n = \alpha_0 + T` and :math:`\beta_n = \beta_0 + n`. This is a Gamma mixture of Poissons: the additional uncertainty in :math:`\lambda` induces overdispersion relative to the Poisson. The predictive mean equals the posterior mean :math:`\alpha_n / \beta_n`, but the predictive variance is :math:`\alpha_n(\beta_n + 1)/\beta_n^2 > \alpha_n / \beta_n` β€” always larger than the Poisson variance equal to the mean. Python Implementation ~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python import numpy as np from scipy import stats from dataclasses import dataclass @dataclass class PoissonGammaModel: """ Conjugate Poisson-Gamma model for rate inference. Uses the rate parameterization of the Gamma distribution: p(lambda) proportional to lambda^(alpha-1) * exp(-beta * lambda). Parameters ---------- alpha0 : float Prior shape (> 0); interpretable as alpha0 prior events. beta0 : float Prior rate (> 0); interpretable as beta0 prior exposure. """ alpha0: float beta0: float def __post_init__(self) -> None: self._alpha_n = self.alpha0 self._beta_n = self.beta0 def update(self, data: np.ndarray) -> "PoissonGammaModel": """ Return updated model after observing count data. Parameters ---------- data : np.ndarray 1D array of non-negative integer counts. Returns ------- PoissonGammaModel Updated model. """ new = PoissonGammaModel(self.alpha0, self.beta0) new._alpha_n = self._alpha_n + data.sum() new._beta_n = self._beta_n + len(data) return new @property def posterior_mean(self) -> float: return self._alpha_n / self._beta_n @property def posterior_mode(self) -> float | None: if self._alpha_n >= 1: return (self._alpha_n - 1) / self._beta_n return None # Mode at 0 when alpha < 1 @property def posterior_variance(self) -> float: return self._alpha_n / self._beta_n ** 2 def credible_interval(self, level: float = 0.95) -> tuple[float, float]: """Equal-tailed credible interval via Gamma quantile function.""" dist = stats.gamma(a=self._alpha_n, scale=1.0 / self._beta_n) alpha = 1 - level return dist.ppf(alpha / 2), dist.ppf(1 - alpha / 2) def posterior_predictive(self) -> stats.rv_discrete: """ Returns a frozen NegBin distribution for predictive count of a new single Poisson observation. Note: scipy NegBin uses parameterization P(X=k) proportional to C(k+n-1,k) * p^n * (1-p)^k where n=alpha_n, p = beta_n/(beta_n+1). """ n_nb = self._alpha_n p_nb = self._beta_n / (self._beta_n + 1) return stats.nbinom(n=n_nb, p=p_nb) .. admonition:: Example πŸ’‘ Hospital Rare Condition Admissions :class: note **Given:** A small regional hospital tracks monthly admissions for a rare neurological condition. Historical records from a neighboring hospital suggest the rate is about 3 admissions per month, with modest uncertainty. This is encoded as a Gamma(3, 1) prior (3 pseudo-events over 1 pseudo-month), giving prior mean :math:`\lambda_0 = 3.0`. Over 12 months of observation, the counts are: :math:`y = (2, 4, 1, 5, 3, 3, 2, 4, 3, 1, 5, 2)` for a total :math:`T = 35`. **Find:** The posterior distribution of the monthly rate :math:`\lambda`, and estimate the probability of exceeding the hospital's capacity limit of 6 admissions next month. **Mathematical approach:** Applying the update rule :eq:`poisson-gamma-posterior`: .. math:: \lambda \mid y \;\sim\; \text{Gamma}(3 + 35,\; 1 + 12) = \text{Gamma}(38, 13) Posterior mean: :math:`38/13 \approx 2.923`. The MLE is :math:`T/n = 35/12 \approx 2.917`. With :math:`n_0 = \beta_0 = 1` pseudo-period and :math:`n = 12` actual months, the data weight is :math:`12/13 \approx 0.923`, so the posterior mean barely differs from the MLE. For the capacity probability, we compute :math:`P(\tilde{y} \geq 7 \mid y)` using the Negative Binomial predictive, subtracting from 1 the probability of observing 6 or fewer admissions: .. code-block:: python import numpy as np from scipy import stats # Data y = np.array([2, 4, 1, 5, 3, 3, 2, 4, 3, 1, 5, 2]) n, T = len(y), y.sum() print(f"n={n}, T={T}, MLE={T/n:.4f}") # Prior and posterior alpha0, beta0 = 3.0, 1.0 model = PoissonGammaModel(alpha0, beta0).update(y) print(f"\nPosterior: Gamma({model._alpha_n:.0f}, {model._beta_n:.0f})") print(f"Posterior mean: {model.posterior_mean:.4f}") lo, hi = model.credible_interval(0.95) print(f"95% CrI for lambda: ({lo:.3f}, {hi:.3f})") # Posterior predictive: P(y_new >= 7) nb_dist = model.posterior_predictive() p_exceed = 1 - nb_dist.cdf(6) print(f"\nP(new month >= 7 admissions): {p_exceed:.4f}") # Compare: plug-in Poisson (ignores uncertainty in lambda) mle = T / n p_poisson_plugin = 1 - stats.poisson.cdf(6, mu=mle) print(f"P(>= 7) under Poisson plug-in: {p_poisson_plugin:.4f}") .. code-block:: text n=12, T=35, MLE=2.9167 Posterior: Gamma(38, 13) Posterior mean: 2.9231 95% CrI for lambda: (2.069, 3.923) P(new month >= 7 admissions): 0.0352 P(>= 7) under Poisson plug-in: 0.0295 **Result:** The probability of exceeding capacity next month is 3.52% under the Bayesian posterior predictive, compared to 2.95% under the Poisson plug-in. The gap β€” about 20% relative difference β€” arises because the Negative Binomial predictive accounts for uncertainty in :math:`\lambda` itself, producing a heavier right tail than the Poisson. This difference matters operationally: a hospital planning capacity based on the plug-in probability systematically underestimates tail risk, and the divergence grows as prior uncertainty in :math:`\lambda` increases relative to the data. .. admonition:: PyMC 5 Implementation: Poisson-Gamma Model :class: note The Poisson-Gamma model in PyMC uses :code:`pm.Gamma` for the prior on the rate and :code:`pm.Poisson` for the likelihood. Note the parameterisation: PyMC's :code:`pm.Gamma` uses :code:`alpha` (shape) and :code:`beta` (rate), matching the convention used throughout this section. .. code-block:: python import pymc as pm import arviz as az import numpy as np # Hospital admissions: n=12 months, counts y_counts = np.array([1, 4, 2, 3, 5, 2, 3, 4, 2, 3, 5, 1]) n, T = len(y_counts), y_counts.sum() # n=12, T=35 # Prior: Gamma(alpha0=3, beta0=1) β€” prior mean = 3 admissions/month alpha0, beta0 = 3.0, 1.0 with pm.Model() as poisson_model: # Rate parameter: lambda ~ Gamma(alpha0, beta0) # PyMC Gamma: alpha=shape, beta=rate (not scale) lam = pm.Gamma('lambda', alpha=alpha0, beta=beta0) # Likelihood: each monthly count ~ Poisson(lambda) y = pm.Poisson('y', mu=lam, observed=y_counts) # Prior predictive: what does the prior imply about counts? prior_idata = pm.sample_prior_predictive(samples=2000, random_seed=42) # Posterior idata = pm.sample( draws=4000, tune=1000, chains=4, random_seed=42, progressbar=False, ) # Posterior predictive: distribution of next month's count pm.sample_posterior_predictive( idata, extend_inferencedata=True, random_seed=42 ) # ── Analytical vs MCMC ────────────────────────────────────────────── alpha_n = alpha0 + T # = 38 beta_n = beta0 + n # = 13 lam_samples = idata.posterior['lambda'].values.flatten() y_new_samples = idata.posterior_predictive['y'].values.flatten() print(f"Analytical posterior: Gamma({alpha_n}, {beta_n})") print(f" mean = {alpha_n/beta_n:.4f}, sd = {(alpha_n/beta_n**2)**0.5:.4f}") print(f"MCMC posterior mean: {lam_samples.mean():.4f}") print() print(az.summary(idata, var_names=['lambda'], hdi_prob=0.95)) print() print(f"P(next month >= 7 | data) [MCMC]: {(y_new_samples >= 7).mean():.4f}") .. code-block:: text Analytical posterior: Gamma(38, 13) mean = 2.9231, sd = 0.4742 mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk ess_tail r_hat lambda 2.923 0.473 2.007 3.869 0.003 0.002 21012.0 11981.0 1.0 P(next month >= 7 | data) [MCMC]: 0.0355 The MCMC posterior mean (2.923) matches the analytical Gamma(38, 13) mean exactly, and the tail probability from posterior predictive samples (3.55%) agrees closely with the Negative Binomial calculation (3.52%). The slight difference arises from Monte Carlo variance in the posterior predictive samples. .. _sec-multinomial-dirichlet: Multinomial-Dirichlet Model ----------------------------- The Beta-Binomial model handles two outcomes (success/failure). Many real applications involve :math:`K > 2` categories: voting intentions in a :math:`K`-candidate election, frequencies of :math:`K` different DNA nucleotides, topic proportions in a text corpus, or multinomial response categories in a survey. The Multinomial-Dirichlet model extends Beta-Binomial conjugacy to this general case. Mathematical Framework ~~~~~~~~~~~~~~~~~~~~~~ Let :math:`n_1, \ldots, n_K` be counts of observations in :math:`K` categories, with :math:`\sum_k n_k = n`. The likelihood under a Multinomial(:math:`n, \theta`) model is: .. math:: p(n_1, \ldots, n_K \mid \theta) \propto \prod_{k=1}^K \theta_k^{n_k}, \qquad \theta_k \geq 0, \quad \sum_k \theta_k = 1 The parameter vector :math:`\theta = (\theta_1, \ldots, \theta_K)` lives on the :math:`(K-1)`-dimensional simplex. The conjugate prior is the Dirichlet distribution: .. math:: :label: dirichlet-prior p(\theta \mid \alpha) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)} \prod_{k=1}^K \theta_k^{\alpha_k - 1}, \qquad \alpha_k > 0 where :math:`\alpha = (\alpha_1, \ldots, \alpha_K)` are the concentration parameters. The Dirichlet is the multivariate generalization of the Beta: when :math:`K=2`, :math:`\text{Dir}(\alpha_1, \alpha_2) = \text{Beta}(\alpha_1, \alpha_2)`. Each :math:`\alpha_k` plays the role of a pseudo-count for category :math:`k`, with prior mean :math:`\mathbb{E}[\theta_k] = \alpha_k / \sum_{j} \alpha_j`. **Conjugate update rule:** Multiplying the Multinomial likelihood by the Dirichlet prior and collecting powers of :math:`\theta_k`: .. math:: p(\theta \mid n) \propto \prod_k \theta_k^{n_k} \cdot \prod_k \theta_k^{\alpha_k - 1} = \prod_k \theta_k^{(\alpha_k + n_k) - 1} .. admonition:: Definition: Multinomial-Dirichlet Conjugate Update :class: note If :math:`\theta \sim \text{Dir}(\alpha_1, \ldots, \alpha_K)` and :math:`(n_1, \ldots, n_K) \mid \theta \sim \text{Multinomial}(n, \theta)`, then: .. math:: :label: dirichlet-posterior \theta \mid n \;\sim\; \text{Dir}(\alpha_1 + n_1,\; \ldots,\; \alpha_K + n_K) The update rule is element-wise addition of observed counts to prior pseudo-counts β€” the same structure as Beta-Binomial, generalized to :math:`K` categories. The posterior mean for category :math:`k` is: .. math:: \mathbb{E}[\theta_k \mid n] = \frac{\alpha_k + n_k}{\sum_j \alpha_j + n} which is a weighted average of the prior proportion :math:`\alpha_k / \sum_j \alpha_j` and the empirical proportion :math:`n_k / n`, with the prior contributing effective sample size :math:`n_0 = \sum_k \alpha_k`. The marginal distribution of any single component :math:`\theta_k` is Beta: :math:`\theta_k \sim \text{Beta}(\alpha_k, \sum_{j \neq k} \alpha_j)`. This means Beta-Binomial analysis of individual categories is always consistent with a full Dirichlet analysis β€” one can zoom in on any single category without loss. .. admonition:: Example πŸ’‘ Three-Candidate Election Forecasting :class: note **Given:** A pre-election poll surveys :math:`n = 500` voters. Results: Candidate A: 215, Candidate B: 180, Candidate C: 105. A prior from the previous election cycle suggests proportions (0.42, 0.38, 0.20), encoded as :math:`\text{Dir}(42, 38, 20)` β€” 100 pseudo-votes reflecting the last election's proportions. **Find:** The posterior distribution of voter proportions and the posterior probability that Candidate A achieves a majority (more than 50% of the vote). .. code-block:: python import numpy as np from scipy import stats # Observed counts and prior counts = np.array([215, 180, 105]) alpha_prior = np.array([42.0, 38.0, 20.0]) # Conjugate update: element-wise addition alpha_post = alpha_prior + counts n0 = alpha_prior.sum() n = counts.sum() print(f"Prior effective sample size: {n0:.0f}") print(f"Data n: {n}") print(f"Data weight: {n/(n + n0):.3f}") print() # Posterior means post_means = alpha_post / alpha_post.sum() prior_means = alpha_prior / alpha_prior.sum() mle = counts / n labels = ['Candidate A', 'Candidate B', 'Candidate C'] print(f"{'':12} {'Prior':>8} {'MLE':>8} {'Post Mean':>10} {'95% CrI':>22}") print("-" * 65) for k in range(3): # Marginal Beta for each candidate beta_dist = stats.beta(alpha_post[k], alpha_post.sum() - alpha_post[k]) lo, hi = beta_dist.ppf(0.025), beta_dist.ppf(0.975) print(f"{labels[k]:12} {prior_means[k]:>8.3f} {mle[k]:>8.3f} " f"{post_means[k]:>10.3f} ({lo:.3f}, {hi:.3f})") # Monte Carlo: P(theta_A > 0.5) rng = np.random.default_rng(42) samples = rng.dirichlet(alpha_post, size=100_000) p_majority = (samples[:, 0] > 0.5).mean() print(f"\nP(Candidate A > 50% majority): {p_majority:.4f}") .. code-block:: text Prior effective sample size: 100 Data n: 500 Data weight: 0.833 Prior MLE Post Mean 95% CrI ----------------------------------------------------------------- Candidate A 0.420 0.430 0.428 (0.389, 0.468) Candidate B 0.380 0.360 0.363 (0.325, 0.402) Candidate C 0.200 0.210 0.208 (0.177, 0.242) P(Candidate A > 50% majority): 0.0002 **Result:** With 83.3% data weight, the posterior proportions track the MLE closely, with the prior pulling candidate A's estimate slightly toward the historical 42% from the poll's 43%. The posterior credible intervals (e.g., 38.9%–46.8% for Candidate A) give a realistic picture of uncertainty given 500 responses. Importantly, the posterior probability of a majority is only 0.02% β€” despite leading the poll, the 95% credible interval does not come close to the 50% threshold. The Dirichlet model makes this uncertainty explicit in a way that a simple point estimate cannot. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig12_conjugate_catalog.png :alt: Visual summary table of the five conjugate families covered in this section, showing for each: likelihood family, conjugate prior, posterior update rule, posterior mean formula as a weighted average, and pseudo-data interpretation :align: center :width: 100% **Figure 5.12**: The Conjugate Family Catalog. Each row summarizes one conjugate pair: likelihood, prior, update rule, posterior mean as a precision-weighted average of prior mean and MLE, and the pseudo-data interpretation of the prior hyperparameters. The pattern is identical across all five families β€” the differences lie only in the parameterization. The weight :math:`w = n/(n + n_0)` quantifies how much the data dominate the prior, with :math:`n_0` defined as the prior's effective sample size in each family. Beyond Conjugacy: The Full Landscape of Prior Specification ------------------------------------------------------------- The five conjugate families give us tractable closed-form posteriors. But they solve only part of the prior specification problem: they tell us *what posterior a given conjugate prior produces*, not *how to choose the prior hyperparameters* or *what to do when no conjugate prior exists*. This section develops the broader landscape of prior specification: invariance-based Jeffreys priors, weakly informative priors for regularization, informative priors grounded in prior data, and the diagnostic tools for evaluating any prior choice. .. _sec-jeffreys-priors: Jeffreys Priors and the Invariance Principle ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A flat prior :math:`p(\theta) \propto 1` seems like a natural "uninformative" choice, but it has a fundamental deficiency: it is not invariant under reparameterization. If :math:`\theta` has a flat prior, then the transformed parameter :math:`\phi = \log\theta` has prior :math:`p(\phi) \propto e^\phi` β€” not flat. The inference about :math:`\theta` would change if the same analysis were conducted in terms of :math:`\phi`. This is not a minor technical inconvenience; it means that "flat" priors depend on an arbitrary parameterization convention. **Harold Jeffreys** proposed a principled solution in 1939: choose the prior to be proportional to the square root of the Fisher information determinant. This prior is *invariant under reparameterization* β€” it produces the same posterior inference regardless of which parameterization is used. Connection to Fisher Information ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Recall from :ref:`ch3_3-sampling-variability` that the Fisher information for a single observation from a model :math:`p(y \mid \theta)` is: .. math:: I_1(\theta) = -\mathbb{E}\!\left[\frac{\partial^2 \log p(y \mid \theta)}{\partial \theta^2}\right] = \mathbb{E}\!\left[\left(\frac{\partial \log p(y \mid \theta)}{\partial \theta}\right)^2\right] For :math:`n` i.i.d. observations, the total Fisher information is :math:`I(\theta) = n \cdot I_1(\theta)`. .. admonition:: Definition: Jeffreys Prior :class: note The **Jeffreys prior** for a single parameter :math:`\theta` is: .. math:: :label: jeffreys-prior p_J(\theta) \propto \sqrt{I_1(\theta)} For a multiparameter model with parameter vector :math:`\theta = (\theta_1, \ldots, \theta_d)`, the Jeffreys prior uses the Fisher information *matrix*: .. math:: p_J(\theta) \propto \sqrt{\det\!\left[\mathcal{I}(\theta)\right]} where :math:`\mathcal{I}_{ij}(\theta) = -\mathbb{E}[\partial^2 \log p(y \mid \theta) / \partial \theta_i \partial \theta_j]`. The invariance property of Jeffreys priors follows immediately from the chain rule of calculus. If :math:`\phi = g(\theta)` is a smooth reparameterization, then by the information transformation formula :math:`I_1^{(\phi)}(\phi) = I_1^{(\theta)}(\theta) / |g'(\theta)|^2`, and the prior density transforms as :math:`p_J(\phi) = p_J(\theta)|d\theta/d\phi| \propto \sqrt{I_1^{(\theta)}(\theta)} / |g'(\theta)| \cdot |d\theta/d\phi| = \sqrt{I_1^{(\phi)}(\phi)}`. The Jeffreys prior in the :math:`\phi` parameterization is exactly what you would obtain by applying the Jeffreys formula directly to the model in the :math:`\phi` parameterization β€” invariance holds exactly. Derived Jeffreys Priors for Standard Models ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Binomial(:math:`n`, :math:`\theta`).** The Fisher information is :math:`I_1(\theta) = 1/[\theta(1-\theta)]`. Therefore: .. math:: :label: jeffreys-binomial p_J(\theta) \propto \left[\theta(1-\theta)\right]^{-1/2} = \theta^{-1/2}(1-\theta)^{-1/2} This is the kernel of a :math:`\text{Beta}(1/2, 1/2)` distribution β€” exactly the Jeffreys prior used in Section 5.1's "Try It Yourself" exercise and in the drug trial example above. The prior concentrates mass at the boundaries :math:`\theta = 0` and :math:`\theta = 1`, reflecting the Fisher information's blow-up near 0 and 1 (these are the parameter values for which an experiment is most informative per unit of Fisher information). **Normal(:math:`\mu`, :math:`\sigma^2`) with both parameters unknown.** The Fisher information matrix for :math:`(\mu, \sigma^2)` is diagonal with entries :math:`I_{\mu\mu} = 1/\sigma^2` and :math:`I_{\sigma^2\sigma^2} = 1/(2\sigma^4)`. The Jeffreys prior is: .. math:: :label: jeffreys-normal p_J(\mu, \sigma^2) \propto \sqrt{I_{\mu\mu} \cdot I_{\sigma^2\sigma^2}} = \frac{1}{\sigma^3} For inference about :math:`\mu` alone (treating :math:`\sigma^2` as a nuisance), the Jeffreys prior is flat: :math:`p_J(\mu) \propto 1`. For the scale parameter :math:`\sigma` (or :math:`\sigma^2`), the Jeffreys prior is :math:`p_J(\sigma) \propto 1/\sigma`, equivalently :math:`p_J(\sigma^2) \propto 1/\sigma^2`. This is the scale-invariant prior used in Section 5.1's 2D grid example. **Poisson(:math:`\lambda`).** The Fisher information is :math:`I_1(\lambda) = 1/\lambda`, giving: .. math:: :label: jeffreys-poisson p_J(\lambda) \propto \lambda^{-1/2} This is the kernel of :math:`\text{Gamma}(1/2, 0)` β€” an improper distribution that cannot be normalized. The Jeffreys prior for the Poisson rate is improper. This is not necessarily a problem: even though the prior does not integrate to a finite value, the posterior after observing any :math:`n \geq 1` observations is proper. However, improper priors must be checked for posterior properness on a case-by-case basis, and they make the marginal likelihood (Bayes factors) undefined. Python Implementation: Computing Jeffreys Priors Numerically ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: python import numpy as np from scipy import stats def jeffreys_prior(family: str, grid: np.ndarray, **params) -> np.ndarray: """ Compute the Jeffreys prior density numerically on a grid, using the Fisher information formula p_J(theta) propto sqrt(I_1(theta)). Parameters ---------- family : str Distribution family: 'binomial', 'poisson', or 'normal_mu'. grid : np.ndarray 1D array of parameter values at which to evaluate the prior. **params : dict Family-specific parameters (e.g., n=50 for binomial). Returns ------- np.ndarray Unnormalized Jeffreys prior evaluated on grid. For proper priors, normalized to sum to 1 over the grid. """ if family == 'binomial': # I_1(theta) = n / (theta * (1 - theta)); Jeffreys: Beta(0.5, 0.5) n_trials = params.get('n', 1) fisher_info = n_trials / (grid * (1 - grid)) prior = np.sqrt(fisher_info) elif family == 'poisson': # I_1(lambda) = 1 / lambda; Jeffreys: Gamma(0.5, 0) improper fisher_info = 1.0 / grid prior = np.sqrt(fisher_info) # propto lambda^{-0.5} elif family == 'normal_mu': # Fisher info for mu with known sigma: I_1(mu) = 1/sigma^2 (const) sigma = params.get('sigma', 1.0) fisher_info = np.full_like(grid, 1.0 / sigma**2) prior = np.sqrt(fisher_info) # constant -> flat prior elif family == 'normal_sigma': # Fisher info for sigma (scale): I_1(sigma) = 2/sigma^2 fisher_info = 2.0 / grid**2 prior = np.sqrt(fisher_info) # propto 1/sigma else: raise ValueError(f"Unknown family: {family}") # Normalize if integrable (improper priors may overflow or underflow) total = prior.sum() if np.isfinite(total) and total > 0: prior = prior / total return prior # --------------------------------------------------------------- # Verification: Jeffreys prior for Binomial matches Beta(0.5, 0.5) # --------------------------------------------------------------- theta_grid = np.linspace(0.001, 0.999, 1000) jp_numerical = jeffreys_prior('binomial', theta_grid, n=1) jp_analytical = stats.beta(0.5, 0.5).pdf(theta_grid) jp_analytical /= jp_analytical.sum() max_err = np.max(np.abs(jp_numerical - jp_analytical)) print(f"Max error vs Beta(0.5, 0.5): {max_err:.2e}") # Jeffreys prior for Poisson (improper, lambda in [0.1, 20]) lambda_grid = np.linspace(0.1, 20.0, 500) jp_poisson = jeffreys_prior('poisson', lambda_grid) print(f"\nPoisson Jeffreys prior is improper: " f"integral over [0.1, 20] = {jp_poisson.sum():.4f} (not normalized)") print("Posterior with Gamma(0.5, 0) prior + n>=1 obs is Gamma(0.5+T, n), proper.") .. code-block:: text Max error vs Beta(0.5, 0.5): 2.63e-06 Poisson Jeffreys prior is improper: integral over [0.1, 20] = 1.0000 (normalized over this finite grid, but not over all lambda > 0) Posterior with Gamma(0.5, 0) prior + n>=1 obs is Gamma(0.5+T, n), proper. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig15_jeffreys_priors.png :alt: Three-panel figure: (left) Jeffreys prior for Binomial parameter theta, shown as Beta(0.5,0.5) density with the Fisher information curve I(theta) overlaid; (center) Jeffreys prior for Normal scale parameter sigma, proportional to 1/sigma, overlaid with a flat prior for comparison; (right) Jeffreys prior for Poisson rate lambda, proportional to lambda^{-0.5}, illustrating the improper prior concentrating mass near zero :align: center :width: 100% **Figure 5.15**: Jeffreys priors for three standard families, derived from Fisher information via :math:`p_J(\theta) \propto \sqrt{I_1(\theta)}`. Left: for the Binomial proportion, the Fisher information diverges at :math:`\theta = 0` and :math:`\theta = 1`, concentrating the prior at the boundaries β€” the resulting Beta(1/2, 1/2) is a proper prior. Center: for the Normal scale :math:`\sigma`, the Jeffreys prior is :math:`1/\sigma` (a scale-invariant improper prior), substantially different from a flat prior on :math:`\sigma`. Right: for the Poisson rate :math:`\lambda`, the Jeffreys prior :math:`\lambda^{-1/2}` is improper, concentrating near zero; it yields a proper posterior for :math:`n \geq 1`. Limitations of Jeffreys Priors ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Despite their theoretical elegance, Jeffreys priors have practical limitations that have motivated alternative approaches. *Multiparameter problems.* In models with several parameters, the joint Jeffreys prior can behave unexpectedly. For the Normal(:math:`\mu, \sigma^2`) model, the joint Jeffreys prior :math:`p_J(\mu, \sigma) \propto 1/\sigma^3` yields a proper posterior only when :math:`n \geq 3` β€” the same joint prior in a slightly different parameterization could be problematic. More seriously, for hierarchical models with many parameters, joint Jeffreys priors are often inadmissible (dominated by other priors in terms of frequentist risk). The Bernardo–Berger *reference prior* framework generalizes Jeffreys priors to the multiparameter setting by separately treating each parameter given the others, but the full development is beyond this course. *Improper posteriors.* An improper Jeffreys prior does not guarantee a proper posterior. The combination of an improper prior with a likelihood that decays insufficiently fast can yield an improper posterior β€” one that does not normalize to a valid probability distribution. This must be verified analytically before using any improper prior. *The "uninformative" misnomer.* Jeffreys priors are not uninformative β€” they encode the assumption of parameterization invariance, which is itself a substantive choice. They place substantial mass at boundaries and near singularities where the Fisher information is large, which can cause numerical difficulties. Modern Bayesian practice often prefers *weakly informative* priors that are proper and regularizing, even at the cost of strict invariance. .. _sec-weakly-informative: Weakly Informative Priors and Prior Predictive Simulation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The modern consensus in applied Bayesian analysis, articulated most clearly by McElreath (2020) and Gelman et al. (2013, Β§2.9), favors *weakly informative priors*: priors that rule out implausible parameter values without imposing strong opinions about the plausible range. The goal is regularization β€” "help the model be sensible without telling it the answer." What Makes a Prior Weakly Informative? ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A weakly informative prior satisfies two criteria: 1. **It excludes the physically impossible.** A regression coefficient for the effect of age on income cannot be :math:`10^{12}`. A standard deviation cannot be negative. A probability cannot exceed one. The prior should assign zero (or negligible) probability to impossible regions. 2. **It does not concentrate mass in a narrow range of plausible values.** The prior should remain broad within the plausible range, allowing the data to identify the parameter's value. A Normal(0, 1) prior on a regression coefficient whose true value might be anywhere from :math:`-10` to :math:`+10` would be too informative β€” it would shrink coefficients toward zero more aggressively than the data warrant. Standard weakly informative priors: - **Regression coefficients (standardized predictors):** :math:`\mathcal{N}(0, 2.5^2)` (Gelman recommendation for logistic regression) or :math:`\mathcal{N}(0, 1^2)` (McElreath default for standardized data) - **Scale parameters (standard deviations):** Half-Normal or Half-:math:`t` with a modest scale; Half-Cauchy(:math:`0, 2.5`) (Gelman 2006) for hierarchical variance components - **Unbounded location parameters:** :math:`\mathcal{N}(0, \sigma^2)` with :math:`\sigma` chosen to cover the plausible range of the parameter Prior Predictive Simulation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The most powerful tool for evaluating a prior is *prior predictive simulation*: simulate data from the model *before* seeing any observations, and ask whether the simulated datasets look plausible. This operationalizes the criterion "does my prior encode reasonable domain knowledge?" in a concrete, falsifiable way. **Algorithm: Prior Predictive Check** .. code-block:: text Input: Prior p(theta); likelihood p(y | theta, x); covariate values x Output: Prior predictive samples {y^(1), ..., y^(S)} 1. For s = 1, 2, ..., S: a. Draw theta^(s) from p(theta) b. Draw y^(s) from p(y | theta^(s), x) 2. Visualize the distribution of {y^(s)} and ask: - Are there simulated outcomes that are physically impossible? - Does the spread of simulated data match domain expectations? - Would a subject-matter expert find these datasets plausible? **Computational Complexity**: :math:`O(S \cdot n)` β€” embarrassingly parallel. .. code-block:: python import numpy as np from scipy import stats def prior_predictive_sample( prior_sampler, likelihood_sampler, n_obs: int, n_sims: int = 1000, rng: np.random.Generator | None = None ) -> np.ndarray: """ Generate prior predictive samples by composing a prior sampler with a likelihood sampler. The two callables define all model assumptions. To compare prior specifications, define multiple prior_sampler functions and call this function once per prior β€” the likelihood_sampler is unchanged. Parameters ---------- prior_sampler : callable Signature (rng, n_sims) -> params, where params is either a dict mapping parameter names to arrays of shape (n_sims,), or an array of shape (n_sims, d_params). likelihood_sampler : callable Signature (params, rng, n_obs) -> array of shape (n_sims, n_obs). Receives the params structure returned by prior_sampler. n_obs : int Number of observations per simulated dataset (or number of predictor points for deterministic spaghetti-plot use). n_sims : int Number of prior predictive datasets. rng : np.random.Generator, optional For reproducibility. Returns ------- np.ndarray, shape (n_sims, n_obs) Prior predictive datasets. """ if rng is None: rng = np.random.default_rng() params = prior_sampler(rng, n_sims) return likelihood_sampler(params, rng, n_obs) .. admonition:: Example πŸ’‘ Prior Predictive Simulation for Normal Regression :class: note **Given:** A simple linear regression model: .. math:: y_i = \alpha + \beta x_i + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \sigma^2) with predictor :math:`x \in [-2, 2]` (standardized). Three candidate priors for the slope :math:`\beta` are compared via prior predictive simulation: - **Prior A**: :math:`\mathcal{N}(0, 100^2)` β€” superficially "vague" - **Prior B**: :math:`\mathcal{N}(0, 10^2)` β€” moderately diffuse - **Prior C**: :math:`\mathcal{N}(0, 1^2)` β€” weakly informative All use :math:`\alpha \sim \mathcal{N}(0, 10^2)` and :math:`\sigma \sim \text{HalfNormal}(1)`. **Find:** Does each prior generate plausible regression lines? What does "vague" N(0, 100Β²) actually imply about the data? **Python implementation:** .. code-block:: python import numpy as np from scipy import stats rng = np.random.default_rng(42) n_sims = 200 x_grid = np.linspace(-2, 2, 100) # Standardized predictor grid # ── Step 1: define the prior sampler ──────────────────────────────── # Returns a dict of parameter arrays, one entry per parameter. # This is the only place sigma_beta appears β€” change it here and # everything downstream updates automatically. def make_prior_sampler(sigma_beta: float): """Return a prior sampler with slope std sigma_beta locked in.""" def sampler(rng: np.random.Generator, n_sims: int) -> dict: return { 'alpha': rng.normal(0, 10, size=n_sims), 'beta': rng.normal(0, sigma_beta, size=n_sims), 'sigma': np.abs(rng.normal(0, 1, size=n_sims)), # half-Normal } return sampler # ── Step 2: define the likelihood sampler ─────────────────────────── # Returns E[y | x] for each draw β€” the deterministic skeleton of each # regression line. n_obs here is len(x_grid), not number of data rows. def regression_lines( params: dict, rng: np.random.Generator, # unused for E[y] but required by API n_obs: int ) -> np.ndarray: """Shape: (n_sims, n_obs) β€” E[y | x] for each prior draw.""" return (params['alpha'][:, np.newaxis] + params['beta'][:, np.newaxis] * x_grid[np.newaxis, :]) # ── Step 3: call the abstraction for each prior ────────────────────── prior_sigmas = {'A: N(0, 100^2)': 100.0, 'B: N(0, 10^2)': 10.0, 'C: N(0, 1^2)': 1.0} for label, sigma_beta in prior_sigmas.items(): y_lines = prior_predictive_sample( prior_sampler = make_prior_sampler(sigma_beta), likelihood_sampler = regression_lines, n_obs = len(x_grid), n_sims = n_sims, rng = rng, ) # y_lines shape: (n_sims, 100) β€” each row is one prior predictive line y_at_x0 = y_lines[:, x_grid.size // 2] # near x = 0 y_at_x2 = y_lines[:, -1] # at x = 2 # Recover beta draws to compute slope statistics beta_draws = (y_lines[:, -1] - y_lines[:, 0]) / (x_grid[-1] - x_grid[0]) print(f"\nPrior {label}:") print(f" E[y] at x=0: mean={y_at_x0.mean():.2f}, " f"sd={y_at_x0.std():.2f}, " f"range=[{y_at_x0.min():.0f}, {y_at_x0.max():.0f}]") print(f" E[y] at x=2: mean={y_at_x2.mean():.2f}, " f"sd={y_at_x2.std():.2f}, " f"range=[{y_at_x2.min():.0f}, {y_at_x2.max():.0f}]") print(f" |slope| > 10: {(np.abs(beta_draws) > 10).mean():.1%}") print(f" |slope| > 50: {(np.abs(beta_draws) > 50).mean():.1%}") .. code-block:: text Prior A: N(0, 100^2): E[y] at x=0: mean=0.27, sd=9.88, range=[-29, 34] E[y] at x=2: mean=0.47, sd=203.66, range=[-559, 557] |slope| > 10: 92.0% |slope| > 50: 66.5% Prior B: N(0, 10^2): E[y] at x=0: mean=0.08, sd=9.20, range=[-29, 35] E[y] at x=2: mean=0.45, sd=21.76, range=[-56, 61] |slope| > 10: 31.5% |slope| > 50: 0.5% Prior C: N(0, 1^2): E[y] at x=0: mean=-0.03, sd=9.62, range=[-27, 27] E[y] at x=2: mean=0.36, sd=9.68, range=[-24, 29] |slope| > 10: 0.0% |slope| > 50: 0.0% **Result:** The prior predictive simulation reveals that the "vague" Prior A (N(0, 100Β²)) generates regression slopes with absolute values exceeding 10 in 92% of simulations and exceeding 50 in 67% β€” implying regression lines that are nearly vertical for a standardized predictor. If the response variable is measured in units where a one-standard-deviation change in the predictor causes a change of at most a few units in the response, Prior A is not vague but wildly overconfident in the wrong direction β€” it concentrates mass on extreme slopes, leaving almost no probability for the moderate slopes that domain knowledge suggests. Prior C (N(0, 1Β²)) generates lines that span the response range without near-vertical or near-horizontal extremes, which is the intended behavior of a weakly informative prior. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig13_prior_predictive.png :alt: Three-panel spaghetti plot showing 200 prior predictive regression lines drawn from each of the three priors (A, B, C) against the standardized predictor x in [-2, 2]; Prior A shows nearly vertical lines spanning thousands of y-units; Prior B shows moderate scatter; Prior C shows realistic shallow-to-moderate slopes :align: center :width: 100% **Figure 5.13**: Prior predictive simulation for linear regression under three priors for the slope :math:`\beta`. Each panel shows 200 prior predictive regression lines :math:`E[y \mid x] = \alpha + \beta x` drawn from the prior. The "vague" Normal(0, 100Β²) prior (left) generates regression lines that are nearly vertical β€” implying changes of hundreds of :math:`y`-units per unit change in :math:`x`, which is physically implausible for most applications. The Normal(0, 10Β²) prior (center) is better but still allows implausibly steep slopes. The weakly informative Normal(0, 1) prior (right) concentrates predictive lines in a plausible range. This demonstration, following McElreath (2020, Ch. 4), shows that prior predictive simulation is the most effective tool for calibrating prior choices before data collection. .. admonition:: Try It Yourself πŸ–₯️ :class: tip The **Prior Predictive Spaghetti Plot** simulation makes the Οƒ_Ξ² effect tangible. Drag the slope prior standard deviation from 0.1 to 6.0 and watch the prior predictive regression lines update in real time. `Launch: Prior Predictive Simulation `_ Key thresholds to notice: - :math:`\sigma_\beta \approx 0.5`–:math:`1.0`: fewer than 5% of lines are implausibly steep β€” weakly informative territory (green indicator) - :math:`\sigma_\beta \approx 2`–:math:`3`: about 30% of lines are steep β€” moderate but still defensible (amber) - :math:`\sigma_\beta \geq 3.5`: over 60% of prior draws imply near-vertical relationships β€” implausibly vague (red) Use the preset buttons to jump directly to the N(0, 1Β²) versus N(0, 100Β²) comparison discussed in the example above. The "Draw New Sample" button reseeds the random number generator so you can verify the pattern holds across different draws. Arrow keys adjust Οƒ_Ξ² in increments of 0.5; press R to reset. .. admonition:: Common Pitfall ⚠️ :class: warning The most common mistake in prior specification is confusing "vague" with "weakly informative." A Normal(0, 100Β²) prior on a standardized regression coefficient is not vague β€” it assigns 92% of its probability mass to slopes with absolute value greater than 10. Before accepting any prior as "uninformative," run a prior predictive simulation and inspect the implied data. If the simulated datasets include outcomes you would describe as "obviously wrong" or "physically impossible," the prior is informative in a harmful way, regardless of how it looks in parameter space. .. _sec-informative-priors: Informative Priors ~~~~~~~~~~~~~~~~~~~ When genuine prior knowledge exists, encoding it in the prior is not only permissible but statistically optimal. The choice is not between "Bayesian (subjective) vs. frequentist (objective)" β€” it is between using available information and discarding it. **Sources of prior information:** - *Previous studies:* If three earlier trials of a related drug found response rates of 0.62, 0.68, and 0.71, the prior for a new trial should reflect this evidence. A Beta prior fitted to these historical rates is a *meta-analytic prior* β€” arguably the most defensible form of informative prior because its basis is explicit and auditable. - *Physical constraints:* A decay constant must be positive; a probability must lie in :math:`[0,1]`; a measurement precision must be positive. These are not opinions but mathematical facts, and the prior should reflect them. - *Expert elicitation:* Domain experts can be asked for predictive judgments β€” "What is the median response rate? What is the 90th percentile?" β€” from which hyperparameters can be solved algebraically. This *percentile matching* approach translates expert knowledge into prior hyperparameters without requiring the expert to reason about probability distributions directly. **Percentile matching for Beta priors.** Suppose an expert states that the response rate is "probably around 0.65, and almost certainly below 0.85." Encoding these as a median of 0.65 and a 90th percentile of 0.85: .. code-block:: python import numpy as np from scipy import optimize, stats def match_beta_to_quantiles(q1: float, p1: float, q2: float, p2: float) -> tuple[float, float]: """ Find Beta(alpha, beta) hyperparameters matching two quantiles. Parameters ---------- q1, q2 : float Quantile probabilities (e.g., 0.50 for median, 0.90 for 90th pctile). p1, p2 : float Corresponding parameter values (the quantile values). Returns ------- tuple[float, float] (alpha, beta) hyperparameters. """ def objective(params: np.ndarray) -> np.ndarray: alpha, beta = np.exp(params) # ensure positivity dist = stats.beta(alpha, beta) return np.array([dist.ppf(q1) - p1, dist.ppf(q2) - p2]) # Initial guess: method of moments from (p1+p2)/2 and 0.1 variance result = optimize.fsolve(objective, x0=[1.0, 1.0], full_output=True) log_ab = result[0] alpha, beta = np.exp(log_ab) return alpha, beta alpha_elicited, beta_elicited = match_beta_to_quantiles( q1=0.50, p1=0.65, q2=0.90, p2=0.85 ) elicited = stats.beta(alpha_elicited, beta_elicited) print(f"Elicited hyperparameters: alpha={alpha_elicited:.2f}, " f"beta={beta_elicited:.2f}") print(f"Verification - median: {elicited.ppf(0.50):.3f} " f"(target: 0.650)") print(f"Verification - 90th pctile: {elicited.ppf(0.90):.3f} " f"(target: 0.850)") print(f"Prior effective sample size: {alpha_elicited+beta_elicited:.1f}") .. code-block:: text Elicited hyperparameters: alpha=3.54, beta=2.04 Verification - median: 0.650 (target: 0.650) Verification - 90th pctile: 0.850 (target: 0.850) Prior effective sample size: 5.6 **Documentation as discipline.** The requirement to specify a prior explicitly imposes a discipline that is itself valuable: it forces the analyst to state their assumptions clearly, where they can be examined, challenged, and revised. A Bayesian analysis with explicit priors is, in this sense, more transparent than a frequentist analysis that makes equally strong implicit assumptions through model specification while hiding them. .. _sec-sensitivity-analysis: Prior Sensitivity Analysis ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Even after careful prior specification, it is good practice to examine how the posterior changes under alternative prior choices. This *sensitivity analysis* quantifies the degree to which the conclusion depends on the prior, which is especially important in high-stakes applications such as clinical trials, policy analysis, or legal proceedings. **Algorithm: Prior Sensitivity Cascade** .. code-block:: text Input: Data y; family of priors {p_s(theta)} indexed by s; model class Output: Table of posterior summaries for each prior 1. For each prior specification s: a. Compute posterior p_s(theta | y) (analytically or via MCMC) b. Summarize: posterior mean, median, 95% CrI 2. Examine: How much do the summaries vary across priors? 3. Check: Does any prior produce a posterior substantially different from the others? 4. Report: All priors and their posterior summaries; flag cases of high sensitivity. **Computational Complexity**: :math:`O(M \times C)` where :math:`M` is the number of prior specifications and :math:`C` is the cost of a single posterior computation. .. code-block:: python import numpy as np import pandas as pd from scipy import stats def prior_sensitivity_beta_binomial( k: int, n: int, prior_specs: list[tuple[float, float]], level: float = 0.95 ) -> pd.DataFrame: """ Compare posterior summaries across multiple Beta prior specifications for a fixed Beta-Binomial dataset. Parameters ---------- k : int Observed successes. n : int Total trials. prior_specs : list of (alpha0, beta0) tuples Prior hyperparameter specifications. level : float Credible interval coverage. Returns ------- pd.DataFrame Summary table with one row per prior specification. """ rows = [] mle = k / n for (a0, b0) in prior_specs: model = BetaBinomialModel(a0, b0).update(k, n) lo, hi = model.credible_interval(level) rows.append({ 'Prior': f'Beta({a0:.1f}, {b0:.1f})', 'n0 (eff. prior n)': a0 + b0, 'w_data': f"{model.prior_data_weight:.3f}", 'Post Mean': f"{model.posterior_mean:.4f}", 'Post Std': f"{np.sqrt(model.posterior_variance):.4f}", f'{int(level*100)}% CrI Lower': f"{lo:.4f}", f'{int(level*100)}% CrI Upper': f"{hi:.4f}", }) df = pd.DataFrame(rows) df.loc[len(df)] = {'Prior': f'MLE (n={n})', 'n0 (eff. prior n)': 0, 'w_data': '1.000', 'Post Mean': f'{mle:.4f}', 'Post Std': 'β€”', '95% CrI Lower': 'β€”', '95% CrI Upper': 'β€”'} return df # --- Sensitivity analysis: drug trial (k=15, n=20) --- prior_specs_small = [ (1.0, 1.0), # flat (0.5, 0.5), # Jeffreys (2.0, 2.0), # mild (5.0, 5.0), # moderate (10.0, 10.0), # skeptical (30.0, 30.0), # strongly skeptical ] df_small = prior_sensitivity_beta_binomial(15, 20, prior_specs_small) print("Drug trial (k=15, n=20) β€” Prior Sensitivity:") print(df_small.to_string(index=False)) # --- Same proportions, large sample --- df_large = prior_sensitivity_beta_binomial(150, 200, prior_specs_small) print("\nLarge-sample analog (k=150, n=200) β€” Prior Sensitivity:") print(df_large.to_string(index=False)) .. code-block:: text Drug trial (k=15, n=20) β€” Prior Sensitivity: Prior n0 (eff. prior n) w_data Post Mean Post Std 95% CrI Lower 95% CrI Upper Beta(1.0, 1.0) 2 0.909 0.7273 0.0929 0.5283 0.8872 Beta(0.5, 0.5) 1 0.952 0.7381 0.0937 0.5358 0.8976 Beta(2.0, 2.0) 4 0.833 0.7083 0.0909 0.5159 0.8679 Beta(5.0, 5.0) 10 0.667 0.6667 0.0847 0.4917 0.8206 Beta(10.0, 10.0) 20 0.500 0.6250 0.0756 0.4718 0.7664 Beta(30.0, 30.0) 60 0.250 0.5625 0.0551 0.4533 0.6688 MLE (n=20) 0 1.000 0.7500 β€” β€” β€” Large-sample analog (k=150, n=200) β€” Prior Sensitivity: Prior n0 (eff. prior n) w_data Post Mean Post Std 95% CrI Lower 95% CrI Upper Beta(1.0, 1.0) 2 0.990 0.7475 0.0305 0.6855 0.8049 Beta(0.5, 0.5) 1 0.995 0.7488 0.0305 0.6867 0.8061 Beta(2.0, 2.0) 4 0.980 0.7451 0.0304 0.6833 0.8024 Beta(5.0, 5.0) 10 0.952 0.7381 0.0303 0.6767 0.7952 Beta(10.0, 10.0) 20 0.909 0.7273 0.0300 0.6667 0.7840 Beta(30.0, 30.0) 60 0.769 0.6923 0.0286 0.6350 0.7468 MLE (n=200) 0 1.000 0.7500 β€” β€” β€” **Interpreting sensitivity results.** With :math:`n = 20`, even the strongly skeptical Beta(30, 30) prior (contributing 60 pseudo-observations) pulls the posterior mean from 0.750 to 0.563, a substantial change. With :math:`n = 200`, the same prior moves the posterior mean only from 0.750 to 0.692 β€” a difference of less than six percentage points that would rarely change a practical conclusion. This is the Bernstein–von Mises effect in action: with enough data, even a strongly misspecified prior is overwhelmed. If a sensitivity analysis reveals that conclusions differ substantially across a plausible range of priors, the analyst should report all results and discuss which prior best represents available knowledge. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig14_sensitivity_cascade.png :alt: Two-column figure showing prior sensitivity cascades for the drug trial. Left column (n=20): five posterior distributions for each prior specification plotted together, illustrating substantial divergence between flat and strongly skeptical priors. Right column (n=200): the same five posteriors are tightly clustered near the MLE, with prior influence nearly invisible. An inset scatter plot shows how posterior mean varies across prior specifications for both sample sizes. :align: center :width: 100% **Figure 5.14**: Prior sensitivity cascade for the drug trial dataset. Left: with :math:`n = 20` observations, priors ranging from flat (Beta(1,1)) to strongly skeptical (Beta(30,30)) produce posterior means spanning the range 0.563–0.738 β€” a clinically meaningful difference. Right: with :math:`n = 200` observations at the same proportion, all six priors yield posteriors concentrated near the MLE (0.750), with means ranging from 0.692 to 0.749. The prior is consequential when :math:`n` is small; it washes out as :math:`n` grows. Bayesian vs. Frequentist Synthesis and the Limits of Conjugacy --------------------------------------------------------------- We have now developed complete analytical machinery for five conjugate families and a principled framework for choosing priors. Before moving forward, it is useful to consolidate what conjugate analysis can and cannot do, and to precisely identify where the remaining chapter sections become necessary. Comparing Bayesian and Frequentist Answers ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The following table summarizes the Bayesian-frequentist comparison across all five conjugate families developed in this section. In each case, the comparison is under the reference (Jeffreys or flat) prior. .. list-table:: Bayesian vs. Frequentist Comparison for Conjugate Families :header-rows: 1 :widths: 18 20 20 22 20 * - Model - Bayes estimate (ref. prior) - Frequentist MLE - Asymptotic agreement? - When they diverge * - Beta-Binomial - :math:`\frac{k+1/2}{n+1}` (Jeffreys) - :math:`k/n` - Yes (:math:`n \to \infty`) - Small :math:`n`; boundary cases :math:`k=0` or :math:`k=n` * - Normal-Normal (known :math:`\sigma^2`) - :math:`\bar{y}` (flat prior) - :math:`\bar{y}` - Exact for all :math:`n` - Never (Bayes CrI = freq. CI under flat prior) * - Normal-Inverse-Gamma - :math:`\bar{y}` (ref. prior); :math:`t_n` marginal - :math:`\bar{y}`; :math:`t_{n-1}` pivotal - Near-exact; gap :math:`O(1/n)` (df :math:`= n` vs :math:`n-1`) - Negligible for :math:`n \geq 30`; diverge under informative prior * - Poisson-Gamma - :math:`T/n` approximately (Jeffreys :math:`\alpha_0=0.5`) - :math:`T/n` - Yes (:math:`n \to \infty`) - Small :math:`n`; zero-count data * - Multinomial-Dirichlet - :math:`(n_k + 0.5)/(n + K/2)` (Jeffreys) - :math:`n_k/n` - Yes (:math:`n \to \infty`) - Small :math:`n`; empty cells The pattern is consistent: under reference priors and with sufficient data, Bayesian and frequentist answers agree. With informative priors or sparse data, the Bayesian answer incorporates the prior information and will generally differ β€” this is precisely the setting where the Bayesian framework adds value. The credible interval and the confidence interval, while numerically similar under reference priors, retain their different interpretations: the credible interval is a probability statement about :math:`\theta` given the observed data; the confidence interval is a statement about the coverage probability of the procedure in repeated sampling. The Limits of Conjugate Analysis ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Conjugate analysis requires the prior to belong to a specific parametric family matched to the likelihood. This is not merely a preference β€” it is an algebraic requirement. When the model contains non-conjugate elements, closed-form posteriors are unavailable. The three most common situations: **Logistic regression.** With outcome :math:`y_i \in \{0, 1\}` and predictors :math:`x_i`, the logistic model is: .. math:: P(y_i = 1 \mid \beta) = \text{logit}^{-1}(x_i^\top \beta) = \frac{\exp(x_i^\top \beta)}{1 + \exp(x_i^\top \beta)} Place a Normal prior :math:`\beta \sim \mathcal{N}(0, \sigma_0^2 I)` on the coefficient vector. The posterior is: .. math:: p(\beta \mid y) \propto \prod_{i=1}^n \left[\text{logit}^{-1}(x_i^\top \beta)\right]^{y_i} \left[1 - \text{logit}^{-1}(x_i^\top \beta)\right]^{1-y_i} \cdot \exp\!\left(-\frac{\|\beta\|^2}{2\sigma_0^2}\right) The logistic function :math:`\text{logit}^{-1}(\cdot)` is not a polynomial or exponential in :math:`\beta`, so the product of the likelihood and the Normal prior has no recognizable kernel. The normalizing constant :math:`\int p(y \mid \beta) p(\beta) \, d\beta` is a multidimensional integral with no closed form. Grid approximation could handle this for one or two predictors, but logistics regression in practice involves tens or hundreds of predictors β€” the curse of dimensionality makes grid approximation infeasible. **Mixture models.** A Gaussian mixture model with :math:`K` components has :math:`3K - 1` parameters (means, variances, mixture weights). Even with conjugate component priors (Normal-Inverse-Gamma for each component) and a Dirichlet prior for the weights, the likelihood involves sums over component assignments β€” log-sums that prevent kernel recognition. **Hierarchical models.** Models with multiple levels of parameter uncertainty (group means drawn from a population distribution, which has its own hyperprior) can be conjugate at each level but not jointly, requiring full joint posterior computation. The gap between "analytically tractable" and "everything else" is the chasm that MCMC bridges. Rather than computing the posterior analytically, MCMC constructs a Markov chain that samples from the posterior without evaluating the normalizing constant β€” making Bayesian inference feasible for essentially any well-specified probabilistic model. We develop the MCMC machinery in :ref:`ch5_4-markov-chains` and :ref:`ch5_5-mcmc-algorithms`, returning to these non-conjugate examples throughout. .. figure:: https://pqyjaywwccbnqpwgeiuv.supabase.co/storage/v1/object/public/STAT%20418%20Images/assets/PartIII/Chapter5/ch5_2_fig16_beyond_conjugacy.png :alt: Conceptual diagram showing three columns: conjugate models (closed-form posteriors, low-dimensional, this section), grid approximation zone (exact but limited to 1-2 parameters, Section 5.1), and MCMC territory (general-purpose, high-dimensional, Sections 5.4-5.5); arrows indicate that logistic regression and hierarchical models fall in the MCMC zone :align: center :width: 100% **Figure 5.16**: The Bayesian computation landscape. Conjugate analysis (left) provides closed-form posteriors for matched prior-likelihood pairs β€” tractable, exact, but limited to specific parametric families. Grid approximation (center) is exact in principle but infeasible beyond two or three parameters. MCMC (right) provides approximate posterior samples for virtually any well-specified model at the cost of convergence diagnostics and longer computation time. The arrows show where logistic regression, mixture models, and hierarchical models fall: entirely in the MCMC zone. Practical Considerations ------------------------ **Hyperparameter specification in practice.** When choosing a conjugate prior, begin with the pseudo-data interpretation: how many observations' worth of prior information do you have? If a pilot study provides 10 data points, encode that as an effective prior sample size of :math:`n_0 = 10`. If you have only a rough sense of the mean and no strong prior precision, choose :math:`n_0 = 1`–2. Start with a weakly informative prior and run a prior predictive simulation to check that it generates plausible data. **Numerical stability.** For the Beta-Binomial marginal likelihood computation, always use :code:`scipy.special.betaln` (log Beta function) rather than computing :math:`B(\alpha, \beta)` directly β€” the Beta function can underflow to zero for large hyperparameters even when the log-Beta is finite. Similarly, for the NIG model, accumulate the posterior sum of squares :math:`\nu_n \sigma_n^2` rather than :math:`\sigma_n^2` separately, since dividing the accumulated quantity by :math:`\nu_n` introduces a single division rather than multiple multiplications and divisions that can magnify rounding errors. **The sequential updating pattern.** All conjugate families admit sequential updating: today's posterior is tomorrow's prior. This is computationally free for conjugate models β€” just pass the posterior hyperparameters as the new prior. In practice, sequential updating can be used to track how the posterior evolves as data accumulate, which provides a natural internal consistency check: the final posterior should be identical regardless of whether data are processed in one batch or sequentially. **When to abandon conjugacy.** Conjugate analysis is a starting point, not an endpoint. If your scientific question requires a model that does not fit any standard conjugate family β€” a Poisson likelihood with covariates, a Normal model with outliers requiring a t-likelihood, a Beta-Binomial with informative covariates β€” you should use the more flexible model rather than forcing the data into a conjugate mold. The methods in Sections 5.4–5.6 make the more general Bayesian analysis computationally feasible. .. admonition:: Common Pitfall ⚠️ :class: warning Improper posteriors can arise silently when improper priors are combined with sparse data. The Jeffreys prior for the Poisson rate, :math:`p(\lambda) \propto \lambda^{-1/2}`, yields a proper posterior for :math:`n \geq 1`, but only because :math:`\int_0^\infty \lambda^{-1/2} \cdot \lambda^T e^{-n\lambda} d\lambda < \infty` when :math:`T + 1/2 > 0`. If you use an improper prior that is *not* a standard Jeffreys prior β€” for example, a flat prior on a variance parameter β€” verify analytically that the posterior is proper before running any computation. Improper posteriors can produce samplers that appear to converge while actually exploring an unnormalized target, yielding nonsensical parameter estimates that show no obvious diagnostic sign of failure. Bringing It All Together ------------------------ This section has developed the Bayesian machinery for the most important class of tractable models β€” conjugate families β€” while simultaneously establishing the broader framework of principled prior specification. The two threads are inseparable: conjugate priors are not merely convenient computational devices but the most carefully studied instances of prior specification, and understanding them quantitatively is the best preparation for choosing priors in the non-conjugate setting. Every conjugate family follows the same pattern: the prior is recognized as the natural exponential family member matched to the likelihood, the posterior update is addition of sufficient statistics to hyperparameters, and the posterior mean is a precision-weighted average of the prior mean and the MLE. The weights quantify how much the prior matters as a function of the sample size and the prior's effective sample size. This pattern β€” and the prior predictive simulation framework for evaluating priors before seeing data β€” gives practitioners concrete tools for making prior choices that are transparent, defensible, and calibrated. The section also clarified exactly where Bayesian computation runs into a wall: the logistic regression posterior, the Gaussian mixture posterior, and hierarchical model posteriors are analytically intractable because the likelihood Γ— prior product has no recognizable kernel. These are not pathological edge cases β€” they are the models that dominate modern applied Bayesian data analysis. The next section (:ref:`ch5_3-credible-intervals`) develops the interval estimation and predictive interval tools that flow from conjugate and approximate posteriors. Sections 5.4–5.5 then develop MCMC to unlock the non-conjugate models where closed-form analysis ends. .. admonition:: Key Takeaways πŸ“ :class: important 1. **Core concept:** Conjugate priors yield closed-form posteriors because the prior's kernel matches the likelihood's sufficient statistic structure. The posterior hyperparameters are obtained by adding observed sufficient statistics to prior hyperparameters β€” a "pseudo-data" interpretation that gives every hyperparameter a concrete meaning in units of observations. [LO 2, 4] 2. **Computational insight:** The posterior mean under any conjugate model is a precision-weighted average of the prior mean and the MLE, with data weight :math:`w = n/(n + n_0)`. This single formula quantifies the prior-data balance: increasing :math:`n` drives :math:`w \to 1`, and the prior washes out asymptotically β€” a finite-sample demonstration of the Bernstein–von Mises theorem. [LO 2, 4] 3. **Invariance and reference priors:** The Jeffreys prior :math:`p_J(\theta) \propto \sqrt{I_1(\theta)}` is the principled parameterization-invariant choice for "uninformative" analysis. For the Binomial it yields Beta(1/2, 1/2); for the Normal scale it yields :math:`1/\sigma`; for the Poisson rate it yields an improper :math:`\lambda^{-1/2}`. Under reference priors, Bayesian credible intervals closely match frequentist confidence intervals in the Normal-Inverse-Gamma model (differing by one degree of freedom, negligible for :math:`n \geq 30`), and agree asymptotically for all conjugate families. [LO 2, 4] 4. **Practical application:** Prior predictive simulation β€” draw parameters from the prior, simulate data from the likelihood, check whether simulated data look plausible β€” is the most effective tool for evaluating prior choices. It reveals that apparently "vague" priors (e.g., Normal(0, 100Β²)) often imply physically implausible data, while weakly informative priors (Normal(0, 1Β²) on standardized predictors) gently regularize without dominating the likelihood. [LO 4] 5. **Limitation and transition:** Conjugate analysis requires the prior to match the likelihood's exponential family structure. Logistic regression, mixture models, hierarchical models, and most real-world applications are non-conjugate: the posterior integral is analytically intractable. Grid approximation handles one or two parameters but scales exponentially in dimensionality. MCMC (Sections 5.4–5.5) resolves this by sampling from the posterior without computing the normalizing constant. [LO 4] Exercises --------- .. admonition:: Exercise 5.2.1: Gamma-Exponential Conjugate Derivation :class: exercise The Exponential distribution :math:`y \sim \text{Exp}(\lambda)` with rate parameter :math:`\lambda > 0` has likelihood proportional to :math:`\lambda^n e^{-\lambda \sum_i y_i}`. Derive the conjugate prior for :math:`\lambda` and identify it as a member of a named parametric family. **(a)** Starting from the Exponential likelihood for :math:`n` i.i.d. observations :math:`y_1, \ldots, y_n`, identify the sufficient statistic :math:`T = T(y_1, \ldots, y_n)`. Write the likelihood as a function of :math:`\lambda` and :math:`T`. **(b)** Propose a prior :math:`p(\lambda) \propto \lambda^{a-1} e^{-b\lambda}` and multiply it by the likelihood. Show that the posterior has the same kernel as the prior. Identify the posterior parameters in terms of the prior hyperparameters and :math:`T`. **(c)** Your posterior is :math:`\text{Gamma}(\alpha_0 + n, \beta_0 + T)`. Verify this matches the Poisson-Gamma update pattern: what is the pseudo-data interpretation of :math:`\alpha_0` and :math:`\beta_0` for the Exponential model? **(d)** **Python verification.** Generate :math:`n = 20` observations from :math:`\text{Exp}(2.0)` using :code:`rng.exponential(scale=0.5, size=20)`. Use a Gamma(2, 1) prior (prior mean :math:`\lambda_0 = 2`). Compute: - The analytical posterior using your derived update rule. - A grid-based posterior approximation on :math:`\lambda \in (0, 10)` with 1,000 grid points. - The posterior mean, mode, and 95% credible interval for both methods. - The maximum absolute difference between the normalized grid posterior and the analytical Gamma pdf on the grid. **Expected result:** The grid and analytical posteriors should agree to at most :math:`10^{-4}` in absolute value. .. dropdown:: Solution :class-container: solution **(a)** The sufficient statistic is :math:`T = \sum_{i=1}^n y_i` (the sum of observations, equivalently the sample mean :math:`\bar{y} \cdot n`). The likelihood is :math:`p(y \mid \lambda) \propto \lambda^n \exp(-\lambda T)`. **(b)** Multiplying: :math:`p(\lambda \mid y) \propto \lambda^n e^{-\lambda T} \cdot \lambda^{\alpha_0 - 1} e^{-\beta_0 \lambda} = \lambda^{(\alpha_0 + n) - 1} e^{-(\beta_0 + T)\lambda}`. This is Gamma(:math:`\alpha_0 + n, \beta_0 + T`). **(c)** :math:`\alpha_0` is the prior count of "pseudo-observations"; :math:`\beta_0` is the total prior "lifetime" (sum of pseudo-observations' times). The prior mean of the rate is :math:`\alpha_0/\beta_0`, which equals the ratio of pseudo-counts to pseudo-total-time β€” exactly the empirical rate formula applied to the pseudo-data. **(d)** Typical output: .. code-block:: python import numpy as np from scipy import stats from scipy.special import logsumexp rng = np.random.default_rng(42) lam_true = 2.0 n = 20 y = rng.exponential(1.0 / lam_true, size=n) T = y.sum() alpha0, beta0 = 2.0, 1.0 alpha_n = alpha0 + n # = 22 beta_n = beta0 + T # = beta0 + sum(y) # Analytical posterior post_dist = stats.gamma(a=alpha_n, scale=1.0/beta_n) post_mean = post_dist.mean() post_mode = (alpha_n - 1) / beta_n lo, hi = post_dist.ppf(0.025), post_dist.ppf(0.975) print(f"Analytical posterior Gamma({alpha_n:.1f}, rate={beta_n:.3f})") print(f" Mean: {post_mean:.4f}, Mode: {post_mode:.4f}") print(f" 95% CrI: ({lo:.4f}, {hi:.4f})") # Grid verification lam_grid = np.linspace(1e-3, 10, 1000) log_lik = n * np.log(lam_grid) - lam_grid * T log_prior = (alpha0 - 1) * np.log(lam_grid) - beta0 * lam_grid log_unnorm = log_lik + log_prior log_post = log_unnorm - logsumexp(log_unnorm) post_grid = np.exp(log_post) post_grid /= post_grid.sum() # Compare with analytical Gamma pdf analytical_pdf = post_dist.pdf(lam_grid) analytical_pdf /= analytical_pdf.sum() max_err = np.max(np.abs(post_grid - analytical_pdf)) print(f"\nGrid vs. Analytical max error: {max_err:.2e}") .. admonition:: Exercise 5.2.2: Effective Prior Sample Size Across Conjugate Families :class: exercise For each of the five conjugate families in this section, the posterior mean can be written as a weighted average of the prior mean and the MLE: .. math:: \mathbb{E}[\theta \mid y] = (1 - w)\mu_0 + w \hat{\theta}_{\text{MLE}}, \qquad w = \frac{n}{n + n_0} where :math:`n_0` is the effective prior sample size. **(a)** Complete the following table by deriving the effective prior sample size :math:`n_0` for each family: .. list-table:: :header-rows: 1 :widths: 20 25 20 35 * - Family - Prior - MLE - :math:`n_0` (your answer) * - Beta-Binomial - Beta(:math:`\alpha_0, \beta_0`) - :math:`k/n` - ? * - Normal-Normal (known :math:`\sigma^2`) - :math:`\mathcal{N}(\mu_0, \sigma_0^2)` - :math:`\bar{y}` - ? * - Poisson-Gamma - Gamma(:math:`\alpha_0, \beta_0`) - :math:`T/n` - ? * - Multinomial-Dirichlet - Dir(:math:`\alpha_1, \ldots, \alpha_K`) - :math:`n_k/n` - ? **(b)** For each family, at what sample size :math:`n^*` does the data weight :math:`w` reach 0.95 (data contributes 95% of the posterior)? Express :math:`n^*` in terms of :math:`n_0`. **(c)** **Python implementation.** Write a function :code:`data_weight(n, n0)` that plots :math:`w(n)` as a function of :math:`n` for a user-supplied :math:`n_0`. Plot five curves for :math:`n_0 \in \{1, 5, 10, 50, 100\}` on the same axes, and mark the :math:`n^*` crossings. .. dropdown:: Solution :class-container: solution **(a)** Effective prior sample sizes: - Beta-Binomial: :math:`n_0 = \alpha_0 + \beta_0` - Normal-Normal: :math:`n_0 = \sigma^2 / \sigma_0^2` (ratio of sampling variance to prior variance) - Poisson-Gamma: :math:`n_0 = \beta_0` (prior exposure) - Multinomial-Dirichlet: :math:`n_0 = \sum_k \alpha_k` **(b)** Setting :math:`w = 0.95`: :math:`n/(n + n_0) = 0.95` gives :math:`n^* = 19 n_0`. To reach 95% data weight, you need 19 times the prior's effective sample size. .. admonition:: Exercise 5.2.3: Jeffreys Prior for the Exponential Distribution :class: exercise The Exponential(:math:`\lambda`) distribution has log-likelihood :math:`\log p(y \mid \lambda) = \log \lambda - \lambda y` for a single observation. **(a)** Compute the Fisher information :math:`I_1(\lambda)` using the score variance formula :math:`I_1(\lambda) = \mathbb{E}[(\partial \log p/\partial \lambda)^2]`. Show that :math:`I_1(\lambda) = 1/\lambda^2`. **(b)** Derive the Jeffreys prior :math:`p_J(\lambda) \propto \sqrt{I_1(\lambda)}`. What named distribution does this correspond to? **(c)** Show that the Jeffreys prior is improper (does not integrate to a finite value over :math:`\lambda > 0`). **(d)** Apply the conjugate update rule from Exercise 5.2.1 with the Jeffreys prior as the "prior" (formally: :math:`\alpha_0 = 1/2, \beta_0 = 0`). Show that the resulting posterior is proper for :math:`n \geq 1`. **(e)** **Invariance demonstration.** Reparameterize by :math:`\phi = 1/\lambda` (the mean of the Exponential distribution). Derive the Jeffreys prior for :math:`\phi` directly by computing the Fisher information in the :math:`\phi` parameterization. Verify that applying the Jacobian transformation :math:`p_J(\phi) = p_J(\lambda(\phi))|d\lambda/d\phi|` to the :math:`\lambda` Jeffreys prior gives the same result. .. dropdown:: Solution :class-container: solution **(a)** Score: :math:`\partial \log p / \partial \lambda = 1/\lambda - y`. Variance: :math:`\text{Var}(1/\lambda - Y) = \text{Var}(Y) = 1/\lambda^2` (since :math:`\text{Var}(Y) = 1/\lambda^2` for :math:`\text{Exp}(\lambda)`). Therefore :math:`I_1(\lambda) = 1/\lambda^2`. **(b)** :math:`p_J(\lambda) \propto \lambda^{-1}`, the scale-invariant improper prior. This is formally :math:`\text{Gamma}(0, 0)` β€” improper. **(c)** :math:`\int_0^\infty 1/\lambda \, d\lambda = [\log\lambda]_0^\infty = \infty`. The integral diverges at both 0 and :math:`\infty`. **(d)** Posterior: :math:`\text{Gamma}(1/2 + n, T)`. For :math:`n \geq 1`, the shape parameter is :math:`\geq 1/2 > 0` and scale parameter :math:`T > 0` a.s., so the Gamma distribution is proper. **(e)** For :math:`\phi = 1/\lambda`, the Exponential mean parameterization gives :math:`p(y|\phi) = (1/\phi)e^{-y/\phi}`. Score: :math:`\partial \log p/\partial\phi = -1/\phi + y/\phi^2`. Fisher: :math:`I_1^{(\phi)}(\phi) = \text{Var}(Y/\phi^2 - 1/\phi) = \text{Var}(Y)/\phi^4 = \phi^2/\phi^4 = 1/\phi^2`. Jeffreys: :math:`p_J(\phi) \propto 1/\phi`. Jacobian check: :math:`\lambda = 1/\phi`, :math:`|d\lambda/d\phi| = 1/\phi^2`, so :math:`p_J(\phi) \propto \sqrt{I_1^{(\lambda)}(1/\phi)} \cdot (1/\phi^2) = (1/\phi)^{-1} \cdot (1/\phi^2) = 1/\phi`. The two derivations agree β€” invariance holds. .. admonition:: Exercise 5.2.4: Prior Predictive Simulation for Logistic Regression :class: exercise Consider a logistic regression with one standardized predictor :math:`x \in [-2, 2]`: .. math:: P(y_i = 1 \mid \alpha, \beta, x_i) = \frac{1}{1 + e^{-(\alpha + \beta x_i)}} **(a)** For :math:`\alpha \sim \mathcal{N}(0, 1.5^2)` and three candidate priors for :math:`\beta`: - :math:`\beta \sim \mathcal{N}(0, 10^2)` (vague) - :math:`\beta \sim \mathcal{N}(0, 1^2)` (weakly informative) - :math:`\beta \sim \mathcal{N}(0, 0.5^2)` (informative) Generate :math:`S = 500` prior predictive probability curves :math:`P(y=1 \mid x)` for each prior. What does each prior imply about the steepness and direction of the relationship between :math:`x` and the probability of success? **(b)** For the vague prior, what fraction of simulated probability curves have :math:`|P(y=1|x=2) - P(y=1|x=-2)| > 0.9` (nearly binary flip across the predictor range)? **(c)** Which prior do you recommend, and why? Justify in terms of what the prior implies about the relative change in probability per unit change in :math:`x`. **(d)** **McElreath comparison.** For each prior, compute the prior predictive mean of :math:`|P(y=1|x=2) - P(y=1|x=-2)|` β€” the total variation in the predicted probability across the predictor range. Organize the results in a table. .. dropdown:: Solution :class-container: solution **(a)** .. code-block:: python import numpy as np from scipy.special import expit # logistic function rng = np.random.default_rng(42) x = np.linspace(-2, 2, 200) S = 500 for label, sigma_beta in [('Vague', 10.0), ('Weakly informative', 1.0), ('Informative', 0.5)]: alpha = rng.normal(0, 1.5, size=S) beta = rng.normal(0, sigma_beta, size=S) # Probability curves: shape (S, 200) eta = alpha[:, None] + beta[:, None] * x[None, :] p_curves = expit(eta) total_var = np.abs(p_curves[:, -1] - p_curves[:, 0]) print(f"{label} (sigma_beta={sigma_beta}):") print(f" Mean |P(x=2)-P(x=-2)|: {total_var.mean():.3f}") print(f" Fraction > 0.9: {(total_var > 0.9).mean():.3f}") Typical output: Vague (sigma=10) mean total variation ~0.85, fraction > 0.9 ~55%; Weakly informative (sigma=1) mean ~0.57, fraction > 0.9 ~10%; Informative (sigma=0.5) mean ~0.37, fraction > 0.9 ~1%. **(c)** The weakly informative Normal(0, 1) prior is appropriate: it allows curves that span most of the probability range when appropriate (:math:`\beta \approx \pm 2` gives total variation ~0.7) without concentrating mass on near-vertical decision boundaries. The vague prior is problematic β€” it generates near-binary step functions (>55% of curves with total variation >0.9), implying almost perfect discrimination in nearly all simulations, which is implausible before seeing data. The informative prior (sigma=0.5) may be too conservative if the predictor is expected to have a moderate-to-large effect. .. admonition:: Exercise 5.2.5: Prior Sensitivity and the Bernstein-von Mises Effect :class: exercise This exercise quantifies the Bernstein-von Mises convergence numerically for the Beta-Binomial model. Fix the observed proportion :math:`k/n = 0.76` (as in the Section 5.1 vaccine trial) and vary :math:`n` over :math:`\{5, 20, 50, 100, 500\}`. For each :math:`n`, compute :math:`k = \text{round}(0.76 \cdot n)`. **(a)** For each :math:`n`, compute the posterior under five priors: Beta(1,1), Beta(2,2), Beta(5,5), Beta(10,10), Beta(30,30). Record the posterior mean, posterior standard deviation, and 95% credible interval lower and upper bounds. Organize results in a table. **(b)** Define the *posterior spread* as the range of posterior means across the five priors (max minus min). Plot the posterior spread as a function of :math:`n`. At what :math:`n` does the spread fall below 0.01 (1 percentage point)? **(c)** Show analytically that the posterior spread converges to 0 as :math:`n \to \infty` for any fixed set of priors with finite effective sample sizes. (Hint: use the weighted average formula.) **(d)** Compare the width of each 95% credible interval (upper - lower) with the corresponding frequentist Wald interval width :math:`2 \times 1.96 \times \sqrt{k/n (1-k/n)/n}`. At what :math:`n` do all five Bayesian credible intervals agree with the Wald interval to within 0.01? .. dropdown:: Solution :class-container: solution **(c)** The posterior mean for prior :math:`i` with effective size :math:`n_0^{(i)}` is :math:`\mu_n^{(i)} = (n_0^{(i)} \mu_0^{(i)} + n\hat{\theta})/(n_0^{(i)} + n)`. The spread between two priors :math:`i` and :math:`j` is: .. math:: |\mu_n^{(i)} - \mu_n^{(j)}| = \left|\frac{n_0^{(i)}(\mu_0^{(i)} - \hat{\theta})}{n_0^{(i)} + n} - \frac{n_0^{(j)}(\mu_0^{(j)} - \hat{\theta})}{n_0^{(j)} + n}\right| \leq \frac{n_0^{(i)}}{n} + \frac{n_0^{(j)}}{n} \to 0 since :math:`n_0^{(i)}` and :math:`n_0^{(j)}` are fixed and :math:`n \to \infty`. .. admonition:: Exercise 5.2.6: Normal-Inverse-Gamma β€” Conjugate vs. Grid vs. Frequentist :class: exercise **(a)** Generate :math:`n = 25` observations from :math:`\mathcal{N}(7.0, 3.0^2)` using :code:`rng = np.random.default_rng(2024)`. Compute: - The analytical NIG posterior using a reference prior (:math:`\kappa_0 = 10^{-6}`, :math:`\nu_0 = 10^{-6}`, :math:`\mu_0 = 0`, :math:`\sigma_0^2 = 1`). - The grid posterior using :code:`grid_posterior_2d` from Section 5.1 on a 300Γ—300 grid. - The marginal posterior for :math:`\mu` from both methods. - The 95% credible interval for :math:`\mu` from both methods. - The frequentist :math:`t`-interval. **(b)** Verify numerically that the NIG marginal posterior for :math:`\mu` (equation :eq:`nig-marginal-mu`) matches a :math:`t_{n-1}` distribution centered at :math:`\bar{y}` under the reference prior. Compute the maximum absolute difference between the analytical :math:`t`-pdf and the NIG marginal :math:`t`-pdf evaluated on a fine grid. **(c)** Now use an **informative** NIG prior: :math:`\mu_0 = 5`, :math:`\kappa_0 = 5`, :math:`\nu_0 = 5`, :math:`\sigma_0^2 = 4`. Recompute the posterior for :math:`\mu`. Compare the informative-prior credible interval with the reference-prior interval. By how much does the posterior mean shift? Is the shift in the expected direction? **(d)** At what sample size :math:`n^*` does the informative-prior NIG posterior become indistinguishable from the reference-prior NIG posterior to within 0.01 in posterior mean? (Increase :math:`n` by resampling from the same Normal distribution.) .. dropdown:: Solution :class-container: solution **(a)** With :code:`rng.normal(7.0, 3.0, 25)`, typical results (:math:`\bar{y} \approx 7.0`, :math:`s \approx 2.9`): NIG posterior for :math:`\mu` is :math:`t_{25}(\bar{y}, s^2/25)`. All three 95% intervals should agree to 3+ decimal places under the reference prior. **(c)** With :math:`\mu_0 = 5 < 7 \approx \bar{y}`, the informative prior pulls the posterior mean toward 5. With :math:`\kappa_0 = 5` pseudo-observations vs. :math:`n = 25` real observations, the data weight is :math:`w = 25/30 \approx 0.833`, so the posterior mean shifts :math:`0.833 \times (\bar{y} - 5) + 5 = 0.167 \times 5 + 5 = 5.83 + ...` β€” approximately 0.33 units toward the prior mean. **(d)** The gap between informative and reference posterior means is :math:`\kappa_0(\bar{y} - \mu_0)/(\kappa_0 + n)`. For this to be less than 0.01: :math:`5 \times |7 - 5|/(\kappa_0 + n) < 0.01 \Rightarrow n > 5 \times 2/0.01 - 5 = 995`. So approximately :math:`n^* \approx 1000` observations are required. References ---------- **Conjugate Analysis and Foundations** .. [RaiffaSchlaifer1961b] Raiffa, H., and Schlaifer, R. (1961). *Applied Statistical Decision Theory*. Harvard Business School. The systematic catalog of conjugate prior families for exponential family likelihoods, providing the algebraic foundation for all five families developed in this section. .. [DeGroot1970] DeGroot, M. H. (1970). *Optimal Statistical Decisions*. McGraw-Hill. (Republished by Wiley, 2004.) Conjugate families derived from the perspective of optimal Bayesian decision theory; the source for the pseudo-data interpretation. .. [GelmanBDA3ch2] Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). *Bayesian Data Analysis* (3rd ed.), Ch. 2–3. Chapman and Hall/CRC. The authoritative modern treatment of conjugate families, including the Normal-Inverse-Gamma derivation in Β§3.2 that the present section follows. .. [Hoff2009] Hoff, P. D. (2009). *A First Course in Bayesian Statistical Methods*. Springer. A concise and rigorous treatment of conjugate families, hierarchical models, and MCMC at a level closely matched to this course. **Noninformative and Jeffreys Priors** .. [Jeffreys1946] Jeffreys, H. (1946). An invariant form for the prior probability in estimation problems. *Proceedings of the Royal Society of London A*, 186(1007), 453–461. The paper introducing the Fisher-information-based invariant prior, with derivations for the Normal, Binomial, and Poisson families. .. [Bernardo1979] Bernardo, J. M. (1979). Reference posterior distributions for Bayesian inference. *Journal of the Royal Statistical Society B*, 41(2), 113–147. Extends Jeffreys priors to the multiparameter case via reference analysis; the theoretical foundation for objective Bayes in modern practice. .. [KassWasserman1996] Kass, R. E., and Wasserman, L. (1996). The selection of prior distributions by formal rules. *Journal of the American Statistical Association*, 91(435), 1343–1370. A comprehensive review of formal rules for choosing noninformative priors, with discussion of Jeffreys, reference, and invariant priors. **Weakly Informative Priors and Prior Predictive Simulation** .. [Gelman2006] Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. *Bayesian Analysis*, 1(3), 515–534. Introduces the half-Cauchy(0, 25) as a weakly informative prior for hierarchical variance components; widely influential in applied Bayesian practice. .. [GelmanSimpsonBetancourt2017] Gelman, A., Simpson, D., and Betancourt, M. (2017). The prior can often only be understood in the context of the likelihood. *Entropy*, 19(10), 555. Argues that prior specification must be evaluated in the context of the full model; motivates prior predictive simulation as the primary evaluation tool. .. [McElreath2020] McElreath, R. (2020). *Statistical Rethinking: A Bayesian Course with Examples in R and Stan* (2nd ed.), Ch. 4. CRC Press. The source of the prior predictive simulation methodology for regression coefficients; the regression-line spaghetti plots in Figure 5.13 follow this text's pedagogical approach. .. [StanPriorWiki] Stan Development Team. (2024). Prior Choice Recommendations. *Stan Wiki*. Available at: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations. A living document of weakly informative prior recommendations for common model types, maintained by practitioners. **Sensitivity Analysis and Robustness** .. [Berger1985] Berger, J. O. (1985). *Statistical Decision Theory and Bayesian Analysis* (2nd ed.), Ch. 4. Springer. The theoretical foundation for Bayesian robustness and prior sensitivity analysis; the framework for prior sensitivity in Block 7d of this section. .. [Gustafson2015] Gustafson, P. (2015). *Bayesian Inference for Partially Identified Models: Exploring the Limits of Limited Data*. CRC Press. Prior sensitivity for non-identified and weakly identified models; relevant for the high-sensitivity regime discussed in Block 7d. **Modern Computational References** .. [Downey2021b] Downey, A. B. (2021). *Think Bayes: Bayesian Statistics in Python* (2nd ed.), Ch. 4–6. O'Reilly Media. Computational-first conjugate analysis using Python; the BetaBinomialModel class design is informed by this text's implementation philosophy. .. [GelmanHill2007] Gelman, A., and Hill, J. (2007). *Data Analysis Using Regression and Multilevel/Hierarchical Models*, Ch. 16. Cambridge University Press. Prior specification for regression models including the Normal weakly informative recommendations implemented in Block 7b. .. [Vehtari2021] Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and BΓΌrkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved :math:`\hat{R}` for assessing convergence of MCMC. *Bayesian Analysis*, 16(2), 667–718. While primarily about MCMC diagnostics, the introduction discusses how prior predictive simulation relates to model checking and the distinction between prior and likelihood information.