Bayesian Binary Probit

Model

The binary probit model links a latent continuous utility $z_i$ to a binary outcome $y_i$ via a threshold rule. The error variance is fixed at 1 for identification (the probit normalisation): only the ratio $\beta/\sigma$ is identified, so $\sigma=1$ is imposed.

• $\beta$ — coefficient vector; $N(b_0, B_0)$ prior (default: $N(0, 100I)$)
• $z_i$ — latent utility; drawn each Gibbs iteration as a truncated Normal
• $\Phi$ — standard Normal CDF

2-Block Gibbs — Albert & Chib (1993)

Albert & Chib (1993) introduce a data-augmentation trick: conditioning on the latent utilities $z_i$ reduces the probit likelihood to a standard Normal regression, making $\beta$ conjugate. The posterior precision $\bar{B}^{-1} = X'X + B_0^{-1}$ is constant across iterations ($\sigma^2=1$ fixed), so the Cholesky factorisation is computed once — making the sampler very fast regardless of $n$.

Bayesian Residuals — Albert & Chib (1995)

Because the Gibbs sampler draws the latent utilities $z_i^{(r)}$ at every iteration, latent residuals $\bar{e}_i = E[z_i - x_i'\beta \mid y]$ are directly available as posterior samples — a diagnostic with no classical analogue. Under the correctly specified model, $e_i \sim N(0,1)$ truncated to the sign of $y_i$. Outliers (e.g. label-flipped observations) are squeezed against the boundary, producing large $|\bar{e}_i|$. The CPO complements this with a proper probabilistic measure: $\text{CPO}_i = p(y_i \mid y_{-i})$ — how probable is $y_i$ given the model trained on all other observations.

Notebooks

Applied to synthetic binary data ($n=400$, $k=4$, 25 label-flipped outliers). Section 1 compares Gibbs posterior with MLE — both recover similar coefficient estimates despite label noise. Section 2 demonstrates Bayesian residual diagnostics: CPO and latent residuals $\bar{e}_i$ recover 13/25 planted outliers at $k=25$ (matching Pearson at equal flag rate). Section 3 shows where CPO strictly outperforms Pearson: with $n=30$ and 5 influential outliers at extreme $x$, an informative prior keeps Gibbs draws near the true coefficient while the MLE collapses — CPO detects 4/5 outliers; Pearson detects 0/5.

Downloads

Gibbs Sampler — Source Code

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

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. — the data-augmentation Gibbs sampler
• Albert, J. H. & Chib, S. (1995). Bayesian residual analysis for binary response regression models. Biometrika 82(4), 747–759. — CPO and latent residuals
• Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. — bayesm::rbprobitGibbs