Bayesian Mixture-of-Normals Probit

Python · R  ·  Download Gibbs sampler

Model

In a binary model the distribution of the disturbance is the link function: P(y=1∣x)=F(x′β)P(y=1\mid x) = F(x'\beta) where FF is the error CDF. Probit takes F=ΦF=\Phi (Normal), logit takes the logistic — both symmetric. Geweke & Keane (1997) let FF be a mixture of normals, which can approximate essentially any shape, so the link is estimated rather than assumed. This is the binary sibling of the mixture-of-normals regression, and it generalizes both Albert–Chib probit and the robit (a Student-t scale mixture).

yt=1[ xt′β+εt>0 ],εt∼∑j=1mpj N(αj, hj−1)y_t = \mathbf{1}[\,x_t'\beta + \varepsilon_t > 0\,], \qquad \varepsilon_t \sim \sum_{j=1}^m p_j\, N(\alpha_j,\ h_j^{-1})

Two variants. A scale mixture (αj=0\alpha_j = 0) gives a symmetric, fat-tailed link that generalizes the robit. A full (location-scale) mixture with free component means αj\alpha_j gives an asymmetric link — the right tool when the dose–response curve is skewed, used in the beetle application below.

Identification: a probit fixes the error location and scale, so each sweep the mixture is recentred to mean 0 (location absorbed into the intercept) and rescaled to Var(ε)=1\mathrm{Var}(\varepsilon)=1 (the probit scale), then components are sorted by mean. The choice probability has a closed form, P(y=1∣x)=∑jpj Φ ⁣(hj (x′β+αj))P(y=1\mid x) = \sum_j p_j\,\Phi\!\big(\sqrt{h_j}\,(x'\beta+\alpha_j)\big) — used for WAIC, the widely applicable information criterion, which scores how well a fitted model would predict observations it has not seen and charges a penalty for effective complexity. Lower is better and only differences between models mean anything.

5-Block Gibbs sampler

BlockDrawFull conditional
(1) y~t∣⋅\tilde y_t \mid \cdotTruncated Normal (Albert–Chib), sd 1/hst1/\sqrt{h_{s_t}}
(2) st∣⋅s_t \mid \cdotCategorical component indicator, P(j)∝pjhj e−12hjrt2P(j) \propto p_j \sqrt{h_j}\, e^{-\frac12 h_j r_t^2}
(3) β∣⋅\beta \mid \cdotWeighted Gaussian (GLS, weights hsth_{s_t})
(4) (αj,hj)∣⋅(\alpha_j, h_j) \mid \cdotNormal-Inverse-Gamma per component (Gamma precision for the scale mixture)
(5) p∣⋅p \mid \cdotDirichlet on the component counts

Notebooks

Section 1 validates the sampler on synthetic data (n=4000n=4000, a skewed two-component error): the full mixture recovers the right-skew, de-biases the slope (1.45→1.331.45 \to 1.33 vs. truth 1.2), and cuts WAIC by ≈109\approx 109 from the symmetric probit. Section 2 fits the Bliss (1935) beetle dose–mortality data — the dataset that birthed the probit, and the textbook case where symmetric links underfit a skewed mortality curve. The flexible-error probit discovers the asymmetry on its own: a 2-component mixture cuts the grouped deviance — a badness-of-fit score measuring how far the fitted probabilities sit from the observed death rates, which is roughly chi-squared under a correct model, so with 6 dose groups a value near 3 is a good fit and one near 10 is a poor one — from 10.1 (conventional probit) to 3.40, essentially equal to the cloglog link (3.45) statisticians historically chose by hand. Cloglog, the complementary log-log, is the standard asymmetric link: it approaches certainty of death faster than it leaves certainty of survival, which is what a dose-response curve for a poison actually looks like. The R notebook confirms it frequentist-side: a symmetric Gosset t-link stays at ≈10.2\approx 10.2 (the problem is skew, not tail weight), while the asymmetric Aranda–Ordaz family — a one-parameter family of links containing the logit at λ=1\lambda=1 and the cloglog at λ=0\lambda=0, so the shape is estimated rather than picked — drives λ→0\lambda \to 0 (cloglog), deviance 3.44. Whether you estimate the link or estimate the error distribution, the data tells one story.

Python Scale & full mixture-of-normals probit Gibbs, synthetic + Bliss beetle, WAIC/deviance (Python) Open → R Estimated flexible links — glm probit/logit/cloglog, Gosset t-link, Aranda–Ordaz (R) Open →

Downloads

mixprobit_gibbs.py From-scratch scale & full mixture-of-normals probit Gibbs samplers + WAIC (NumPy / scipy) bliss_beetle.csv Bliss (1935) beetle dose–mortality — 8 doses of CS₂, killed / total per dose (481 beetles)

Gibbs Sampler — Source Code

"""
Mixture-of-normals probit (Geweke & Keane 1997) — from-scratch Gibbs (NumPy + scipy.norm).

Binary probit with a SCALE MIXTURE OF NORMALS disturbance (zero-mean components, mixing on
precision — symmetric, leptokurtic; the variant GK found best for PSID LFP):

    y_t = 1[ x_t' beta + eps_t > 0 ],   eps_t ~ sum_j p_j N(0, h_j^{-1})

Extends Albert-Chib (single-normal probit) and robit (t = a scale mixture) to a flexible,
semiparametric symmetric error. Gibbs blocks:
  1. latent y*_t           -- truncated normal (Albert-Chib), sd = 1/sqrt(h_{s_t})
  2. component indicator s_t -- categorical, P(j) ∝ p_j sqrt(h_j) exp(-.5 h_j r_t^2)
  3. beta                  -- weighted Gaussian (GLS with weights h_{s_t})
  4. h_j                   -- Gamma (conjugate precision per component)
  5. p                     -- Dirichlet
Identification: rescale each sweep so Var(eps)=sum_j p_j/h_j = 1 (the probit scale, as in
conventional probit); components sorted by h to fix labeling.

Choice probability has closed form:  P(y=1|x) = sum_j p_j * Phi( sqrt(h_j) * x'beta ).
"""
import numpy as np
from scipy.stats import norm


def _rtruncnorm(mu, sd, d, rng):
    """Truncated normal: d=1 -> (0, inf), d=0 -> (-inf, 0). Inverse-CDF, vectorized."""
    thr = norm.cdf(-mu / sd)                       # Phi at the 0-boundary (standardized)
    a = np.where(d == 1, thr, 0.0)
    b = np.where(d == 1, 1.0, thr)
    u = np.clip(a + (b - a) * rng.random(len(mu)), 1e-12, 1 - 1e-12)
    return mu + sd * norm.ppf(u)


def simulate_mixprobit(n, beta, error='t', df=4, seed=0):
    """Synthetic probit data with a chosen (variance-1) error. error: 'normal','t','mix2'."""
    rng = np.random.default_rng(seed)
    beta = np.asarray(beta, float); k = len(beta)
    X = np.column_stack([np.ones(n)] + [rng.normal(size=n) for _ in range(k - 1)])
    if error == 'normal':
        eps = rng.standard_normal(n)
    elif error == 't':
        eps = rng.standard_t(df, n) / np.sqrt(df / (df - 2))      # standardized to var 1
    elif error == 'mix2':                                          # leptokurtic 2-component
        z = rng.random(n) < 0.85
        raw = np.where(z, rng.normal(0, 0.6, n), rng.normal(0, 2.2, n))
        eps = raw / raw.std()
    ystar = X @ beta + eps
    return X, (ystar > 0).astype(float), eps


def mixprobit_gibbs(y, X, m=2, R=6000, burn=2000, seed=0,
                    prior_sd_beta=10.0, a0=2.0, b0=1.0, alpha_dir=1.0):
    """Scale-mixture-of-normals probit via Gibbs. m = number of mixture components."""
    rng = np.random.default_rng(seed)
    n, k = X.shape
    Hb = np.eye(k) / prior_sd_beta ** 2
    beta = np.zeros(k)
    h = np.linspace(0.6, 1.8, m) if m > 1 else np.array([1.0])
    p = np.ones(m) / m
    s = rng.integers(0, m, n)

    keep = R - burn
    B = np.zeros((keep, k)); P = np.zeros((keep, m)); H = np.zeros((keep, m)); LL = np.zeros((keep, n))
    for g in range(R):
        mu = X @ beta
        # 1. latent utilities (Albert-Chib), per-obs sd from its component
        yl = _rtruncnorm(mu, 1.0 / np.sqrt(h[s]), y, rng)
        # 2. component indicators
        r = yl - mu
        logpj = np.log(p)[None, :] + 0.5 * np.log(h)[None, :] - 0.5 * h[None, :] * r[:, None] ** 2
        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)
        # 3. beta | . (weighted Gaussian)
        w = h[s]
        V = np.linalg.inv((X * w[:, None]).T @ X + Hb)
        beta = V @ ((X * w[:, None]).T @ yl) + np.linalg.cholesky(V) @ rng.standard_normal(k)
        # 4. h_j | . (Gamma)
        r = yl - X @ beta
        for j in range(m):
            mk = s == j; nj = int(mk.sum())
            h[j] = rng.gamma(a0 + nj / 2.0, 1.0 / (b0 + 0.5 * np.sum(r[mk] ** 2))) if nj > 0 \
                else rng.gamma(a0, 1.0 / b0)
        # 5. p | . (Dirichlet)
        p = rng.dirichlet(alpha_dir + np.bincount(s, minlength=m))
        # identification: fix Var(eps)=1, then sort components by precision
        v = np.sum(p / h); beta = beta / np.sqrt(v); h = h * v
        o = np.argsort(h); h = h[o]; p = p[o]; s = np.argsort(o)[s]
        if g >= burn:
            gg = g - burn; B[gg] = beta; P[gg] = p; H[gg] = h
            P1 = np.clip((p[None, :] * norm.cdf(np.sqrt(h)[None, :] * (X @ beta)[:, None])).sum(1), 1e-12, 1 - 1e-12)
            LL[gg] = np.where(y == 1, np.log(P1), np.log(1 - P1))
    return dict(beta=B, p=P, h=H, ll=LL, m=m)


def mixprobit_full_gibbs(y, X, m=2, R=8000, burn=3000, seed=0,
                         prior_sd_beta=10.0, a0=4.0, kappa0=0.05, alpha_dir=2.0):
    """FULL (location-scale) mixture-of-normals probit:  eps ~ sum_j p_j N(alpha_j, h_j^{-1}).

    Unlike the scale mixture above (symmetric), the free component means alpha_j let the
    disturbance be ASYMMETRIC -- a flexible/skewed link, the right tool for skewed dose-response
    curves (e.g. the Bliss beetle data, where probit/logit underfit but cloglog fits).

    Identification each sweep: recentre the mixture to mean 0 (location -> intercept) AND rescale
    to Var(eps)=sum_j p_j(alpha_j^2 + 1/h_j)=1 (the probit scale); components sorted by mean.
    Choice prob:  P(y=1|x) = sum_j p_j * Phi( sqrt(h_j) * (x'beta + alpha_j) )."""
    rng = np.random.default_rng(seed)
    n, k = X.shape
    Hb = np.eye(k) / prior_sd_beta ** 2
    b0 = 0.5 * (a0 - 1)                                   # prior mean component var ~ 0.5
    beta = np.zeros(k)
    alpha = np.linspace(-0.5, 0.5, m) if m > 1 else np.zeros(1)
    h = np.ones(m); p = np.ones(m) / m
    s = rng.integers(0, m, n)

    keep = R - burn
    B = np.zeros((keep, k)); P = np.zeros((keep, m)); A = np.zeros((keep, m))
    H = np.zeros((keep, m)); LL = np.zeros((keep, n))
    for g in range(R):
        mu = X @ beta
        yl = _rtruncnorm(mu + alpha[s], 1.0 / np.sqrt(h[s]), y, rng)
        r = yl - mu
        logpj = np.log(p)[None, :] + 0.5 * np.log(h)[None, :] - 0.5 * h[None, :] * (r[:, None] - alpha[None, :]) ** 2
        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)
        # beta | . (weighted GLS on yl - alpha_s)
        w = h[s]
        V = np.linalg.inv((X * w[:, None]).T @ X + Hb)
        beta = V @ ((X * w[:, None]).T @ (yl - alpha[s])) + np.linalg.cholesky(V) @ rng.standard_normal(k)
        # (alpha_j, h_j) | . (Normal-Inverse-Gamma)
        r = yl - X @ beta
        for j in range(m):
            mk = s == j; nj = int(mk.sum()); rj = r[mk]
            kn = kappa0 + nj; mn = rj.sum() / kn; 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
            h[j] = rng.gamma(an, 1.0 / bn)
            alpha[j] = mn + np.sqrt(1.0 / (kn * h[j])) * rng.standard_normal()
        p = rng.dirichlet(alpha_dir + np.bincount(s, minlength=m))
        # identification: recentre to mean 0, rescale to Var=1, sort by mean
        mbar = np.sum(p * alpha); beta[0] += mbar; alpha -= mbar
        sv = np.sqrt(np.sum(p * (alpha ** 2 + 1.0 / h)))
        beta = beta / sv; alpha = alpha / sv; h = h * sv ** 2
        o = np.argsort(alpha); alpha = alpha[o]; h = h[o]; p = p[o]; s = np.argsort(o)[s]
        if g >= burn:
            gg = g - burn; B[gg] = beta; P[gg] = p; A[gg] = alpha; H[gg] = h
            P1 = np.clip((p[None, :] * norm.cdf(np.sqrt(h)[None, :] * ((X @ beta)[:, None] + alpha[None, :]))).sum(1),
                         1e-12, 1 - 1e-12)
            LL[gg] = np.where(y == 1, np.log(P1), np.log(1 - P1))
    return dict(beta=B, p=P, alpha=A, h=H, ll=LL, m=m)


def prob_curve(out, X):
    """Posterior-mean P(y=1|X) for a full-mixture-probit fit, evaluated at rows of X."""
    b = out['beta'].mean(0); p = out['p'].mean(0); a = out['alpha'].mean(0); h = out['h'].mean(0)
    eta = X @ b
    return (p[None, :] * norm.cdf(np.sqrt(h)[None, :] * (eta[:, None] + a[None, :]))).sum(1)


def waic(out):
    """WAIC from pointwise log-likelihood draws (lower = better)."""
    ll = out['ll']; mx = ll.max(0)
    lppd = (np.log(np.exp(ll - mx).mean(0)) + mx)
    pw = ll.var(0)
    return float(-2 * (lppd.sum() - pw.sum()))


def error_density(out, grid):
    """Posterior-mean fitted disturbance density on `grid` (the estimated shock distribution)."""
    p = out['p'].mean(0); h = out['h'].mean(0)
    return sum(p[j] * norm.pdf(grid, 0, 1 / np.sqrt(h[j])) for j in range(out['m']))


def error_density_full(out, grid):
    """Fitted disturbance density for a FULL-mixture fit (components have free means alpha_j)."""
    p = out['p'].mean(0); a = out['alpha'].mean(0); h = out['h'].mean(0)
    return sum(p[j] * norm.pdf(grid, a[j], 1 / np.sqrt(h[j])) for j in range(out['m']))

References