""" binary_probit_gibbs.py Bayesian binary probit model via Albert & Chib (1993) data augmentation. Model ----- y_i = 1{z_i > 0}, z_i = x_i'β + ε_i, ε_i ~ N(0, 1) Pr(y_i = 1 | x_i) = Φ(x_i'β) The error variance is fixed at 1 for identification (probit normalisation). Gibbs sampler — two blocks (Albert & Chib 1993) ------------------------------------------------ (1) z_i | y_i, β ~ TN(x_i'β, 1; 0, ∞) if y_i = 1 ~ TN(x_i'β, 1; -∞, 0) if y_i = 0 (2) β | z ~ N(b̄, B̄) B̄⁻¹ = X'X + B₀⁻¹ b̄ = B̄ (X'z + B₀⁻¹ b₀) Default prior: β ~ N(0, 100·I) (weakly informative) Bayesian residuals for outlier detection (Albert & Chib 1995) -------------------------------------------------------------- The data-augmentation representation makes latent residuals directly available as posterior draws: e_i = z_i - x_i'b posterior latent residual Under the model: e_i ~ N(0,1) truncated to the sign of y_i. Outliers show large |e_bar_i|: z_i is squeezed against the boundary. p_bar_i = E[Phi(x_i'b) | y] posterior predictive probability CPO_i = p(y_i | y_{-i}) conditional predictive ordinate ~ {(1/R) sum_r 1/p(y_i|b^r)}^{-1} (harmonic mean estimator) Interpretation: Low CPO_i -> observation i is surprising given the rest of the data. For y_i=1 with x_i'b_bar << 0: z_i squeezed just above 0 -> e_bar_i >> 0. For y_i=0 with x_i'b_bar >> 0: z_i squeezed just below 0 -> e_bar_i << 0. References ---------- Albert, J. H., & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association, 88(422), 669-679. [Gibbs sampler via data augmentation] Albert, J. H., & Chib, S. (1995). Bayesian residual analysis for binary response regression models. Biometrika, 82(4), 747-759. [Latent residual diagnostics and outlier detection] """ import time import numpy as np import pandas as pd from scipy.special import ndtr, ndtri # ── public API ──────────────────────────────────────────────────────────────── def binary_probit_gibbs( y, X, b0=None, B0=None, R=10000, burn=2000, seed=42, nprint=1000, ): """ Albert & Chib (1993) Gibbs sampler for the binary probit model. Parameters ---------- y : (n,) binary outcomes (0/1) X : (n, k) design matrix (include intercept column if desired) b0 : (k,) prior mean for β [default: zeros] B0 : (k, k) prior covariance for β [default: 100·I] R, burn : int total and burn-in iterations seed : int nprint : int print progress every nprint iterations (0 = silent) Returns ------- dict with keys 'beta' : (R_kept, k) posterior draws for β 'z' : (R_kept, n) posterior draws for latent utilities z_i """ rng = np.random.default_rng(seed) y = np.asarray(y, float).ravel() X = np.asarray(X, float) n, k = X.shape y_bool = y.astype(bool) # ── prior ───────────────────────────────────────────────────────────────── 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)) # ── precompute (σ²=1 is fixed, so Bbar is constant across iterations) ──── XtX = X.T @ X Prec = XtX + B0_inv # posterior precision (constant) Bbar = np.linalg.inv(Prec) # posterior covariance (constant) L_chol = np.linalg.cholesky(Bbar) B0_inv_b0 = B0_inv @ b0 # ── initial values ──────────────────────────────────────────────────────── beta = np.linalg.lstsq(X, y - 0.5, rcond=None)[0] z = np.where(y_bool, np.abs(rng.standard_normal(n)), -np.abs(rng.standard_normal(n))) # ── storage ─────────────────────────────────────────────────────────────── 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 | y_i, β (truncated Normal draw) ─────────────────── mu = X @ beta # (n,) u = rng.uniform(size=n) Phi_neg = ndtr(-mu) # Φ(−μᵢ) = P(zᵢ < 0 | μᵢ) gap_pos = 1.0 - Phi_neg # Φ(μᵢ) = P(zᵢ > 0) — mass for y=1 gap_neg = Phi_neg # Φ(−μᵢ) = P(zᵢ < 0) — mass for y=0 # y=1: draw from TN(μᵢ, 1; 0, ∞) # p = Φ(−μᵢ) + u·Φ(μᵢ) in (Φ(−μᵢ), 1) p1 = Phi_neg + u * gap_pos # Extreme case (μᵢ ≪ −37): gap_pos underflows → use Exp(|μ|) approx # E[z | z>0, N(μ,1)] ≈ 1/|μ| when μ→−∞ → z ~ Exp(|μ|) z1 = np.where( gap_pos < 1e-300, -np.log(np.clip(u, 1e-15, 1.0)) / np.maximum(-mu, 1e-6), ndtri(np.clip(p1, 1e-15, 1.0 - 1e-15)) + mu, ) # y=0: draw from TN(μᵢ, 1; −∞, 0) # p = u·Φ(−μᵢ) in (0, Φ(−μᵢ)) p0 = u * gap_neg # Extreme case (μᵢ ≫ 37): gap_neg underflows → z ~ −Exp(μ) z0 = np.where( gap_neg < 1e-300, np.log(np.clip(u, 1e-15, 1.0)) / np.maximum(mu, 1e-6), ndtri(np.clip(p0, 1e-15, 1.0 - 1e-15)) + mu, ) z = np.where(y_bool, z1, z0) # ── Block 2: β | z (Normal posterior — Bbar is precomputed) ───────── bbar = Bbar @ (X.T @ z + B0_inv_b0) beta = bbar + L_chol @ rng.standard_normal(k) # ── store ───────────────────────────────────────────────────────────── if r > burn: beta_draws[kept] = beta z_draws[kept] = z kept += 1 if nprint > 0 and r % nprint == 0: print(f' Iter {r:5d}/{R} ({time.time()-t0:.1f}s elapsed)') elapsed = time.time() - t0 print(f'Done in {elapsed/60:.1f} min | kept draws: {R_kept}') return {'beta': beta_draws, 'z': z_draws} def probit_mle(y, X): """ Maximum likelihood estimation for the binary probit model via statsmodels. Returns ------- dict with keys 'beta' : (k,) MLE parameter estimates 'se' : (k,) asymptotic standard errors 'vcov' : (k, k) estimated covariance matrix 'loglik' : float log-likelihood at the MLE 'result' : statsmodels ProbitResults object """ import statsmodels.api as sm res = sm.Probit(y, X).fit(disp=False) return { 'beta': res.params, 'se': res.bse, 'vcov': res.cov_params(), 'loglik': res.llf, 'result': res, } def bayesian_residuals(y, X, beta_draws, z_draws): """ Compute Bayesian residuals and CPO-based outlier diagnostics. Parameters ---------- y : (n,) binary outcomes X : (n, k) design matrix beta_draws : (R, k) posterior β draws z_draws : (R, n) posterior latent-z draws Returns ------- pandas.DataFrame with columns: y observed outcome (0/1) p_bar posterior predictive probability E[Φ(x'β) | y] CPO conditional predictive ordinate (harmonic mean estimator) log_CPO log CPO (contribution to LPML) e_bar posterior mean latent residual E[z − x'β | y] e_sd posterior SD of latent residual surprise |y − p_bar| (classical-style misfit; 1 = complete surprise) """ y = np.asarray(y, float).ravel() n = len(y) # Posterior predictive probabilities p̄ᵢ = E[Φ(xᵢ'β) | y] mu_draws = X @ beta_draws.T # (n, R) p_mat = ndtr(mu_draws) # (n, R) Φ(xᵢ'β^(r)) p_bar = p_mat.mean(axis=1) # (n,) # CPO: harmonic mean of p(yᵢ | β^(r)) lik = np.where(y[:, None], p_mat, 1.0 - p_mat) # (n, R) lik = np.clip(lik, 1e-10, 1.0) cpo = 1.0 / np.mean(1.0 / lik, axis=1) # (n,) # Latent residuals ēᵢ = E[zᵢ − xᵢ'β | y] e_draws = z_draws - mu_draws.T # (R, n) e_bar = e_draws.mean(axis=0) # (n,) e_sd = e_draws.std(axis=0) # (n,) return pd.DataFrame({ 'y': y.astype(int), 'p_bar': p_bar.round(4), 'CPO': cpo.round(6), 'log_CPO': np.log(np.clip(cpo, 1e-15, 1)).round(4), 'e_bar': e_bar.round(4), 'e_sd': e_sd.round(4), 'surprise': np.abs(y - p_bar).round(4), }) def posterior_summary(draws, param_names=None, true_vals=None): """ Posterior mean, SD, and 95% credible interval for each column. Parameters ---------- draws : (R, p) posterior draws param_names : list of p strings [optional] true_vals : (p,) true values; adds 'true' and 'covered' columns Returns ------- pandas.DataFrame """ p = draws.shape[1] if param_names is None: param_names = [f'b{j}' for j in range(p)] lo = np.quantile(draws, 0.025, axis=0) hi = np.quantile(draws, 0.975, axis=0) df = pd.DataFrame({ 'mean': draws.mean(0), 'sd': draws.std(0), 'q025': lo, 'q975': hi, }, index=param_names) if true_vals is not None: tv = np.asarray(true_vals) df['true'] = tv df['covered'] = (tv >= lo) & (tv <= hi) return df.round(4)