Bayesian Binary Logit

Model

The binary logit model links the linear index $\eta_i = x_i'\beta$ to a binary outcome $y_i$ through the logistic CDF $\sigma$ . Unlike the probit, the logit posterior has no conjugate form under a Normal prior, so it cannot be sampled with a plain Gibbs step. This example implements two from-scratch routes — a tuned random-walk Metropolis and a Pólya–Gamma data-augmentation Gibbs sampler — and cross-checks both against bambi, PyMC, R's bayesm, and the MLE.

y_i \sim \text{Bernoulli}(p_i), \quad p_i = \sigma(x_i'\beta) = \frac{1}{1+e^{-x_i'\beta}}, \quad \beta \sim N(0, \sigma_\beta^2 I)

$\beta$ — coefficient vector; diffuse Normal prior (default $\sigma_\beta = 10$ )
$\sigma$ — logistic CDF (the inverse-logit link)
$\eta_i = x_i'\beta$ — linear index; $\text{softplus}(\eta) = \log(1+e^\eta)$

Random-Walk Metropolis

Random-walk Metropolis proposes $\beta' = \beta + N(0, \Sigma_{\text{prop}})$ and accepts with probability $\min\{1, \exp[\ell(\beta') - \ell(\beta)]\}$ . The proposal $\Sigma_{\text{prop}} = c^2\hat{\Sigma}$ uses the inverse Fisher information at the MLE $\hat{\Sigma} = (X'\hat{W}X)^{-1}$ so the proposal matches the posterior's local shape, with the Roberts–Gelman–Gilks optimal scale $c = 2.38/\sqrt{k}$ (acceptance $\approx 0.234$ in high dimension). It works, but mixes slowly — lag-1 autocorrelation $\approx 0.86$ .

\ell(\beta) = \sum_i \big[y_i\,\eta_i - \text{softplus}(\eta_i)\big] - \tfrac{1}{2\sigma_\beta^2}\beta'\beta

Pólya–Gamma Gibbs — the logit Albert & Chib

Pólya–Gamma augmentation (Polson, Scott & Windle 2013) makes the logit likelihood conditionally Gaussian — the exact logit analog of Albert–Chib for probit. Introducing a latent $\omega_i \sim \mathrm{PG}(1,\, x_i'\beta)$ per observation renders $\beta \mid \omega$ a Normal full conditional, giving a tuning-free Gibbs sampler with no accept/reject. The rpg draw uses the truncated sum-of-Gammas series representation. Lag-1 autocorrelation drops to $\approx 0.28\text{–}0.51$ . This is the probit/logit symmetry made concrete: probit augments with truncated Normals (Albert–Chib), logit with Pólya–Gamma variables — which is exactly why bayesm does probit by Gibbs but logit by Metropolis.

Block	Draw	Distribution
(1)	$\omega_i \mid \beta$	$\mathrm{PG}(1,\, x_i'\beta)$ — Pólya–Gamma augmentation, one per observation
(2)	$\beta \mid \omega$	Normal conjugate $N(m, V)$ with $V = (X'\Omega X + B_0^{-1})^{-1}$ , $m = V(X'\kappa + B_0^{-1}b_0)$ , $\kappa = y - \tfrac12$

Samplers compared

Method	Type	Mechanism	Tuning
MLE (IRLS)	frequentist	Newton–Raphson	—
RW-Metropolis	Bayesian, scratch	direct MH, inverse-Fisher proposal	proposal scale
Pólya–Gamma Gibbs	Bayesian, scratch	PG augmentation → conjugate	none
`bambi` / PyMC	Bayesian, package	NUTS (HMC)	auto (warmup)
bayesm (R)	Bayesian, package	independence Metropolis, MNL(2)	none

Notebooks

Section 1 validates both from-scratch samplers on synthetic data ( $n=2000$ , $k=4$ ): posterior mean $\approx$ MLE $\approx$ truth to within 0.017 — the Bernstein–von Mises result under a diffuse prior. Section 2 adds the Pólya–Gamma Gibbs and shows the same posterior with far better mixing. Section 3 cross-checks against bambi and PyMC (NUTS) — five routes agree to Monte-Carlo error. Section 4 applies all three from-scratch routes to Mroz (1987) women's labor-force participation ( $n=753$ ; Wooldridge Example 17.1): each child under 6 multiplies the odds of participation by $e^{-1.44}\approx 0.23$ ; education raises it ( $\times 1.25$ /yr); school-age children have no credible effect. The R bayesm notebook reproduces the same results via independence Metropolis — two languages, six samplers, all in agreement.

Python From-scratch RW-Metropolis and Pólya–Gamma Gibbs — synthetic validation against the MLE, a bambi/PyMC cross-check, and the Mroz labor-force application with odds ratios and average partial effects (Python) Open → R R cross-check — bayesm::rmnlIndepMetrop (logit as a 2-alternative MNL by independence Metropolis), on the same synthetic and Mroz data (R) Open →

Downloads

binary_logit_mcmc.py From-scratch samplers — IRLS MLE, RW-Metropolis, Pólya–Gamma Gibbs (NumPy only) mroz_lfp.csv Mroz (1987) labor-force participation — n=753 married women

Sampler Module — Source Code

"""
Bayesian binary logistic regression — from-scratch samplers (NumPy only).

Model:   y_i ~ Bernoulli(p_i),  p_i = sigmoid(x_i' beta),  beta ~ N(0, prior_sd^2 I)

Implemented here:
  * mle_logit       -- frequentist MLE by IRLS / Newton-Raphson (comparison baseline)
  * rw_metropolis   -- random-walk Metropolis on the logit posterior
                       (proposal = c^2 * Sigma_hat, Sigma_hat from the MLE Hessian;
                        c = 2.38/sqrt(k) is the Roberts-Rosenthal optimal scale)

Polya-Gamma Gibbs (the logit analog of Albert-Chib) will be added alongside.
Companion to binary_probit_gibbs.py.
"""
import numpy as np


def sigmoid(eta):
    """Numerically stable logistic function."""
    out = np.empty_like(eta, dtype=float)
    pos = eta >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-eta[pos]))
    e = np.exp(eta[~pos])
    out[~pos] = e / (1.0 + e)
    return out


def _softplus(eta):
    """log(1 + exp(eta)), stable."""
    return np.maximum(eta, 0.0) + np.log1p(np.exp(-np.abs(eta)))


def loglik(beta, X, y):
    eta = X @ beta
    return np.sum(y * eta - _softplus(eta))


def logpost(beta, X, y, prior_prec):
    return loglik(beta, X, y) - 0.5 * prior_prec * np.dot(beta, beta)


def simulate_logit(n, beta, seed=0, intercept=True):
    """Synthetic binary logit data. Returns (X, y, p). beta includes intercept if intercept=True."""
    rng = np.random.default_rng(seed)
    beta = np.asarray(beta, float)
    k = len(beta)
    kx = k - 1 if intercept else k
    Xc = rng.normal(size=(n, kx))
    X = np.column_stack([np.ones(n), Xc]) if intercept else Xc
    p = sigmoid(X @ beta)
    y = (rng.random(n) < p).astype(float)
    return X, y, p


def mle_logit(X, y, tol=1e-10, maxit=100):
    """Logistic-regression MLE by IRLS. Returns (beta_hat, cov, se)."""
    n, k = X.shape
    beta = np.zeros(k)
    for _ in range(maxit):
        p = sigmoid(X @ beta)
        W = p * (1.0 - p)
        grad = X.T @ (y - p)
        H = (X * W[:, None]).T @ X                 # observed information X'WX
        step = np.linalg.solve(H, grad)
        beta = beta + step
        if np.max(np.abs(step)) < tol:
            break
    p = sigmoid(X @ beta)
    cov = np.linalg.inv((X * (p * (1 - p))[:, None]).T @ X)
    return beta, cov, np.sqrt(np.diag(cov))


def rw_metropolis(X, y, R=30000, burn=5000, prior_sd=10.0, scale=None, seed=0):
    """Random-walk Metropolis for Bayesian logit. Proposal cov from the MLE Hessian."""
    rng = np.random.default_rng(seed)
    n, k = X.shape
    prior_prec = 1.0 / prior_sd ** 2
    bhat, cov, se = mle_logit(X, y)
    if scale is None:
        scale = 2.38 / np.sqrt(k)                  # Roberts-Rosenthal optimal RW scale
    L = np.linalg.cholesky(cov) * scale            # proposal sqrt-cov

    beta = bhat.copy()
    lp = logpost(beta, X, y, prior_prec)
    keep = R - burn
    draws = np.zeros((keep, k))
    acc = 0
    for g in range(R):
        prop = beta + L @ rng.standard_normal(k)
        lpp = logpost(prop, X, y, prior_prec)
        if np.log(rng.random()) < lpp - lp:
            beta, lp = prop, lpp
            acc += 1
        if g >= burn:
            draws[g - burn] = beta
    return dict(draws=draws, accept=acc / R, bhat=bhat, mle_cov=cov, mle_se=se,
                scale=scale, prior_sd=prior_sd)


def summarize(draws, names=None):
    """Posterior mean, sd, and 95% CI per coefficient."""
    m = draws.mean(0); s = draws.std(0)
    lo = np.percentile(draws, 2.5, axis=0); hi = np.percentile(draws, 97.5, axis=0)
    k = draws.shape[1]
    names = names or [f'beta_{j}' for j in range(k)]
    return {names[j]: dict(mean=m[j], sd=s[j], q025=lo[j], q975=hi[j]) for j in range(k)}


# ----------------------------------------------------------------------------
# Polya-Gamma data augmentation (the logit analog of Albert-Chib for probit)
# ----------------------------------------------------------------------------
def rpg(b, c, rng, n_terms=128):
    """Sample Polya-Gamma PG(b, c) via the truncated sum-of-Gammas series:

        PG(b, c) =d  (1 / 2 pi^2) * sum_{j=1}^inf  g_j / ((j-1/2)^2 + c^2/(4 pi^2)),
        g_j ~ Gamma(b, 1) iid.

    Vectorized over observations; `c` is an array of shape (n,), `b` scalar/array.
    Truncating at n_terms gives a tiny, negligible bias (dropped tail ~ 1/n_terms).
    (Devroye's alternating-series method is the exact finite-time alternative.)
    """
    c = np.asarray(c, float)
    b = np.broadcast_to(b, c.shape)
    j = np.arange(1, n_terms + 1)
    denom = (j - 0.5) ** 2 + (c[:, None] ** 2) / (4.0 * np.pi ** 2)   # (n, n_terms)
    g = rng.gamma(b[:, None], 1.0, size=(c.shape[0], n_terms))         # Gamma(b, 1)
    return (1.0 / (2.0 * np.pi ** 2)) * np.sum(g / denom, axis=1)


def pg_gibbs(X, y, R=6000, burn=1000, prior_sd=10.0, seed=0, n_terms=128):
    """Polya-Gamma Gibbs for Bayesian logit (Polson, Scott & Windle 2013).

    Augment with omega_i ~ PG(1, x_i'beta); then beta | omega is conjugate Gaussian:
        Omega = diag(omega),  kappa = y - 1/2,
        V = (X' Omega X + B0^{-1})^{-1},   m = V (X' kappa + B0^{-1} b0),
        beta | . ~ N(m, V).
    No accept/reject -- every draw is kept.
    """
    rng = np.random.default_rng(seed)
    k = X.shape[1]
    B0inv = np.eye(k) / prior_sd ** 2
    b0 = np.zeros(k)
    kappa = y - 0.5
    XtKappa = X.T @ kappa
    beta, *_ = mle_logit(X, y)                      # initialise at the MLE

    keep = R - burn
    draws = np.zeros((keep, k))
    for g in range(R):
        omega = rpg(1.0, X @ beta, rng, n_terms)    # PG(1, eta) per observation
        V = np.linalg.inv((X * omega[:, None]).T @ X + B0inv)
        m = V @ (XtKappa + B0inv @ b0)
        beta = m + np.linalg.cholesky(V) @ rng.standard_normal(k)
        if g >= burn:
            draws[g - burn] = beta
    return dict(draws=draws, accept=1.0)            # Gibbs: acceptance is 1 by construction

References

Polson, N. G., Scott, J. G. & Windle, J. (2013). Bayesian inference for logistic models using Pólya–Gamma latent variables. Journal of the American Statistical Association 108(504), 1339–1349. — the Pólya–Gamma Gibbs sampler
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 probit antecedent of the augmentation trick
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics 21(6), 1087–1092. — the original Metropolis algorithm
Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57(1), 97–109. — its generalization — the Metropolis–Hastings sampler used in the random-walk fit
Roberts, G. O., Gelman, A. & Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Annals of Applied Probability 7(1), 110–120. — the 2.38/√d proposal scaling
Mroz, T. A. (1987). The sensitivity of an empirical model of married women's hours of work to economic and statistical assumptions. Econometrica 55(4), 765–799. — the labor-force participation data
Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. — bayesm::rmnlIndepMetrop