Bayesian Mixture-of-Normals Regression

Python · R  ·  Download Gibbs sampler

Model

Earnings are not symmetric. Most people sit inside a reasonably tight band, and a long thin group sits well below it — so a model whose errors form one bell curve is describing a shape the data does not have. What that mismatch actually costs is the question: it turns out to change very little about the average worker and a great deal about who ends up in the tail.

An ordinary linear regression assumes a single Normal error. This model replaces it with a mixture of normals on the disturbance — which can approximate essentially any density, removing the normality assumption and capturing the skew and fat tails a single Normal cannot (earnings are the textbook case). It is the continuous sibling of Geweke & Keane (1997): GK put a mixture of normals on the disturbance of a binary probit; here we observe yy directly, so the same flexible error applies to a plain regression.

yt=xt′β+εt,εt∼∑j=1mpj N(μj, σj2)y_t = x_t'\beta + \varepsilon_t, \qquad \varepsilon_t \sim \sum_{j=1}^m p_j\, N(\mu_j,\ \sigma_j^2)

Two different "mixtures of normals" — don't confuse them. Here (and in Geweke–Keane) the mixture is on the error εt\varepsilon_t → a flexible distribution in one equation, no hierarchy. In the coefficient-mixture examples (hierarchical MNL, hierarchical-LM mixtures) the mixture is on the coefficients βi\beta_i across units → heterogeneity / latent segments, which needs panel data. Different targets, same tool.

Identification: the error mixture is recentred to mean 0 each sweep (its location is absorbed into the intercept), so β\beta is the usual mean function and the mixture describes only the shape of the disturbance around it. Model comparison is by WAIC (GK used Bayes factors) — the widely applicable information criterion, which estimates how well a fitted model would predict observations it has not seen, by scoring the log-likelihood observation by observation and then deducting a penalty for effective complexity. Lower is better, and only the difference between two models carries meaning; the absolute level does not.

4-Block Gibbs sampler

BlockDrawFull conditional
(1) st∣⋅s_t \mid \cdotCategorical component indicator, P(j)∝pj N(rt; μj,σj2)P(j) \propto p_j\, N(r_t;\, \mu_j, \sigma_j^2)
(2) β∣⋅\beta \mid \cdotWeighted Gaussian (GLS, weights 1/σst21/\sigma^2_{s_t}) on y−μsty - \mu_{s_t}
(3) (μj,σj2)∣⋅(\mu_j, \sigma_j^2) \mid \cdotNormal-Inverse-Gamma per component
(4) p∣⋅p \mid \cdotDirichlet on the component counts

Closed-form predictive density p(y∣x)=∑jpj N(y; x′β+μj, σj2)p(y \mid x) = \sum_j p_j\, N(y;\, x'\beta + \mu_j,\, \sigma_j^2) — used for WAIC.

Notebooks

Section 1 validates the sampler on synthetic data (n=3000n=3000, true β=(1.0,0.5,−0.7)\beta=(1.0, 0.5, -0.7), a skewed two-component error): β\beta is recovered by every mm (it is the conditional mean), but WAIC drops ≈929\approx 929 from m=1m=1 to m=2m=2 and then flattens — correctly identifying the two-component truth. Section 2 fits PSID log-earnings (Geweke 2005; 2,733 positive earners) on age, age², and education: the OLS residual has kurtosis ≈23\approx 23, and the mixture wins decisively — WAIC falls ≈1,580\approx 1{,}580 from the normal. The conditional mean is stable (return to schooling ≈0.27\approx 0.27 per SD), but the error resolves into an ≈86%\approx 86\% tight majority component plus a wide ≈14%\approx 14\% low-earnings left tail a single Normal cannot represent. The R notebook confirms it two ways — flexmix BIC (selects k=2k=2) and bayesm's rnmixGibbs. A PyMC companion approaches the Normal-vs-mixture question from the other model-comparison angle: instead of WAIC (predictive accuracy, ≈\approx leave-one-out), it computes each model's marginal likelihood via Sequential Monte Carlo and forms a Bayes factor and posterior model probabilities. Three terms worth unpacking. A model's marginal likelihood is the probability it assigns to the observed data after averaging over every value its parameters could have taken — so a model that fits well only under finely tuned parameter values is penalised automatically, and complexity is charged for without any explicit penalty term. A Bayes factor is simply the ratio of two models' marginal likelihoods: the odds the data alone place on one against the other. Sequential Monte Carlo is how the marginal likelihood is obtained here — a population of particles is walked from the prior to the posterior through a sequence of intermediate, progressively tempered distributions, and the marginal likelihood falls out as a by-product of the reweighting at each step. SMC is validated against the closed-form Normal evidence, then the synthetic data give log⁡BF≈41\log\text{BF}\approx 41 (P(mixture∣y)≈1P(\text{mixture}\mid y)\approx 1) and PSID a log⁡BF\log\text{BF} in the hundreds. It also reports the honest caveat a mixture Bayes factor demands — the magnitude is prior-sensitive (the Bartlett–Lindley Occam penalty grows with the prior SD on the extra component — spreading a parameter's prior over a wider range dilutes the probability it puts near the values the data actually favour, so a vaguer prior on the extra component makes the mixture score worse no matter what the evidence says), but the conclusion stays robust across two orders of magnitude in the prior, and agrees emphatically with the WAIC verdict.

Python Mixture-of-normals error Gibbs sampler, synthetic validation + PSID earnings, WAIC (Python) Open → Python Bayes factors via SMC — marginal likelihood, posterior model probabilities, prior sensitivity (PyMC) Open → R bayesm::rnmixGibbs + flexmix BIC — R cross-check on the earnings residuals (R) Open →

Downloads

mixreg_gibbs.py From-scratch mixture-of-normals error Gibbs sampler + WAIC (NumPy / scipy) psid_earnings.csv PSID extract (Geweke 2005) — age, education, earnings

Gibbs Sampler — Source Code

"""
Mixture-of-normals regression — from-scratch Gibbs (NumPy + scipy.norm).

The continuous sibling of the Geweke-Keane (1997) mixture-of-normals probit: instead of a
binary outcome with a latent-utility step, we observe y directly and let the regression error
follow a flexible mixture of normals (capturing skew AND fat tails in, e.g., earnings):

    y_t = x_t' beta + eps_t,   eps_t ~ sum_j p_j N(mu_j, sigma_j^2)

Identification: the error mixture is recentred to mean 0 each sweep (its location goes into the
intercept), so beta is the usual regression mean function and the mixture describes the
*shape* of the disturbance around it.

Gibbs blocks:
  1. s_t         -- component indicator, categorical  P(j) ∝ p_j N(r_t; mu_j, sigma_j^2)
  2. beta        -- weighted Gaussian (GLS, weights 1/sigma_{s_t}^2) on y - mu_{s_t}
  3. (mu_j,sig2_j)-- Normal-Inverse-Gamma per component
  4. p           -- Dirichlet
Closed-form density:  p(y|x) = sum_j p_j N(y; x'beta + mu_j, sigma_j^2)  -> used for WAIC.
"""
import numpy as np
from scipy.stats import norm


def mixreg_gibbs(y, X, m=2, R=6000, burn=2000, seed=0,
                 prior_sd_beta=10.0, a0=3.0, kappa0=0.01, alpha_dir=1.0):
    rng = np.random.default_rng(seed)
    n, k = X.shape
    Hb = np.eye(k) / prior_sd_beta ** 2
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    resid = y - X @ beta
    vy = float(np.var(resid)); b0 = 0.5 * vy * (a0 - 1)        # prior mean sigma2 ~ var(resid)
    sig2 = np.full(m, vy)
    mu = np.linspace(-1.0, 1.0, m) * np.std(resid) if m > 1 else np.zeros(1)
    p = np.ones(m) / m
    s = rng.integers(0, m, n)

    keep = R - burn
    B = np.zeros((keep, k)); P = np.zeros((keep, m)); MU = np.zeros((keep, m)); S2 = np.zeros((keep, m))
    LL = np.zeros((keep, n))
    for g in range(R):
        r = y - X @ beta
        # 1. component indicators
        logpj = (np.log(p)[None, :] - 0.5 * np.log(sig2)[None, :]
                 - 0.5 * (r[:, None] - mu[None, :]) ** 2 / sig2[None, :])
        logpj -= logpj.max(1, keepdims=True)
        pp = np.exp(logpj); pp /= pp.sum(1, keepdims=True)
        s = (rng.random(n)[:, None] > np.cumsum(pp, 1)).sum(1).clip(0, m - 1)
        # 2. beta | . (weighted Gaussian on y - mu_{s})
        w = 1.0 / sig2[s]
        yc = y - mu[s]
        V = np.linalg.inv((X * w[:, None]).T @ X + Hb)
        beta = V @ ((X * w[:, None]).T @ yc) + np.linalg.cholesky(V) @ rng.standard_normal(k)
        # 3. (mu_j, sigma2_j) | . (Normal-Inverse-Gamma)
        r = y - X @ beta
        for j in range(m):
            mk = s == j; nj = int(mk.sum()); rj = r[mk]
            kn = kappa0 + nj
            mn = rj.sum() / kn                                  # prior mean 0
            an = a0 + nj / 2.0
            bn = b0 + 0.5 * (np.sum((rj - rj.mean()) ** 2) if nj > 0 else 0.0) \
                 + 0.5 * (kappa0 * nj / kn) * (rj.mean() if nj > 0 else 0.0) ** 2
            sig2[j] = 1.0 / rng.gamma(an, 1.0 / bn)
            mu[j] = mn + np.sqrt(sig2[j] / kn) * rng.standard_normal()
        # 4. p | .
        p = rng.dirichlet(alpha_dir + np.bincount(s, minlength=m))
        # identification: recentre mixture to mean 0 (location -> intercept), then sort by mu
        mbar = np.sum(p * mu); beta[0] += mbar; mu -= mbar
        o = np.argsort(mu); mu = mu[o]; sig2 = sig2[o]; p = p[o]; s = np.argsort(o)[s]
        if g >= burn:
            gg = g - burn; B[gg] = beta; P[gg] = p; MU[gg] = mu; S2[gg] = sig2
            dens = (p[None, :] * norm.pdf(r[:, None], mu[None, :], np.sqrt(sig2)[None, :])).sum(1)
            LL[gg] = np.log(np.clip(dens, 1e-300, None))
    return dict(beta=B, p=P, mu=MU, sig2=S2, ll=LL, m=m)


def waic(out):
    ll = out['ll']; mx = ll.max(0)
    lppd = np.log(np.exp(ll - mx).mean(0)) + mx
    return float(-2 * (lppd.sum() - ll.var(0).sum()))


def error_density(out, grid):
    """Posterior-mean fitted error density on `grid`."""
    p = out['p'].mean(0); mu = out['mu'].mean(0); sd = np.sqrt(out['sig2'].mean(0))
    return sum(p[j] * norm.pdf(grid, mu[j], sd[j]) for j in range(out['m']))

References