Bayesian Robust Binary Probit

Model

Robit regression (Liu 2004) replaces the standard binary probit's normal link $\Phi$ with a Student-$t_\nu$ CDF, making the model robust to label noise and heavy-tailed errors. The key representation (Geweke 1993): each observation gets a latent precision weight $\lambda_i \sim \text{Gamma}(\nu/2, \nu/2)$ that scales its variance. Outliers acquire small $\lambda_i$ — large local variance — and are automatically downweighted in the estimation of $\beta$. As $\nu \to \infty$, $T_\nu \to \Phi$ and the model reduces to standard binary probit.

• $\beta$ — coefficient vector; prior $N(b_0, B_0)$
• $\lambda_i$ — per-observation precision weight; $E[\lambda_i] = 1$, small value = outlier
• $\nu$ — degrees of freedom; $\nu \in [3,7]$ heavy tails; $\nu \to \infty$ collapses to standard probit
• $T_\nu$ — standard Student-$t$ CDF with $\nu$ d.f.

3-Block Gibbs — Liu (2004) / Geweke (1993)

The 3-block Gibbs sampler extends Albert & Chib (1993) with a third block for the precision weights. The critical difference from standard probit: in the robit Block 2, the posterior precision $\bar{B}^{-1} = X'\Lambda X + B_0^{-1}$ depends on $\Lambda = \text{diag}(\lambda_i)$ which changes every iteration, so $\bar{B}$ must be recomputed at each step at cost $O(nk^2)$ (vs. the one-time Cholesky of standard probit).

Scale Correction

Robit estimates $\beta$ in the $t_\nu$ scale, where the marginal variance of $z_i^*$ is $\nu/(\nu-2)$ rather than 1. To compare with standard probit estimates multiply by the scale correction $\sqrt{(\nu-2)/\nu}$ — this maps every robit estimate back to the probit-equivalent value regardless of $\nu$.

Results

Applied to three scenarios. Section 1 — Leptokurtic synthetic data ($n=400$, true $t_5$ errors): standard probit converges to $\beta_{\text{true}} \times \sqrt{3/5}$ — systematically attenuated by 23%; correctly-specified robit ($\nu=5$) recovers the structural $\beta$ with near-zero bias (mean $|$bias$|$ roughly halved, 0.15 → 0.07). The predictive gap, by contrast, is modest — LPML (log pseudo-marginal likelihood, the leave-one-out predictive score $\sum_i \log \text{CPO}_i$; higher is better) favours robit by only $\approx +0.5$ nat and the CPO-by-residual comparison is a wash: robit's Section 1 advantage is in parameter recovery, not label prediction, because a binary outcome barely encodes the latent tail shape. The precision weights $\bar{\lambda}_i$ are negatively correlated ($r = -0.45$) with true residual magnitude. Section 2 — Sensitivity to $\nu$: on this $n=400$ sample the LPML curve is nearly flat (total spread $\approx 0.5$ nat) — a binary outcome barely encodes the tail index, so LPML cannot identify $\nu$ (its grid maximum sits at the edge $\nu=3$, not the true 5); the one robust finding is that every robit $\nu$ beats standard probit. The raw slope, by contrast, moves materially — declining monotonically with $\nu$ and crossing the true 1.2 near $\nu=7$ — so $\nu$ matters for parameter recovery even though it is invisible to model-selection by LPML. Section 3 — Finney (1947) vasoconstriction ($n=39$, $\nu=7$): LPML is essentially flat ($\Delta = -0.02$ nats) at this sample size, but $\bar{\lambda}_i$ correctly identifies the most anomalous observations — including obs 4 ($\bar{\lambda} = 0.667$, not flagged by Liu) and Liu's discordant obs 18 ($\bar{\lambda} = 0.725$) — without any analyst input.

Notebook

Downloads

Gibbs Sampler — Source Code

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

References

• Liu, C. (2004). Robit regression: a simple robust alternative to logistic and probit regression. In A. Gelman & X.-L. Meng (eds.), Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives, 227–238. Wiley. — coined robit; the Student-$t$ link and the $\lambda_i$ diagnostic implemented here
• Geweke, J. (1993). Bayesian treatment of the independent Student-t linear model. Journal of Applied Econometrics 8(S1), S19–S40. — the Normal scale-mixture representation the Gibbs sampler is built on
• 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 probit that Block 1–2 extend
• Holmes, C. C. & Held, L. (2006). Bayesian auxiliary variable models for binary and multinomial regression. Bayesian Analysis 1(1), 145–168. — auxiliary-variable Gibbs for binary/multinomial models with scale mixtures
• Finney, D. J. (1947). Probit Analysis. Cambridge University Press. — Table 1 vasoconstriction data (Beall's 1929 experiment) used in Section 3