Bayesian Mixture-of-Normals Probit

Model

In a binary model the distribution of the disturbance is the link function: $P(y=1\mid x) = F(x'\beta)$ where $F$ is the error CDF. Probit takes $F=\Phi$ (Normal), logit takes the logistic — both symmetric. Geweke & Keane (1997) let $F$ 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).

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})

$\beta$ — coefficients (the mean function); Normal prior
$(p_j, \alpha_j, h_j)$ — weight, location, precision of mixture component $j$
$m$ — number of components; chosen by WAIC / deviance

Two variants. A scale mixture ( $\alpha_j = 0$ ) gives a symmetric, fat-tailed link that generalizes the robit. A full (location-scale) mixture with free component means $\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 $\mathrm{Var}(\varepsilon)=1$ (the probit scale), then components are sorted by mean. The choice probability has a closed form, $P(y=1\mid x) = \sum_j p_j\,\Phi\!\big(\sqrt{h_j}\,(x'\beta+\alpha_j)\big)$ — used for WAIC.

5-Block Gibbs sampler

Block	Draw	Full conditional
(1)	$\tilde y_t \mid \cdot$	Truncated Normal (Albert–Chib), sd $1/\sqrt{h_{s_t}}$
(2)	$s_t \mid \cdot$	Categorical component indicator, $P(j) \propto p_j \sqrt{h_j}\, e^{-\frac12 h_j r_t^2}$
(3)	$\beta \mid \cdot$	Weighted Gaussian (GLS, weights $h_{s_t}$ )
(4)	$(\alpha_j, h_j) \mid \cdot$	Normal-Inverse-Gamma per component (Gamma precision for the scale mixture)
(5)	$p \mid \cdot$	Dirichlet on the component counts

Notebooks

Section 1 validates the sampler on synthetic data ( $n=4000$ , a skewed two-component error): the full mixture recovers the right-skew, de-biases the slope ( $1.45 \to 1.33$ vs. truth 1.2), and cuts WAIC by $\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 from 10.1 (conventional probit) to 3.40 — essentially equal to the cloglog link (3.45) statisticians historically chose by hand. The R notebook confirms it frequentist-side: a symmetric Gosset t-link stays at $\approx 10.2$ (the problem is skew, not tail weight), while the asymmetric Aranda–Ordaz family drives $\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

Geweke, J. & Keane, M. (1997). Mixture of Normals Probit Models. Federal Reserve Bank of Minneapolis Staff Report 237. — the flexible mixture-of-normals link this example implements
Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88(422), 669–679. — the data-augmentation truncated-normal step (Gibbs block 1)
Bliss, C. I. (1935). The calculation of the dosage–mortality curve. Annals of Applied Biology 22(1), 134–167. — the beetle dose–mortality data, and the origin of probit analysis
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research 11, 3571–3594. — WAIC, used to compare the number of mixture components
Aranda-Ordaz, F. J. (1981). On two families of transformations to additivity for binary response data. Biometrika 68(2), 357–363. — the asymmetric link family (λ→0 gives cloglog) used as the R frequentist cross-check
Koenker, R. & Yoon, J. (2009). Parametric links for binary choice models: a Fisherian–Bayesian colloquy. Journal of Econometrics 152(2), 120–130. — the Gosset and Pregibon link families behind R's glmx