""" robust_probit_gibbs.py Bayesian robust binary probit ("robit regression") via data-augmentation Gibbs sampling. Errors follow a Student-t distribution implemented via a Normal scale-mixture (Geweke 1993). Robit model (Liu 2004) ----------------------- y_i = 1(z_i* > 0) z_i* | λ_i ~ N(x_i'β, 1/λ_i) latent utility, t-distributed noise λ_i ~ Gamma(ν/2, ν/2) precision weight (E[λ_i] = 1) β ~ N(b₀, B₀) prior Marginal over λ_i: z_i* ~ t_ν(x_i'β, 1) [Student-t link] P(y_i=1 | β) = T_ν(x_i'β) t-CDF at linear predictor As ν → ∞, T_ν → Φ and the model reduces to standard binary probit. Small ν (3–7) yields heavy tails robust to label noise / outliers. Outliers get small λ_i → large variance 1/λ_i → little influence on β. Gibbs sampler — 3 blocks ------------------------- Block 1. z_i | β, λ_i, y_i ~ TN(x_i'β, 1/λ_i ; lo_i, hi_i) lo_i = 0, hi_i = +∞ if y_i = 1 lo_i = -∞, hi_i = 0 if y_i = 0 Block 2. β | z, Λ ~ N(b̄, B̄) B̄⁻¹ = X'ΛX + B₀⁻¹, b̄ = B̄(X'Λz + B₀⁻¹b₀) (B̄ changes every iteration because Λ = diag(λ) changes) Block 3. λ_i | z_i, β ~ Gamma((ν+1)/2, (ν + (z_i − x_i'β)²)/2) rate parameterisation; large |residual| → small λ_i Package landscape ----------------- No standard Python or R package ships a ready-made Bayesian t-link binary probit. This script fills the gap. Closest alternatives: R: rstan / brms (custom Stan model), LaplacesDemon (user-coded) robustbase::glmrob (frequentist M-estimation, different mechanism) Python: PyMC (scale-mixture coded manually), pystan References ---------- Liu, C. (2004). Robit Regression: A Simple Robust Alternative to Logistic and Probit Regression. In Gelman & Meng (Eds.), Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives, pp. 227-238. Wiley. Geweke, J. (1993). Bayesian Treatment of the Independent Student-t Linear Model. Journal of Applied Econometrics 8(S1), S19-S40. Albert, J. H., & Chib, S. (1993). Bayesian Analysis of Binary and Polychotomous Response Data. JASA 88(422), 669-679. Holmes, C. C., & Held, L. (2006). Bayesian Auxiliary Variable Models for Binary and Multinomial Regression. Bayesian Analysis 1(1), 145-168. """ import time import numpy as np from scipy.special import ndtr, ndtri from scipy.stats import t as t_dist import pandas as pd # ── Truncated Normal helpers ────────────────────────────────────────────────── def _rtn(lo, hi, mu, rng): """ Vectorised draw from N(mu_i, 1) truncated to (lo_i, hi_i). Inverse-CDF; midpoint fallback when tail mass < 1e-300. """ p_lo = ndtr(lo - mu) p_hi = ndtr(hi - mu) gap = p_hi - p_lo u = rng.uniform(size=len(mu)) p = p_lo + u * gap tiny = gap < 1e-300 p = np.where(tiny, 0.5 * (p_lo + p_hi), p) return ndtri(np.clip(p, 1e-15, 1.0 - 1e-15)) + mu def _rtn_scaled(lo, hi, mu, sd, rng): """ Vectorised draw from N(mu_i, sd_i²) truncated to (lo_i, hi_i). Per-observation standard deviation sd — required for robit Block 1. """ a = (lo - mu) / sd b = (hi - mu) / sd p_lo = ndtr(a) p_hi = ndtr(b) gap = p_hi - p_lo u = rng.uniform(size=len(mu)) p = p_lo + u * gap tiny = gap < 1e-300 p = np.where(tiny, 0.5 * (p_lo + p_hi), p) return ndtri(np.clip(p, 1e-15, 1.0 - 1e-15)) * sd + mu # ── Standard binary probit (Albert & Chib 1993) ─────────────────────────────── def probit_gibbs( y, X, b0=None, B0=None, R=10000, burn=2000, seed=42, nprint=1000, ): """ Albert & Chib (1993) 2-block Gibbs for standard binary probit. Equivalent to robit_gibbs with ν → ∞ (all λ_i fixed at 1). Parameters ---------- y : (n,) int binary outcomes in {0, 1} X : (n, k) design matrix (include intercept column if desired) b0 : (k,) prior mean [default zeros] B0 : (k, k) prior covariance [default 100·I] R, burn, seed, nprint : int Returns ------- dict 'beta' (R_kept, k) | 'z' (R_kept, n) """ rng = np.random.default_rng(seed) y = np.asarray(y, int).ravel() X = np.asarray(X, float) n, k = X.shape if b0 is None: b0 = np.zeros(k) if B0 is None: B0 = 100.0 * np.eye(k) b0 = np.asarray(b0, float) B0_inv = np.linalg.inv(np.asarray(B0, float)) B0_inv_b0 = B0_inv @ b0 # Posterior covariance is constant (Λ = I throughout) Bbar_inv = X.T @ X + B0_inv Bbar = np.linalg.inv(Bbar_inv) L = np.linalg.cholesky(Bbar) NINF, PINF = -1e10, 1e10 lo = np.where(y == 1, 0.0, NINF) hi = np.where(y == 1, PINF, 0.0) beta = np.zeros(k) z = np.where(y == 1, 0.5, -0.5) R_kept = R - burn beta_draws = np.empty((R_kept, k)) z_draws = np.empty((R_kept, n)) kept = 0 t0 = time.time() for r in range(1, R + 1): # Block 1: z_i | beta, y_i ~ TN(x_i'beta, 1) mu = X @ beta z = _rtn(lo, hi, mu, rng) # Block 2: beta | z ~ N(bbar, Bbar) [Bbar constant] bbar = Bbar @ (X.T @ z + B0_inv_b0) beta = bbar + L @ rng.standard_normal(k) if r > burn: beta_draws[kept] = beta z_draws[kept] = z kept += 1 if nprint > 0 and r % nprint == 0: print(f' [probit std] iter {r:6d}/{R} ({time.time()-t0:.1f}s)') print(f'[probit std] Done in {time.time()-t0:.1f}s | kept: {R_kept}') return {'beta': beta_draws, 'z': z_draws} # ── Robit regression (Liu 2004, Geweke 1993) ────────────────────────────────── def robit_gibbs( y, X, nu=7.0, b0=None, B0=None, R=10000, burn=2000, seed=42, nprint=1000, ): """ Robust binary probit ("robit") via data-augmentation Gibbs sampling. Errors are Student-t with ν degrees of freedom via Normal scale-mixture. Parameters ---------- y : (n,) int binary outcomes in {0, 1} X : (n, k) design matrix nu : float > 2 degrees of freedom 3–7 heavy tails, robust to outliers 15+ near-normal; 1e6+ collapses to probit_gibbs b0 : (k,) prior mean [default zeros] B0 : (k, k) prior covariance [default 100·I] R, burn, seed, nprint : int Returns ------- dict 'beta' : (R_kept, k) coefficient draws 'lambda' : (R_kept, n) precision weights (small = influential outlier) 'z' : (R_kept, n) latent utility draws """ rng = np.random.default_rng(seed) y = np.asarray(y, int).ravel() X = np.asarray(X, float) n, k = X.shape if b0 is None: b0 = np.zeros(k) if B0 is None: B0 = 100.0 * np.eye(k) b0 = np.asarray(b0, float) B0_inv = np.linalg.inv(np.asarray(B0, float)) B0_inv_b0 = B0_inv @ b0 NINF, PINF = -1e10, 1e10 lo = np.where(y == 1, 0.0, NINF) hi = np.where(y == 1, PINF, 0.0) beta = np.zeros(k) lam = np.ones(n) # precision weights; init at 1 (probit) z = np.where(y == 1, 0.5, -0.5) shape_lam = (nu + 1.0) / 2.0 # Gamma shape for Block 3 (constant) R_kept = R - burn beta_draws = np.empty((R_kept, k)) lam_draws = np.empty((R_kept, n)) z_draws = np.empty((R_kept, n)) kept = 0 t0 = time.time() for r in range(1, R + 1): # Block 1: z_i | beta, lambda_i, y_i ~ TN(x_i'beta, 1/lambda_i) mu = X @ beta sd = 1.0 / np.sqrt(lam) z = _rtn_scaled(lo, hi, mu, sd, rng) # Block 2: beta | z, Lambda ~ N(bbar, Bbar) # Bbar^{-1} = X'ΛX + B0^{-1}, bbar = Bbar(X'Λz + B0^{-1}b0) XL = X * lam[:, None] # (n, k): row-i scaled by lam_i Bbar_inv = XL.T @ X + B0_inv # (k, k) Bbar = np.linalg.inv(Bbar_inv) L = np.linalg.cholesky(Bbar) bbar = Bbar @ (X.T @ (lam * z) + B0_inv_b0) beta = bbar + L @ rng.standard_normal(k) # Block 3: lambda_i | z_i, beta ~ Gamma((nu+1)/2, (nu + resid_i^2) / 2) # numpy Gamma uses scale = 1/rate resid = z - X @ beta rate_lam = (nu + resid ** 2) / 2.0 lam = rng.gamma(shape_lam, 1.0 / rate_lam) if r > burn: beta_draws[kept] = beta lam_draws[kept] = lam z_draws[kept] = z kept += 1 if nprint > 0 and r % nprint == 0: print(f' [robit ν={nu:.0f}] iter {r:6d}/{R} ({time.time()-t0:.1f}s)') print(f'[robit ν={nu:.0f}] Done in {time.time()-t0:.1f}s | kept: {R_kept}') return {'beta': beta_draws, 'lambda': lam_draws, 'z': z_draws} # ── CPO / LPML ──────────────────────────────────────────────────────────────── def compute_cpo_probit(y, X, beta_draws): """ LPML / CPO for standard probit using normal link Φ. CPO_i = 1 / mean_r(1 / p(y_i | β^r)) LPML = Σ_i log CPO_i Returns ------- lpml : float cpo : (n,) array """ y = np.asarray(y, int) xb = beta_draws @ X.T # (R, n) s = np.where(y == 1, 1.0, -1.0) lik = ndtr(s[None, :] * xb) # (R, n) lik = np.clip(lik, 1e-10, 1.0) cpo = 1.0 / np.mean(1.0 / lik, axis=0) lpml = float(np.sum(np.log(np.clip(cpo, 1e-15, 1.0)))) return lpml, cpo def compute_cpo_robit(y, X, beta_draws, nu): """ LPML / CPO for robit using marginal t-link T_ν. Integrates out λ_i analytically via the Student-t CDF. p(y_i=1 | β^r) = T_ν(x_i'β^r) [standard t CDF with ν d.f.] Using the marginal rather than the conditional CPO avoids conditioning on the observation-specific nuisance λ_i that would be held out with y_i. Returns ------- lpml : float cpo : (n,) array """ y = np.asarray(y, int) xb = beta_draws @ X.T # (R, n) s = np.where(y == 1, 1.0, -1.0) lik = t_dist.cdf(s[None, :] * xb, df=nu) # (R, n) lik = np.clip(lik, 1e-10, 1.0) cpo = 1.0 / np.mean(1.0 / lik, axis=0) lpml = float(np.sum(np.log(np.clip(cpo, 1e-15, 1.0)))) return lpml, cpo # ── Predicted probabilities ─────────────────────────────────────────────────── def pred_probs_probit(X, beta_draws): """Posterior-mean P(y=1 | x_i) using standard normal link Φ.""" xb = beta_draws @ X.T # (R, n) return ndtr(xb).mean(axis=0) # (n,) def pred_probs_robit(X, beta_draws, nu): """Posterior-mean P(y=1 | x_i) using t-link T_ν.""" xb = beta_draws @ X.T # (R, n) return t_dist.cdf(xb, df=nu).mean(axis=0) # (n,) # ── Posterior summary ───────────────────────────────────────────────────────── def posterior_summary(draws, param_names=None): """(R, p) draws → DataFrame with mean, sd, q025, q975.""" p = draws.shape[1] if param_names is None: param_names = [f'p{j}' for j in range(p)] return pd.DataFrame({ 'mean': draws.mean(0).round(4), 'sd': draws.std(0).round(4), 'q025': np.quantile(draws, 0.025, axis=0).round(4), 'q975': np.quantile(draws, 0.975, axis=0).round(4), }, index=param_names)