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$.

Which means there are three different numbers in play whenever the two models are compared, and it is worth naming them once. The structural $\beta$ lives on the $t_\nu$ scale and is what generated the data; the probit-equivalent $\beta$ is that same quantity shrunk by $\sqrt{(\nu-2)/\nu}$, and is the most a correctly-working probit can hope to recover; and the corrected robit estimate is the raw robit draw multiplied by that factor. The pairing to remember is that the raw robit estimate should match the structural truth, while the corrected robit estimate should match what standard probit returns — the correction exists to put robit and probit on a common footing, not to recover the structural value. On the synthetic data below, the four numbers for the main slope are: structural truth $1.200$, probit-equivalent truth $0.930$, standard probit $1.049$, robit raw $1.277$, robit corrected $0.989$.

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%; the correctly-specified robit ($\nu=5$) recovers the structural $\beta$ with near-zero bias — mean $|$bias$|$ roughly halved, 0.15 → 0.07, both measured against the $t_5$-scale truth and using the raw robit estimate, per the three-number convention above. Applying the scale correction here would be a mistake: it would return $0.989$, which is the right answer to a different question (what would a probit see?) and looks like a 18% failure against the structural 1.2. 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 as $\nu$ rises — so $\nu$ matters for parameter recovery even though it is invisible to model-selection by LPML. One caveat on reading that curve: the raw slope is on the $t_\nu$ scale, and that scale itself depends on $\nu$, so the fixed structural value 1.2 is the correct target only at the true $\nu=5$ — where the raw slope is already $1.277$, slightly above it. The curve crosses 1.2 further along, near $\nu=7$, and that crossing is an artefact of comparing a moving scale against a fixed number rather than an indication that $\nu=7$ is preferred. 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