Bayesian Student-t Linear Model

Python · R · Download Gibbs sampler

Model

The standard linear model assumes Gaussian errors. When data contain occasional large shocks — outliers from economic crises, measurement error, fat-tailed processes — the Normal assumption overstates the influence of extreme observations on β\beta. Replacing N(0,σ2)\mathcal{N}(0, \sigma^2) errors with a Student-tt distribution provides automatic robustness: the heavier tails absorb outliers without requiring any manual flagging.

yi=xi′β+εi,εi∼i.i.d.t(ν, 0, σ2)y_i = x_i'\beta + \varepsilon_i, \qquad \varepsilon_i \stackrel{\text{i.i.d.}}{\sim} t(\nu,\, 0,\, \sigma^2)

Scale-Mixture Representation

The key to tractable Bayesian inference is the scale-mixture-of-normals representation of the Student-tt. Introducing a latent precision weight λi\lambda_i per observation makes all four full conditionals conjugate. Outlier observations receive low λi\lambda_i (high local variance σ2/λi\sigma^2/\lambda_i), limiting their leverage on β\beta automatically.

yi∣λi∼N ⁣(xi′β,  σ2λi),λi∣ν∼Gamma ⁣(ν2, ν2)y_i \mid \lambda_i \sim N\!\left(x_i'\beta,\; \frac{\sigma^2}{\lambda_i}\right), \qquad \lambda_i \mid \nu \sim \text{Gamma}\!\left(\tfrac{\nu}{2},\, \tfrac{\nu}{2}\right)

The degrees-of-freedom parameter ν\nu controls tail heaviness: ν→∞\nu \to \infty recovers the Normal; ν=4\nu = 4 has substantially heavier tails. Rather than reparametrising ν\nu for continuous HMC, Geweke places a discrete uniform prior on a grid {2,3,…,40}\{2, 3, \ldots, 40\} and draws ν\nu as a categorical variable each iteration — fast and numerically stable near the ν=2\nu = 2 boundary.

4-Block Gibbs — Geweke (1993)

BlockDrawDistribution
(1) β∣λ,σ2,y\beta \mid \lambda, \sigma^2, yNormal conjugate — WLS posterior N(bˉ,Bˉ)N(\bar{b}, \bar{B}) with Bˉ−1=X′ΛX/σ2\bar{B}^{-1} = X'\Lambda X / \sigma^2
(2) σ2∣β,λ,y\sigma^2 \mid \beta, \lambda, yInverse-Gamma conjugate IG(an/2,  sn/2)\text{IG}(a_n/2,\; s_n/2) where sn=s0+∑λiei2s_n = s_0 + \textstyle\sum \lambda_i e_i^2
(3) λi∣β,σ2,ν,y\lambda_i \mid \beta, \sigma^2, \nu, yGamma conjugate Gamma ⁣(ν+12,  ν+ei2/σ22)\text{Gamma}\!\left(\tfrac{\nu+1}{2},\; \tfrac{\nu + e_i^2/\sigma^2}{2}\right) independently
(4) ν∣λ\nu \mid \lambdaDiscrete draw on ν\nu-grid via Gamma log-likelihood

Notebooks

Applied to two datasets. Section 1: synthetic contaminated data (n=120n=120, ν⋆=4\nu^\star=4), where the sampler recovers the true degrees of freedom and correctly identifies downweighted outliers. Section 2: Nelson-Plosser (1982) 14 annual U.S. macroeconomic series (1909–1970), where low posterior ν\nu for series like Real GNP, Industrial Production, and Employment corresponds to Depression- and WWII-era shocks — observations the model down-weights rather than attributes to a stochastic trend. Section 3: PyMC comparison with continuous ν\nu via NUTS on both datasets. Two companion notebooks cross-check the model against general-purpose tools. The PyMC notebook treats ν\nu as continuous with a weakly-informative Gamma prior and samples (β,σ,ν)(\beta, \sigma, \nu) jointly by NUTS — the β\beta posteriors match the discrete-ν Gibbs to within a posterior SD, with the trade-off that NUTS marginalises out the latent λi\lambda_i (so it loses the per-observation outlier weights the scale-mixture Gibbs exposes directly). The R notebook adds the frequentist counterparts — hett::tlm (Student-t by maximum likelihood, recovering ν≈3.9\nu\approx 3.9 on the synthetic t(4)t(4) data), MASS::rlm (Huber M-estimation), and lm (the non-robust OLS baseline that gets dragged toward the outliers). Gibbs, NUTS, and ML all agree the macro series are heavy-tailed.

Python Student-t Gibbs sampler (discrete ν), Nelson-Plosser application (Python) Open → Python PyMC cross-check — continuous ν via pm.StudentT / NUTS, vs the discrete-ν Gibbs (Python) Open → R R cross-check — hett::tlm (Student-t ML), MASS::rlm (robust), lm (OLS baseline) (R) Open →

Downloads

student_t_lm_gibbs.py 4-block Gibbs sampler module nelplo.csv Nelson-Plosser (1982) macroeconomic series — 14 annual U.S. series, 1909–1970

Gibbs Sampler — Source Code

"""
student_t_lm_gibbs.py

Bayesian inference for the linear model with i.i.d. Student-t errors.
Implements Geweke (1993) "Bayesian treatment of the independent Student-t
linear model."

Model
-----
    y_i = x_i' β + ε_i,    ε_i ~ i.i.d. t(ν, 0, σ²)

Scale-mixture representation (Student-t = Normal scale-mixture)
---------------------------------------------------------------
    y_i | λ_i  ~  N(x_i' β,  σ²/λ_i)
    λ_i | ν    ~  Gamma(ν/2, ν/2)          [shape / rate]
    ε_i  marginally  ~  t(ν, 0, σ²)

Four-block Gibbs sampler
------------------------
    (1)  β    | λ, σ², y  ~  N(b̄, B̄)          WLS posterior
    (2)  σ²   | β, λ, y   ~  IG(aₙ/2, sₙ/2)
    (3)  λᵢ   | β, σ², ν, y  ~  Gamma((ν+1)/2,  (ν + eᵢ²/σ²)/2)
    (4)  ν    | λ          ~  discrete on ν_grid  (via Gamma log-likelihood)

Default priors
--------------
    β   ~  flat  (or N(b0, B0) if supplied)
    σ²  ~  IG(a0/2, s0/2)       flat when a0=s0=0
    ν   ~  discrete uniform on {2, 3, …, 40}

Reference
---------
    Geweke, J. (1993). "Bayesian treatment of the independent Student-t
    linear model." Journal of Applied Econometrics, 8(S1), S19-S40.
"""

import time
import numpy as np
from scipy.special import gammaln


# ── public API ────────────────────────────────────────────────────────────────

def student_t_gibbs(
    y, X,
    nu_grid=None,
    b0=None, B0=None,
    a0=0.0, s0=0.0,
    R=10000, burn=2000,
    seed=42, nprint=1000,
):
    """
    Four-block Gibbs sampler for the Student-t linear model.

    Parameters
    ----------
    y       : (n,)        observations
    X       : (n, k)      design matrix (include intercept column if desired)
    nu_grid : (G,) array  grid of candidate ν values  [default: 2, 3, …, 40]
    b0      : (k,)        prior mean for β             [flat if omitted]
    B0      : (k, k)      prior covariance for β       [flat if omitted]
    a0, s0  : float       IG(a0/2, s0/2) prior for σ²  [flat if both 0]
    R, burn : int         total and burn-in iterations
    seed    : int
    nprint  : int         print progress every nprint iterations (0 = silent)

    Returns
    -------
    dict with keys
      'beta'   : (R_kept, k)  posterior draws for β
      'sigma2' : (R_kept,)    posterior draws for σ²
      'lam'    : (R_kept, n)  posterior draws for the mixing weights λᵢ
      'nu'     : (R_kept,)    posterior draws for ν  (discrete)
    """
    rng   = np.random.default_rng(seed)
    y     = np.asarray(y, float).ravel()
    X     = np.asarray(X, float)
    n, k  = X.shape

    if nu_grid is None:
        nu_grid = np.arange(2.0, 41.0)
    nu_grid = np.asarray(nu_grid, float)
    G = len(nu_grid)

    # ── prior ─────────────────────────────────────────────────────────────────
    if b0 is not None and B0 is not None:
        B0_inv = np.linalg.inv(np.asarray(B0, float))
        b0     = np.asarray(b0, float)
    else:
        B0_inv = None

    # ── initial values ────────────────────────────────────────────────────────
    beta   = np.linalg.lstsq(X, y, rcond=None)[0]
    e      = y - X @ beta
    sigma2 = max(float(e @ e) / max(n - k, 1), 1e-8)
    lam    = np.ones(n)
    nu     = float(nu_grid[G // 2])

    # ── pre-compute ν-dependent log-normaliser (constant across iterations) ──
    # log p(λ | ν) = n × [(ν/2) log(ν/2) − logΓ(ν/2)]
    #              + (ν/2 − 1) Σ log λᵢ  −  ν/2 Σ λᵢ
    half_nu    = nu_grid / 2.0
    log_norm   = n * (half_nu * np.log(half_nu) - gammaln(half_nu))  # (G,)

    # ── storage ───────────────────────────────────────────────────────────────
    R_kept       = R - burn
    beta_draws   = np.empty((R_kept, k))
    sigma2_draws = np.empty(R_kept)
    lam_draws    = np.empty((R_kept, n))
    nu_draws     = np.empty(R_kept)
    kept = 0
    t0   = time.time()

    for r in range(1, R + 1):

        # ── Block 1: β | λ, σ², y  (WLS posterior) ───────────────────────────
        XtLX = X.T @ (lam[:, None] * X)   # (k, k)  X'ΛX
        XtLy = X.T @ (lam * y)            # (k,)    X'Λy

        if B0_inv is not None:
            Prec = XtLX / sigma2 + B0_inv
            rhs  = XtLy / sigma2 + B0_inv @ b0
        else:
            Prec = XtLX / sigma2
            rhs  = XtLy / sigma2

        Bbar   = np.linalg.inv(Prec)
        bbar   = Bbar @ rhs
        L_chol = np.linalg.cholesky(Bbar)
        beta   = bbar + L_chol @ rng.standard_normal(k)

        # ── Block 2: σ² | β, λ, y  (IG conjugate) ───────────────────────────
        e      = y - X @ beta
        SS_w   = float(lam @ (e * e))     # Σ λᵢ eᵢ²
        a_n    = a0 + n
        s_n    = s0 + SS_w
        sigma2 = s_n / rng.chisquare(a_n)

        # ── Block 3: λᵢ | β, σ², ν, y  ~ Gamma((ν+1)/2, (ν+eᵢ²/σ²)/2) ─────
        shape_lam = (nu + 1.0) / 2.0
        rate_lam  = (nu + e * e / sigma2) / 2.0   # (n,)
        lam       = rng.gamma(shape_lam, 1.0 / rate_lam)

        # ── Block 4: ν | λ  (discrete draw) ──────────────────────────────────
        sum_log_lam = float(np.log(lam).sum())
        sum_lam     = float(lam.sum())
        log_p = log_norm + (half_nu - 1.0) * sum_log_lam - half_nu * sum_lam
        log_p -= log_p.max()       # shift for numerical stability
        p = np.exp(log_p)
        p /= p.sum()
        nu = float(rng.choice(nu_grid, p=p))

        # ── store ─────────────────────────────────────────────────────────────
        if r > burn:
            beta_draws[kept]   = beta
            sigma2_draws[kept] = sigma2
            lam_draws[kept]    = lam
            nu_draws[kept]     = nu
            kept += 1

        if nprint > 0 and r % nprint == 0:
            elapsed = time.time() - t0
            print(f'  Iter {r:5d}/{R}  ({elapsed:.1f}s elapsed)')

    elapsed = time.time() - t0
    print(f'Done in {elapsed / 60:.1f} min  |  kept draws: {R_kept}')
    return {
        'beta':   beta_draws,
        'sigma2': sigma2_draws,
        'lam':    lam_draws,
        'nu':     nu_draws,
    }


def normal_gibbs(y, X, b0=None, B0=None, a0=0.0, s0=0.0,
                 R=10000, burn=2000, seed=42, nprint=0):
    """
    Gibbs sampler for the Normal linear model (ν → ∞ limit of Student-t).

    Identical to student_t_gibbs but with λᵢ ≡ 1 (no mixing weights).
    Used as a baseline for model comparison.

    Returns
    -------
    dict with keys 'beta' and 'sigma2'.
    """
    rng   = np.random.default_rng(seed)
    y     = np.asarray(y, float).ravel()
    X     = np.asarray(X, float)
    n, k  = X.shape

    if b0 is not None and B0 is not None:
        B0_inv = np.linalg.inv(np.asarray(B0, float))
        b0     = np.asarray(b0, float)
    else:
        B0_inv = None

    beta   = np.linalg.lstsq(X, y, rcond=None)[0]
    e      = y - X @ beta
    sigma2 = max(float(e @ e) / max(n - k, 1), 1e-8)

    R_kept       = R - burn
    beta_draws   = np.empty((R_kept, k))
    sigma2_draws = np.empty(R_kept)
    kept = 0
    t0   = time.time()

    XtX = X.T @ X
    Xty = X.T @ y

    for r in range(1, R + 1):
        if B0_inv is not None:
            Prec = XtX / sigma2 + B0_inv
            rhs  = Xty / sigma2 + B0_inv @ b0
        else:
            Prec = XtX / sigma2
            rhs  = Xty / sigma2

        Bbar   = np.linalg.inv(Prec)
        bbar   = Bbar @ rhs
        L_chol = np.linalg.cholesky(Bbar)
        beta   = bbar + L_chol @ rng.standard_normal(k)

        e      = y - X @ beta
        a_n    = a0 + n
        s_n    = s0 + float(e @ e)
        sigma2 = s_n / rng.chisquare(a_n)

        if r > burn:
            beta_draws[kept]   = beta
            sigma2_draws[kept] = sigma2
            kept += 1

        if nprint > 0 and r % nprint == 0:
            print(f'  Iter {r:5d}/{R}  ({time.time()-t0:.1f}s elapsed)')

    if nprint > 0:
        print(f'Done in {(time.time()-t0)/60:.1f} min  |  kept draws: {R_kept}')
    return {'beta': beta_draws, 'sigma2': sigma2_draws}


def posterior_summary(draws, param_names=None, true_vals=None):
    """
    Posterior mean, SD, and 95% credible interval for each column.

    Parameters
    ----------
    draws       : (R, p) array of posterior draws
    param_names : list of p strings       [optional]
    true_vals   : (p,) true values; adds 'true' and 'covered' columns

    Returns
    -------
    pandas.DataFrame
    """
    import pandas as pd
    p = draws.shape[1]
    if param_names is None:
        param_names = [f'b{j}' for j in range(p)]
    lo = np.quantile(draws, 0.025, axis=0)
    hi = np.quantile(draws, 0.975, axis=0)
    df = pd.DataFrame({
        'mean': draws.mean(0),
        'sd':   draws.std(0),
        'q025': lo,
        'q975': hi,
    }, index=param_names)
    if true_vals is not None:
        tv = np.asarray(true_vals)
        df['true']    = tv
        df['covered'] = (tv >= lo) & (tv <= hi)
    return df.round(4)


def log_lik_t(y, X, beta_draws, sigma2_draws, nu_draws):
    """
    Posterior mean log-likelihood under the Student-t model.

    Evaluates log p(y | βᵣ, σ²ᵣ, νᵣ) for each draw r and returns the mean.
    Used for informal model comparison (not a proper Bayes factor).
    """
    from scipy.stats import t as t_dist
    R = len(nu_draws)
    lls = np.empty(R)
    for r in range(R):
        e = y - X @ beta_draws[r]
        lls[r] = t_dist.logpdf(
            e, df=nu_draws[r], loc=0, scale=np.sqrt(sigma2_draws[r])
        ).sum()
    return float(lls.mean()), float(lls.std())


def log_lik_normal(y, X, beta_draws, sigma2_draws):
    """Posterior mean log-likelihood under the Normal model."""
    from scipy.stats import norm
    R = len(sigma2_draws)
    lls = np.empty(R)
    for r in range(R):
        e = y - X @ beta_draws[r]
        lls[r] = norm.logpdf(e, loc=0, scale=np.sqrt(sigma2_draws[r])).sum()
    return float(lls.mean()), float(lls.std())

References