Mixture-of-normals probit (Geweke & Keane 1997) — a flexible link for the Bliss beetle data¶
References: Geweke & Keane (1997), Mixture of Normals Probit Models (FRB Minneapolis SR 237); Albert & Chib (1993); Bliss (1935), The calculation of the dosage-mortality curve.
| Section | Content |
|---|---|
| Model | Binary probit with a mixture-of-normals disturbance = a flexible (possibly asymmetric) link |
| Section 1 | Synthetic validation — recover an asymmetric error; WAIC favors the mixture |
| Section 2 | Bliss (1935) beetle dose–mortality — the mixture matches the cloglog link the symmetric probit/logit miss |
From-scratch sampler: mixprobit_gibbs.py (mixprobit_full_gibbs).
Model — the probit's error is the link¶
$$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})$$
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 let $F$ be a mixture of normals, which can approximate essentially any shape, so the link is estimated rather than assumed.
- Scale mixture ($\alpha_j=0$): symmetric, fat-tailed (generalizes the robit).
mixprobit_gibbs. - Full mixture (free $\alpha_j$): asymmetric — needed when the dose–response curve is skewed.
mixprobit_full_gibbs(used here).
Identification: a probit fixes the error location and scale; each sweep we recentre the mixture to mean 0 (location $\to$ intercept) and rescale to $\mathrm{Var}(\varepsilon)=1$, then sort components by mean. Choice probability: $P(y=1\mid x)=\sum_j p_j\,\Phi\!\big(\sqrt{h_j}\,(x'\beta+\alpha_j)\big)$.
Gibbs blocks: (1) latent $\tilde y$ — truncated normal (Albert–Chib); (2) component $s_t$ — categorical; (3) $\beta$ — weighted Gaussian; (4) $(\alpha_j,h_j)$ — Normal-Inverse-Gamma; (5) $p$ — Dirichlet. Comparison by WAIC, plus the grouped deviance to benchmark against the classic links.
What is WAIC?¶
WAIC (Widely Applicable Information Criterion; Watanabe 2010) measures a model's out-of-sample predictive accuracy directly from the posterior draws — lower is better. With draws $\theta^{(1)},\dots,\theta^{(S)}$ from the posterior,
$$\mathrm{WAIC} = -2\Big(\underbrace{\textstyle\sum_{i} \log \tfrac1S\textstyle\sum_s p(y_i\mid\theta^{(s)})}_{\text{lppd — fit}} \;-\; \underbrace{\textstyle\sum_i \mathrm{Var}_s\big[\log p(y_i\mid\theta^{(s)})\big]}_{p_\mathrm{WAIC}\ \text{— complexity penalty}}\Big).$$
- lppd (log pointwise predictive density): each observation's likelihood averaged over the posterior, then logged and summed — rewards fit.
- $p_\mathrm{WAIC}$ (effective number of parameters): the posterior variance of each observation's log-likelihood, summed — penalizes models that fit only by being flexible.
Unlike AIC/DIC (which plug in a single point estimate), WAIC uses the whole posterior and asymptotically approximates leave-one-out cross-validation. For a binary model the per-observation likelihood is just the closed-form choice probability $P(y_t{=}1\mid x_t)=\sum_j p_j\,\Phi(\sqrt{h_j}(x_t'\beta+\alpha_j))$ (and $1-P$ for $y_t{=}0$), computed by waic() in mixprobit_gibbs.py. It is the predictive-accuracy analog of Geweke–Keane's Bayes factors. Caveat for binary data: a 0/1 outcome carries little information about error shape, so WAIC moves only modestly here — which is why we also report the grouped deviance (a sharper, aggregated fit measure that lines up with the classic links).
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from mixprobit_gibbs import mixprobit_full_gibbs, waic, prob_curve, error_density_full
print('mixprobit_full_gibbs loaded')
mixprobit_full_gibbs loaded
1. Synthetic validation (asymmetric error)¶
$y=\mathbf{1}[\beta_0+\beta_1 x+\varepsilon>0]$, $n=4000$, true $\beta=(0.3,1.2)$, with a skewed error (80% $N(-0.4,0.5^2)$ + 20% $N(1.6,1^2)$, standardized). A symmetric probit ($m=1$) cannot represent the skew; the full mixture should recover it and improve WAIC.
rng = np.random.default_rng(1)
n = 4000
X = np.column_stack([np.ones(n), rng.normal(size=n)])
beta_true = np.array([0.3, 1.2])
zc = rng.random(n) < 0.8
e = np.where(zc, rng.normal(-0.4, 0.5, n), rng.normal(1.6, 1.0, n)); e = (e - e.mean()) / e.std()
y = ((X @ beta_true + e) > 0).astype(float)
sfits = {m: mixprobit_full_gibbs(y, X, m=m, R=6000, burn=2500, seed=2) for m in (1, 2, 3)}
print(f"{'m':>3}{'WAIC':>10} beta mixture (alpha, p)")
for m, o in sfits.items():
extra = '' if m == 1 else f" alpha={np.round(o['alpha'].mean(0),2)} p={np.round(o['p'].mean(0),2)}"
print(f"{m:>3}{waic(o):>10.1f} {np.round(o['beta'].mean(0),3)}{extra}")
print(f"\ntruth beta = {beta_true}")
m WAIC beta mixture (alpha, p) 1 3192.4 [0.296 1.45 ]
2 3083.4 [0.324 1.329] alpha=[-0.47 1.23] p=[0.71 0.29] 3 3081.1 [0.323 1.299] alpha=[-0.62 0.12 1.94] p=[0.47 0.39 0.14] truth beta = [0.3 1.2]
grid = np.linspace(e.min(), e.max(), 300)
plt.figure(figsize=(8, 4.6))
plt.hist(e, bins=60, density=True, color='lightgray', label='true errors (sample)')
from scipy.stats import norm as _n
plt.plot(grid, _n.pdf(grid, 0, 1), 'b--', lw=2, label='symmetric probit (m=1)')
plt.plot(grid, error_density_full(sfits[2], grid), 'r-', lw=2.2, label='full mixture (m=2)')
plt.xlabel('error'); plt.ylabel('density')
plt.title('Synthetic: full mixture recovers the skew a symmetric probit cannot')
plt.legend(); plt.tight_layout()
plt.savefig('mixprobit_synth.png', dpi=120, bbox_inches='tight'); plt.show()
1. Results — synthetic¶
| m | WAIC | β | mixture |
|---|---|---|---|
| 1 (probit) | 3192.4 | [0.30, 1.45] | — (symmetric) |
| 2 | 3083.4 | [0.32, 1.33] | α=[−0.47, 1.23], p=[0.71, 0.29] |
| 3 | 3081.1 | [0.32, 1.30] | — |
WAIC drops ≈109 from the symmetric probit to the full mixture, which recovers the right-skewed error (a dominant low component + a small high one) and de-biases the slope (1.45 → 1.33 vs truth 1.2). Sampler validated for asymmetric errors.
2. Bliss (1935) beetle dose–mortality¶
The data that birthed the probit: flour beetles exposed to 8 doses of gaseous CS₂, number killed at each. It is the textbook case where symmetric links underfit — the mortality curve is asymmetric, and the complementary-log-log (cloglog) link fits far better than probit or logit. Can a mixture-of-normals probit discover that asymmetry on its own?
import pandas as pd
beetle = pd.read_csv('bliss_beetle.csv') # Bliss (1935) grouped dose-mortality
dose = beetle['dose'].values
ntot = beetle['ntot'].values.astype(int)
kill = beetle['kill'].values.astype(int)
print(f"{ntot.sum()} beetles, {kill.sum()} killed; mortality {np.round(kill/ntot,3)}")
# classic GLM links on grouped data (residual df = 6)
Xg = sm.add_constant(dose); yg = np.column_stack([kill, ntot - kill])
links = [(sm.families.links.Probit(), 'probit'),
(sm.families.links.Logit(), 'logit'),
(sm.families.links.CLogLog(), 'cloglog')]
print("\nclassic link deviances:")
for lk, nm in links:
r = sm.GLM(yg, Xg, family=sm.families.Binomial(link=lk)).fit()
print(f" {nm:<8} deviance = {r.deviance:.2f}")
481 beetles, 291 killed; mortality [0.102 0.217 0.29 0.5 0.825 0.898 0.984 1. ] classic link deviances: probit deviance = 10.12 logit deviance = 11.23 cloglog deviance = 3.45
# expand grouped -> individual binary, fit mixture probit
xi = np.repeat(dose, ntot)
yi = np.concatenate([[1.0] * k + [0.0] * (nt - k) for k, nt in zip(kill, ntot)])
xbar = xi.mean()
Xi = np.column_stack([np.ones(len(xi)), xi - xbar])
Xg2 = np.column_stack([np.ones(8), dose - xbar])
def grouped_deviance(pred):
p = np.clip(pred, 1e-12, 1 - 1e-12); obs = kill / ntot
with np.errstate(divide='ignore', invalid='ignore'): # 0*log0 := 0
t1 = np.where(kill > 0, kill * np.log(obs / p), 0.0)
t2 = np.where(ntot - kill > 0, (ntot - kill) * np.log((1 - obs) / (1 - p)), 0.0)
return 2 * np.sum(t1 + t2)
bfits = {m: mixprobit_full_gibbs(yi, Xi, m=m, R=8000, burn=3000, seed=3) for m in (1, 2, 3)}
print(f"{'m':>3}{'WAIC':>10}{'grouped dev':>14} mixture")
for m, o in bfits.items():
dev = grouped_deviance(prob_curve(o, Xg2))
extra = '' if m == 1 else f" alpha={np.round(o['alpha'].mean(0),2)} p={np.round(o['p'].mean(0),2)}"
print(f"{m:>3}{waic(o):>10.1f}{dev:>14.2f}{extra}")
print("\n(conventional probit deviance 10.12; cloglog 3.45 -- see previous cell)")
m WAIC grouped dev mixture 1 375.3 10.13 2 371.0 3.40 alpha=[-0.46 1.34] p=[0.72 0.28] 3 370.7 3.60 alpha=[-0.68 0.01 1.52] p=[0.45 0.33 0.21] (conventional probit deviance 10.12; cloglog 3.45 -- see previous cell)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.8))
dg = np.linspace(dose.min() - 0.01, dose.max() + 0.01, 200)
Xgrid = np.column_stack([np.ones(len(dg)), dg - xbar])
probit_fit = sm.GLM(yg, Xg, family=sm.families.Binomial(link=sm.families.links.Probit())).fit()
cll_fit = sm.GLM(yg, Xg, family=sm.families.Binomial(link=sm.families.links.CLogLog())).fit()
ax[0].plot(dose, kill / ntot, 'ko', ms=7, label='observed mortality')
ax[0].plot(dg, probit_fit.predict(sm.add_constant(dg)), 'b--', lw=2, label='conventional probit (dev 10.1)')
ax[0].plot(dg, cll_fit.predict(sm.add_constant(dg)), 'g:', lw=2, label='cloglog (dev 3.45)')
ax[0].plot(dg, prob_curve(bfits[2], Xgrid), 'r-', lw=2.2, label='mixture probit m=2 (dev 3.4)')
ax[0].set_xlabel('log dose (CS2)'); ax[0].set_ylabel('P(killed)')
ax[0].set_title('Beetle dose-response: mixture matches cloglog, beats probit'); ax[0].legend()
gx = np.linspace(-3, 4, 300)
from scipy.stats import norm as _n
ax[1].plot(gx, _n.pdf(gx, 0, 1), 'b--', lw=2, label='normal (probit error)')
ax[1].plot(gx, error_density_full(bfits[2], gx), 'r-', lw=2.2, label='estimated mixture error')
ax[1].axvline(0, color='gray', lw=.6)
ax[1].set_xlabel('disturbance'); ax[1].set_ylabel('density')
ax[1].set_title('Estimated error is right-skewed (not the symmetric normal)'); ax[1].legend()
plt.tight_layout()
plt.savefig('mixprobit_beetle.png', dpi=120, bbox_inches='tight'); plt.show()
2. Results — beetle¶
| link / model | deviance |
|---|---|
| logit | 11.23 |
| conventional probit | 10.12 |
| cloglog | 3.45 |
| mixture probit, m=2 | 3.40 |
| mixture probit, m=3 | 3.60 |
The flexible-error probit reproduces the asymmetry by itself. A 2-component mixture cuts the grouped deviance from 10.1 (conventional probit) to 3.40 — essentially equal to the cloglog link (3.45), the link statisticians historically had to choose by hand for this data. (m=1 reproduces the textbook probit deviance, confirming the comparison; WAIC on the 481 individual outcomes moves the same way, 375→371, more modestly because a binary outcome carries limited information about error shape.)
The estimated disturbance is right-skewed (α = [−0.46, +1.34], p = [0.72, 0.28]) — a dominant lower component plus a smaller upper one — which is exactly the skew the symmetric normal (and logistic) cannot produce, and why they misfit the tails of the mortality curve.
Takeaway: in a binary model the error distribution is the link. Geweke–Keane's mixture turns link choice into estimation: instead of testing probit vs. logit vs. cloglog, you let the data shape the disturbance — and here it independently recovers the cloglog-quality fit. (Contrast: on the PSID 90/10 participation data the binary signal was too weak to detect any non-normality; the beetle data, with mortality sweeping 0→1 across the dose range, is far more informative about the link shape.)