Bayesian binary logistic regression — from scratch¶

From-scratch samplers: Random-walk Metropolis & Pólya–Gamma Gibbs¶

References: Metropolis et al. (1953); Hastings (1970); Roberts, Gelman & Gilks (1997).

Section Content
Model Logit model, random-walk Metropolis, MLE baseline, references
Section 1 Synthetic data — RW-Metropolis vs. MLE (correctness check)
Section 2 Pólya–Gamma Gibbs (logit Albert–Chib) — same posterior, better mixing
Section 3 Package cross-check — bambi & PyMC (NUTS)
Section 4 Real data — Mroz (1987) women’s labor-force participation

R / bayesm version in binary_logit_bayesm.ipynb.

From-scratch sampler module: binary_logit_mcmc.py — IRLS MLE, random-walk Metropolis and Pólya–Gamma Gibbs, NumPy only.

Model — Bayesian binary logit¶

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

with a diffuse prior ($\sigma_\beta = 10$). The posterior $\propto \prod_i p_i^{y_i}(1-p_i)^{1-y_i}\, N(\beta;0,\sigma_\beta^2 I)$ has no conjugate form, so we sample it by Metropolis.

Random-walk Metropolis¶

Each step: propose $\beta' = \beta + N(0,\Sigma_{\text{prop}})$, accept with probability $\min\{1,\ \exp[\ell(\beta') - \ell(\beta)]\}$, where $\ell$ is the log-posterior

$$\ell(\beta) = \sum_i \big[y_i\,\eta_i - \text{softplus}(\eta_i)\big] - \tfrac{1}{2\sigma_\beta^2}\beta'\beta, \qquad \eta = X\beta,\ \ \text{softplus}(\eta)=\log(1+e^{\eta}).$$

  • Proposal $\Sigma_{\text{prop}} = c^2\,\hat\Sigma$, with $\hat\Sigma = (X'\hat W X)^{-1}$ the inverse Fisher information at the MLE ($\hat W = \text{diag}(\hat p_i(1-\hat p_i))$) — so the proposal matches the posterior's local shape.
  • Scale $c = 2.38/\sqrt{k}$ — the Roberts–Gelman–Gilks optimal scaling (acceptance $\approx 0.234$ in high dimension).

MLE baseline¶

The frequentist MLE (IRLS / Newton–Raphson) gives $\hat\beta$ and $\text{se}=\sqrt{\text{diag}(\hat\Sigma)}$. Under a diffuse prior the Bernstein–von Mises theorem says the posterior mean $\approx\hat\beta$ and posterior sd $\approx\text{se}$ — so matching the MLE is a correctness check on the sampler (it is not the same as proving the prior matters; with informative priors they would differ).

References¶

  • Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (1953); Hastings (1970)
  • Roberts, Gelman & Gilks (1997) — optimal scaling of random-walk Metropolis ($2.38/\sqrt d$)
  • Gelman, Carlin, Stern, Dunson, Vehtari & Rubin — Bayesian Data Analysis (3rd ed.)
  • Rossi, Allenby & McCulloch (2005) — Bayesian Statistics and Marketing (bayesm)
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from binary_logit_mcmc import simulate_logit, mle_logit, rw_metropolis, pg_gibbs, summarize

print('binary_logit_mcmc loaded')

binary_logit_mcmc loaded

1. Synthetic data: posterior vs. MLE¶

True $\beta = (0.5,\,-1.0,\,1.5,\,-0.75)$ (intercept + 3 covariates), $n=2000$, $x\sim N(0,1)$. Expected: posterior mean $\approx$ MLE $\approx$ truth, acceptance $\approx 0.3$, and 95% credible intervals matching the MLE standard errors.

In [2]:
BETA_TRUE = np.array([0.5, -1.0, 1.5, -0.75])      # intercept + 3 covariates
N = 2000
X, y, p = simulate_logit(N, BETA_TRUE, seed=7)
bhat, cov, se = mle_logit(X, y)
print(f'Synthetic logit: n={N}, k={len(BETA_TRUE)}, mean(y)={y.mean():.3f}')
print('MLE   :', np.round(bhat, 3))
print('MLE se:', np.round(se, 3))

Synthetic logit: n=2000, k=4, mean(y)=0.569
MLE   : [ 0.445 -1.005  1.535 -0.619]
MLE se: [0.059 0.066 0.08  0.061]
In [3]:
out = rw_metropolis(X, y, R=40000, burn=10000, prior_sd=10.0, seed=0)
D = out['draws']
names = ['const', 'x1', 'x2', 'x3']
pm, psd = D.mean(0), D.std(0)

print(f"RW-Metropolis: acceptance={out['accept']:.3f}  (target ~0.23-0.40)\n")
hdr = f"{'param':>6}{'truth':>8}{'MLE':>9}{'MLE se':>8}{'postmean':>10}{'postsd':>8}"
print(hdr); print('-'*len(hdr))
for j, nm in enumerate(names):
    print(f"{nm:>6}{BETA_TRUE[j]:8.2f}{bhat[j]:9.3f}{se[j]:8.3f}{pm[j]:10.3f}{psd[j]:8.3f}")
print(f"\nmax|postmean - MLE| = {np.max(np.abs(pm-bhat)):.4f}")
print(f"max|postsd  - MLEse| = {np.max(np.abs(psd-se)):.4f}")

RW-Metropolis: acceptance=0.303  (target ~0.23-0.40)

 param   truth      MLE  MLE se  postmean  postsd
-------------------------------------------------
 const    0.50    0.445   0.059     0.446   0.058
    x1   -1.00   -1.005   0.066    -1.008   0.067
    x2    1.50    1.535   0.080     1.542   0.080
    x3   -0.75   -0.619   0.061    -0.623   0.060

max|postmean - MLE| = 0.0067
max|postsd  - MLEse| = 0.0007
In [4]:
k = D.shape[1]
fig, axes = plt.subplots(2, k, figsize=(3.2*k, 5))
for j in range(k):
    axes[0, j].plot(D[:, j], lw=0.3, color='steelblue')
    axes[0, j].set_title(names[j]); axes[0, j].set_xlabel('iteration')
    axes[1, j].hist(D[:, j], bins=50, density=True, color='steelblue', alpha=0.6)
    axes[1, j].axvline(BETA_TRUE[j], color='black', ls='--', lw=1.5, label='truth')
    axes[1, j].axvline(bhat[j], color='firebrick', ls='-', lw=1.2, label='MLE')
    axes[1, j].legend(fontsize=7); axes[1, j].set_xlabel(names[j])
axes[0, 0].set_ylabel('trace'); axes[1, 0].set_ylabel('posterior density')
fig.suptitle('RW-Metropolis logit: traces (top) and posteriors vs. truth & MLE (bottom)', fontsize=11)
plt.tight_layout()
plt.savefig(r'C:\Users\user\project\binary_logit_rwm.png', dpi=120, bbox_inches='tight')
plt.show()
No description has been provided for this image

1. Results — RW-Metropolis vs. MLE¶

param truth MLE MLE se post. mean post. sd
const 0.50 0.445 0.059 0.446 0.058
x1 −1.00 −1.005 0.066 −1.008 0.067
x2 1.50 1.535 0.080 1.542 0.080
x3 −0.75 −0.619 0.061 −0.623 0.060
  • Posterior ≈ MLE to within 0.007 on the means and 0.001 on the sds — exactly the Bernstein–von Mises result under a diffuse prior, which validates the from-scratch sampler. (β₃'s gap from truth, −0.62 vs −0.75, is in the data: the MLE shows the same −0.62, i.e. ordinary sampling variation at n=2000, not a sampler error.)
  • Acceptance ≈ 0.30, squarely in the optimal band — the proposal is the inverse-Fisher covariance scaled by $2.38/\sqrt{k}$.
  • Lag-1 autocorrelation ≈ 0.86 per coefficient: the usual cost of random-walk Metropolis. With 30k retained draws the effective sample size is still ample, but this is precisely the inefficiency that the Pólya–Gamma Gibbs sampler (Section 2, next) removes — it draws every iteration with no accept/reject and far lower autocorrelation.

Takeaway: a hand-written random-walk Metropolis recovers the logit posterior and reproduces the MLE under a flat prior. Next we add the Pólya–Gamma data-augmentation Gibbs (the logit Albert–Chib), then cross-check both against bambi/PyMC.

2. Pólya–Gamma Gibbs (the logit Albert–Chib)¶

Random-walk Metropolis (Section 1) works but mixes slowly (lag-1 autocorr ≈ 0.86) and needs a tuned proposal. Pólya–Gamma augmentation (Polson, Scott & Windle 2013) makes the logit likelihood conditionally Gaussian — the exact logit analog of Albert–Chib for probit — giving a tuning-free Gibbs sampler.

Introduce a latent $\omega_i \sim \mathrm{PG}(1, x_i'\beta)$ per observation. Conditional on $\omega$, the likelihood in $\beta$ is Gaussian, so the two blocks are:

$$\omega_i \mid \beta \sim \mathrm{PG}(1,\ x_i'\beta), \qquad \beta \mid \omega \sim N(m, V),\quad V=(X'\Omega X + B_0^{-1})^{-1},\ \ m=V(X'\kappa + B_0^{-1}b_0),$$

with $\Omega=\mathrm{diag}(\omega)$ and $\kappa = y - \tfrac12$. No accept/reject — every draw is kept.

The PG draw, from scratch. rpg samples $\mathrm{PG}(1,c)$ from its series representation

$$\mathrm{PG}(1,c)\ \stackrel{d}{=}\ \frac{1}{2\pi^2}\sum_{j=1}^{\infty}\frac{g_j}{(j-\tfrac12)^2 + c^2/(4\pi^2)},\qquad g_j\sim\mathrm{Gamma}(1,1),$$

truncated at many terms (negligible bias). Devroye's alternating-series method is the exact finite-time alternative used by BayesLogit / pypolyagamma.

References: Polson, Scott & Windle (2013, JASA 108, 1339–1349); Albert & Chib (1993) — the probit antecedent.

In [5]:
pg = pg_gibbs(X, y, R=6000, burn=1000, prior_sd=10.0, seed=0)   # ~20s
Dp = pg['draws']; pm_pg, psd_pg = Dp.mean(0), Dp.std(0)

print('Pólya–Gamma Gibbs (no accept/reject)\n')
hdr = f"{'param':>6}{'truth':>8}{'MLE':>9}{'RWM mean':>10}{'PG mean':>9}{'PG sd':>8}"
print(hdr); print('-'*len(hdr))
for j, nm in enumerate(names):
    print(f"{nm:>6}{BETA_TRUE[j]:8.2f}{bhat[j]:9.3f}{D[:,j].mean():10.3f}{pm_pg[j]:9.3f}{psd_pg[j]:8.3f}")

ac = lambda Z: [np.corrcoef(Z[:-1, j], Z[1:, j])[0, 1] for j in range(Z.shape[1])]
print('\nlag-1 autocorr  RWM:', np.round(ac(D), 2))
print('lag-1 autocorr  PG :', np.round(ac(Dp), 2))
print('max|PG mean - MLE| =', round(np.max(np.abs(pm_pg - bhat)), 4))
print('max|PG mean - RWM mean| =', round(np.max(np.abs(pm_pg - D.mean(0))), 4))

Pólya–Gamma Gibbs (no accept/reject)

 param   truth      MLE  RWM mean  PG mean   PG sd
--------------------------------------------------
 const    0.50    0.445     0.446    0.449   0.059
    x1   -1.00   -1.005    -1.008   -1.016   0.066
    x2    1.50    1.535     1.542    1.552   0.081
    x3   -0.75   -0.619    -0.623   -0.625   0.061

lag-1 autocorr  RWM: [0.86 0.86 0.86 0.86]
lag-1 autocorr  PG : [0.28 0.39 0.51 0.36]
max|PG mean - MLE| = 0.0167
max|PG mean - RWM mean| = 0.01
In [6]:
k = Dp.shape[1]
fig, axes = plt.subplots(2, k, figsize=(3.2*k, 5))
for j in range(k):
    # row 0: posterior densities, PG vs RWM, with truth & MLE
    axes[0, j].hist(D[:, j], bins=50, density=True, alpha=0.45, color='steelblue', label='RWM')
    axes[0, j].hist(Dp[:, j], bins=50, density=True, alpha=0.45, color='seagreen', label='PG')
    axes[0, j].axvline(BETA_TRUE[j], color='black', ls='--', lw=1.3)
    axes[0, j].axvline(bhat[j], color='firebrick', ls='-', lw=1.0)
    axes[0, j].set_title(names[j]); axes[0, j].legend(fontsize=7)
    # row 1: first 800 iterations of each chain (mixing)
    axes[1, j].plot(D[:800, j], lw=0.4, color='steelblue', alpha=0.8, label='RWM')
    axes[1, j].plot(Dp[:800, j], lw=0.4, color='seagreen', alpha=0.8, label='PG')
    axes[1, j].set_xlabel('iteration'); axes[1, j].legend(fontsize=7)
axes[0, 0].set_ylabel('posterior density'); axes[1, 0].set_ylabel('trace (first 800)')
fig.suptitle('Pólya–Gamma Gibbs vs. RW-Metropolis: same posterior, better mixing', fontsize=11)
plt.tight_layout()
plt.savefig(r'C:\Users\user\project\binary_logit_pg.png', dpi=120, bbox_inches='tight')
plt.show()
No description has been provided for this image

2. Results — Pólya–Gamma Gibbs vs. RW-Metropolis¶

param truth MLE RWM mean PG mean PG sd
const 0.50 0.445 0.446 0.449 0.059
x1 −1.00 −1.005 −1.008 −1.016 0.066
x2 1.50 1.535 1.542 1.552 0.081
x3 −0.75 −0.619 −0.623 −0.625 0.061
  • Same posterior, three ways. PG-Gibbs matches the MLE to within 0.017 and the RW-Metropolis means to within 0.01 — the data-augmentation sampler and the direct Metropolis sampler agree, which validates both.
  • Much better mixing. Lag-1 autocorrelation drops from ≈0.86 (RWM) to ≈0.28–0.51 (PG), and there is no accept/reject and no proposal to tune — the practical payoff of the Pólya–Gamma trick. (RWM is cheaper per iteration here because the from-scratch PG draw sums many Gammas; PG wins on robustness and tuning-free use, and scales better as the model grows.)
  • The probit/logit symmetry. This is the logit counterpart of your Albert–Chib probit sampler: probit augments with latent truncated normals, logit augments with latent Pólya–Gamma variables — in both cases the augmented model is conditionally Gaussian and β is drawn from a normal full conditional. It is exactly why bayesm does probit by Gibbs but logit by Metropolis (no clean Gibbs for logit without Pólya–Gamma).

Next: cross-check both from-scratch samplers against bambi and PyMC (Section 3), then the R bayesm version, then a real-data application.

3. Package cross-check: bambi and PyMC¶

A final check: fit the same synthetic logit with two mainstream Bayesian tools and confirm they land on the posterior our two from-scratch samplers (and the MLE) already pinned down.

  • bambi — formula API (y ~ x1 + x2 + x3, family="bernoulli") on top of PyMC; the R-style one-liner.
  • PyMC — the model written out directly: β ~ N(0, 10), Bernoulli(logit_p = Xβ).

Both use NUTS (gradient-based HMC). Backend note: PyMC's default sampler compiles via PyTensor's C backend, which needs a C++ compiler (g++); if that's unavailable you'll see a g++ not available warning. To avoid it entirely we use the NumPyro/JAX backend (nuts_sampler="numpyro" / inference_method="numpyro"), which needs no compiler. Switch to the default ("mcmc") if your environment has gxx installed.

In [7]:
import bambi as bmb
df = pd.DataFrame({'y': y, 'x1': X[:, 1], 'x2': X[:, 2], 'x3': X[:, 3]})

bmodel = bmb.Model("y ~ x1 + x2 + x3", df, family="bernoulli")
# NumPyro backend (no C compiler needed); use inference_method="mcmc" if you have gxx
bidata = bmodel.fit(inference_method="numpyro", draws=1000, tune=1000,
                    chains=2, random_seed=0, progressbar=False)
po = bidata.posterior
bambi_mean = np.array([float(po['Intercept'].mean())] +
                      [float(po[v].mean()) for v in ['x1', 'x2', 'x3']])
bambi_sd = np.array([float(po['Intercept'].std())] +
                    [float(po[v].std()) for v in ['x1', 'x2', 'x3']])
print('bambi posterior mean:', np.round(bambi_mean, 3))
print('bambi posterior sd  :', np.round(bambi_sd, 3))

Modeling the probability that y==1
c:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
NUTS[numpyro]: [Intercept, x1, x2, x3]
We recommend running at least 4 chains for robust computation of convergence diagnostics

bambi posterior mean: [ 0.444 -1.     1.527 -0.616]
bambi posterior sd  : [0.057 0.065 0.078 0.058]
In [8]:
import pymc as pm

with pm.Model() as pmodel:
    beta = pm.Normal('beta', 0.0, 10.0, shape=X.shape[1])
    pm.Bernoulli('y_obs', logit_p=pm.math.dot(X, beta), observed=y)
    # NumPyro/JAX backend -> no C compiler needed (default "mcmc" needs g++)
    pidata = pm.sample(1000, tune=1000, chains=2, random_seed=0,
                       nuts_sampler="numpyro", progressbar=False,
                       compute_convergence_checks=False)
pymc_mean = pidata.posterior['beta'].mean(('chain', 'draw')).values
pymc_sd = pidata.posterior['beta'].std(('chain', 'draw')).values

print('\n      param   truth      MLE   RWM mean   PG mean  bambi   PyMC')
for j, nm in enumerate(names):
    print(f"{nm:>11}{BETA_TRUE[j]:8.2f}{bhat[j]:9.3f}{D[:,j].mean():10.3f}"
          f"{Dp[:,j].mean():10.3f}{bambi_mean[j]:8.3f}{pymc_mean[j]:8.3f}")
print('\nmax|bambi - MLE| =', round(np.max(np.abs(bambi_mean - bhat)), 4))
print('max|PyMC  - MLE| =', round(np.max(np.abs(pymc_mean - bhat)), 4))

c:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
NUTS[numpyro]: [beta]

      param   truth      MLE   RWM mean   PG mean  bambi   PyMC
      const    0.50    0.445     0.446     0.449   0.444   0.444
         x1   -1.00   -1.005    -1.008    -1.016  -1.000  -1.005
         x2    1.50    1.535     1.542     1.552   1.527   1.539
         x3   -0.75   -0.619    -0.623    -0.625  -0.616  -0.620

max|bambi - MLE| = 0.0075
max|PyMC  - MLE| = 0.0045

3. Results — all roads lead to the same posterior¶

Five routes to the binary-logit posterior — two from scratch, two packaged, plus the MLE baseline — agree to within Monte-Carlo error (posterior means ≈ const 0.45, x1 −1.01, x2 1.54, x3 −0.62; the printed table fills the exact bambi/PyMC columns on your run):

Method Type Discrete/aux handling Tuning
MLE (IRLS) frequentist — —
RW-Metropolis Bayesian, from scratch direct MH on the posterior proposal scale
Pólya–Gamma Gibbs Bayesian, from scratch PG augmentation → conjugate none
bambi Bayesian, package NUTS (formula API) auto (warmup)
PyMC Bayesian, package NUTS (explicit model) auto (warmup)

Reading: the from-scratch samplers are validated against production tools (and vice-versa) — convergence of independent implementations is the strongest correctness evidence. Practically:

  • bambi — fastest to write (one formula), great default priors, ideal for routine work.
  • PyMC — full control of priors/structure; the path to hierarchical or non-standard logits.
  • RW-Metropolis / PG-Gibbs (scratch) — for understanding the machinery, and for settings where you want a dependency-free, transparent sampler. PG-Gibbs in particular is the logit twin of your Albert–Chib probit.

Backend caveat: these two cells need the pymc-env kernel (where bambi/pymc/numpyro live) and use the NumPyro backend to avoid the PyTensor C-compiler requirement. Sections 1–2 (pure NumPy) run anywhere.

Next: the R bayesm version (binary_logit_bayesm.ipynb, logit by independence Metropolis), then a real-data application.

4. Real data: Mroz (1987) women's labor-force participation¶

The canonical binary-choice dataset (Mroz 1987; Wooldridge, Introductory Econometrics, Example 17.1). $n=753$ married women; inlf $=1$ if in the labor force. Model:

$$\text{inlf} \sim \text{nwifeinc} + \text{educ} + \text{exper} + \text{exper}^2 + \text{age} + \text{kidslt6} + \text{kidsge6}$$

We fit it with all three from-scratch routes (MLE, RW-Metropolis, Pólya–Gamma Gibbs) and check against Wooldridge's published logit estimates. Data: mroz_lfp.csv.

Variables (Mroz 1987; 1976 PSID, 1975 calendar year; married white women aged 30–60, $n=753$):

Variable Definition Units mean range exp. sign
inlf (outcome) =1 if in the labor force in 1975 0/1 0.568 {0,1} —
nwifeinc Non-wife income = (family income − wife's earnings) $1000s 20.13 −0.03…96 − (income effect)
educ Years of schooling years 12.29 5…17 + (↑ market wage)
exper Actual labor-market experience years 10.63 0…45 + (attachment)
expersq exper² years² 178.0 0…2025 − (concave)
age Age years 42.54 30…60 − (toward retirement)
kidslt6 Number of children < 6 count 0.238 0…3 − (childcare)
kidsge6 Number of children 6–18 count 1.353 0…8 ≈ 0

nwifeinc is the household budget available without the wife working; exper+expersq give a concave experience profile. Coefficients below are partial effects (each controlling for the rest).

In [9]:
mroz = pd.read_csv('mroz_lfp.csv')
mcols = ['nwifeinc', 'educ', 'exper', 'expersq', 'age', 'kidslt6', 'kidsge6']
Xm = np.column_stack([np.ones(len(mroz))] + [mroz[c].values for c in mcols])
ym = mroz['inlf'].values.astype(float)
mnames = ['const'] + mcols

bm, covm, sem = mle_logit(Xm, ym)
rwm_m = rw_metropolis(Xm, ym, R=40000, burn=10000, seed=0)
pg_m  = pg_gibbs(Xm, ym, R=6000, burn=1000, seed=0)
Drm, Dpm = rwm_m['draws'], pg_m['draws']

wool = {'const':0.425,'nwifeinc':-0.021,'educ':0.221,'exper':0.206,
        'expersq':-0.003,'age':-0.088,'kidslt6':-1.443,'kidsge6':0.060}
print(f"n={len(mroz)}  in-labor-force rate={ym.mean():.3f}  RWM acceptance={rwm_m['accept']:.3f}\n")
hdr = f"{'param':>9}{'Wooldridge':>11}{'MLE':>9}{'RWM':>9}{'PG':>9}{'post sd':>9}"
print(hdr); print('-'*len(hdr))
for j, nm in enumerate(mnames):
    print(f"{nm:>9}{wool[nm]:11.3f}{bm[j]:9.3f}{Drm[:,j].mean():9.3f}{Dpm[:,j].mean():9.3f}{Dpm[:,j].std():9.3f}")

n=753  in-labor-force rate=0.568  RWM acceptance=0.274

    param Wooldridge      MLE      RWM       PG  post sd
--------------------------------------------------------
    const      0.425    0.425    0.451    0.441    0.870
 nwifeinc     -0.021   -0.021   -0.021   -0.022    0.009
     educ      0.221    0.221    0.224    0.225    0.044
    exper      0.206    0.206    0.207    0.208    0.032
  expersq     -0.003   -0.003   -0.003   -0.003    0.001
      age     -0.088   -0.088   -0.090   -0.090    0.015
  kidslt6     -1.443   -1.443   -1.457   -1.467    0.208
  kidsge6      0.060    0.060    0.060    0.059    0.076
In [10]:
# Forest plot: PG posterior mean +/- 95% CI vs MLE (slopes only; const has a huge CI)
sl = list(range(1, len(mnames)))            # drop const for readability
lo = np.percentile(Dpm, 2.5, 0); hi = np.percentile(Dpm, 97.5, 0); mean = Dpm.mean(0)
yy = np.arange(len(sl))
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.errorbar(mean[sl], yy, xerr=[mean[sl]-lo[sl], hi[sl]-mean[sl]], fmt='o',
            color='seagreen', capsize=3, label='PG posterior (95% CI)')
ax.scatter(bm[sl], yy+0.14, color='firebrick', marker='s', zorder=5, label='MLE')
ax.axvline(0, color='gray', ls=':')
ax.set_yticks(yy); ax.set_yticklabels([mnames[j] for j in sl]); ax.invert_yaxis()
ax.set_xlabel('coefficient (log-odds)'); ax.legend(fontsize=8)
ax.set_title('Mroz LFP logit — Pólya–Gamma posterior vs. MLE (slopes)')
plt.tight_layout()
plt.savefig(r'C:\Users\user\project\binary_logit_mroz.png', dpi=120, bbox_inches='tight')
plt.show()
No description has been provided for this image

4. Results — Mroz LFP¶

param Wooldridge MLE RWM PG odds ratio $e^\beta$
const 0.425 0.425 0.451 0.441 —
nwifeinc −0.021 −0.021 −0.021 −0.022 0.98 per $1k
educ 0.221 0.221 0.224 0.225 1.25 per year
exper 0.206 0.206 0.207 0.208 1.23 (at 0 exper)
expersq −0.003 −0.003 −0.003 −0.003 concave in experience
age −0.088 −0.088 −0.090 −0.090 0.91 per year
kidslt6 −1.443 −1.443 −1.457 −1.467 0.23 per child < 6
kidsge6 0.060 0.060 0.060 0.059 1.06 (n.s.)
  • All three from-scratch routes reproduce Wooldridge's published logit estimates (MLE to 3 decimals; RWM/PG posteriors ≈ MLE under the diffuse prior), validating the samplers on real data. RWM acceptance ≈ 0.27.
  • Economics (the classic story): a young child (kidslt6) is the dominant force — each one multiplies the odds of labor-force participation by $e^{-1.44}\approx 0.23$. Education raises participation ($\times 1.25$ per year), experience helps at a decreasing rate (positive exper, negative expersq), and both non-wife income and age reduce it modestly. Children over 6 (kidsge6) have no significant effect.
  • The Bayesian payoff: we get the full posterior of every coefficient — and hence of any odds ratio or predicted probability — with credible intervals, not just point estimates ± SEs. The forest plot shows the PG-Gibbs 95% credible intervals against the MLE points.

This closes the binary-logit arc: from-scratch RW-Metropolis & Pólya–Gamma Gibbs, cross-checked against bambi/PyMC/bayesm, validated on synthetic data and now on the canonical Mroz LFP dataset.

In [11]:
# --- Analysis of findings: significance, odds ratios, predicted probabilities, partial effects ---
from binary_logit_mcmc import sigmoid
lo = np.percentile(Dpm, 2.5, 0); hi = np.percentile(Dpm, 97.5, 0); pmean = Dpm.mean(0)

print('coefficient     post mean        95% CI         signif    odds ratio e^beta [95% CI]')
for j, nm in enumerate(mnames):
    sig = 'yes' if (lo[j] > 0 or hi[j] < 0) else 'no '
    print(f"{nm:>9}   {pmean[j]:8.3f}   [{lo[j]:7.3f},{hi[j]:7.3f}]   {sig}    "
          f"{np.exp(pmean[j]):6.3f} [{np.exp(lo[j]):.3f}, {np.exp(hi[j]):.3f}]")

xbar = Xm.mean(0)
P = lambda xrow: sigmoid(Dpm @ xrow)
def prof(**kw):
    x = xbar.copy()
    for k, v in kw.items(): x[mnames.index(k)] = v
    return x
print(f"\nPredicted P(inlf) at sample means = {P(xbar).mean():.3f}")
print(f"  educ 12 -> 16 (others at mean):  {P(prof(educ=12)).mean():.3f} -> {P(prof(educ=16)).mean():.3f}")
print(f"  kids<6  0 -> 1 (others at mean):  {P(prof(kidslt6=0)).mean():.3f} -> {P(prof(kidslt6=1)).mean():.3f}")

# average partial effects: change averaged over all women AND over the posterior
base = sigmoid(Xm @ Dpm.T)
def ame(col, delta):
    X2 = Xm.copy(); X2[:, mnames.index(col)] += delta
    e = (sigmoid(X2 @ Dpm.T) - base).mean(0)
    return e.mean(), np.percentile(e, 2.5), np.percentile(e, 97.5)
print()
for col, delta, lab in [('kidslt6', 1, '+1 child < 6'), ('educ', 1, '+1 year educ'), ('age', 10, '+10 yrs age')]:
    m, l, h = ame(col, delta)
    print(f"  AME {lab:>13}: {m:+.3f}  [{l:+.3f}, {h:+.3f}]")
acc = ((sigmoid(Xm @ pmean) > 0.5).astype(int) == ym).mean()
print(f"\nin-sample accuracy (0.5 cutoff) = {acc:.3f}   (base rate {max(ym.mean(), 1-ym.mean()):.3f})")

coefficient     post mean        95% CI         signif    odds ratio e^beta [95% CI]
    const      0.441   [ -1.262,  2.126]   no      1.555 [0.283, 8.382]
 nwifeinc     -0.022   [ -0.039, -0.006]   yes     0.978 [0.962, 0.994]
     educ      0.225   [  0.138,  0.315]   yes     1.253 [1.148, 1.370]
    exper      0.208   [  0.146,  0.272]   yes     1.231 [1.158, 1.313]
  expersq     -0.003   [ -0.005, -0.001]   yes     0.997 [0.995, 0.999]
      age     -0.090   [ -0.119, -0.061]   yes     0.914 [0.888, 0.941]
  kidslt6     -1.467   [ -1.892, -1.075]   yes     0.231 [0.151, 0.341]
  kidsge6      0.059   [ -0.089,  0.209]   no      1.061 [0.915, 1.233]

Predicted P(inlf) at sample means = 0.584
  educ 12 -> 16 (others at mean):  0.568 -> 0.763
  kids<6  0 -> 1 (others at mean):  0.665 -> 0.316

  AME  +1 child < 6: -0.257  [-0.316, -0.197]
  AME  +1 year educ: +0.039  [+0.025, +0.053]
  AME   +10 yrs age: -0.160  [-0.205, -0.112]

in-sample accuracy (0.5 cutoff) = 0.736   (base rate 0.568)

4. Analysis of findings — Mroz LFP¶

What drives labor-force participation? With the full posterior in hand we can read significance, odds ratios, predicted probabilities, and average partial effects — each with credible intervals.

Significance (95% credible interval excludes 0): six of seven predictors are credibly nonzero — nwifeinc, educ, exper, expersq, age, kidslt6. Only kidsge6 (children 6–18) is not credible (β 95% CI ≈ [−0.09, 0.21], OR ≈ 1.06) — school-age children don't measurably affect the mother's participation once everything else is controlled. (const's wide CI just reflects extrapolation to all-zero covariates.)

Odds ratios $e^\beta$:

  • kidslt6 0.23 — each child under 6 multiplies the odds of working by ~0.23 (a ~77% drop). By far the strongest effect.
  • educ 1.25 / year, exper 1.23 (concave: expersq 0.997), age 0.91 / year, nwifeinc 0.98 / $1,000.

Predicted probabilities (more interpretable than odds; others held at sample means):

  • At the average profile, $\widehat{\Pr}(\text{inlf}) \approx 0.58$.
  • High school → college (educ 12→16): 0.57 → 0.76 (+19 pp).
  • First child under 6 (kidslt6 0→1): 0.67 → 0.32 (−35 pp at this profile).

Average partial effects (the discrete/continuous change averaged over all women and the posterior — the policy-relevant magnitudes, since logit coefficients aren't probabilities):

  • +1 child < 6: −0.26 [−0.32, −0.20] — a young child lowers participation probability by ~26 percentage points on average.
  • +1 year of education: +0.039 [+0.025, +0.053] (~4 pp per year).
  • +10 years of age: −0.16 [−0.21, −0.11] — about a 16-percentage-point fall in participation (the move toward retirement).

Fit: 73.6% in-sample accuracy at a 0.5 cutoff, vs. a 56.8% base rate.

The Bayesian payoff: every quantity above — odds ratios, predicted probabilities, and average partial effects — carries a credible interval, because coefficient uncertainty is propagated automatically through the posterior. A frequentist logit gives the same point estimates but needs the delta method (or bootstrapping) for CIs on these nonlinear functions; here they fall out of the same posterior draws.