Bayesian Binary Probit Model¶

Albert & Chib (1993) Gibbs sampler — Bayesian vs MLE comparison¶

Model: $$y_i = \mathbf{1}\{z_i > 0\}, \quad z_i = x_i'\beta + \varepsilon_i, \quad \varepsilon_i \sim N(0,1)$$ $$P(y_i = 1 \mid x_i) = \Phi(x_i'\beta)$$

Error variance fixed at 1 for identification (probit normalisation).

Gibbs sampler — two blocks (Albert & Chib 1993):

Block Update
1. $z_i \mid y_i, \beta$ $TN(x_i'\beta, 1;\; 0, \infty)$ if $y_i=1$; $TN(x_i'\beta, 1;\; -\infty, 0)$ if $y_i=0$
2. $\beta \mid z$ $N(\bar{b}, \bar{B})$ — Normal regression posterior with $\bar{B}^{-1} = X'X + B_0^{-1}$

Note: $\bar{B}^{-1}$ is constant across iterations (σ²=1 fixed) — the Cholesky factorisation is computed once, making the sampler very fast.

Structure:

Section Topic
Section 1 Synthetic data ($n=400$): Gibbs posterior vs MLE — β estimates and credible/confidence intervals
Section 2 Bayesian residuals (AC 1995): CPO and latent residuals for outlier detection on random label flips
Section 3 When CPO outperforms Pearson: high-leverage masking with small $n=30$

References:

  • Albert & Chib (1993), JASA 88(422), 669–679 — Gibbs sampler via data augmentation
  • Albert & Chib (1995), Biometrika 82(4), 747–759 — Bayesian residual analysis
In [1]:
import os, sys

conda_lib = os.path.join(os.environ.get('CONDA_PREFIX', ''), 'Library', 'bin')
if os.path.isdir(conda_lib) and conda_lib not in os.environ.get('PATH', ''):
    os.environ['PATH'] = conda_lib + ';' + os.environ.get('PATH', '')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.special import ndtr

sys.path.insert(0, r'c:\Users\user\project')
from binary_probit_gibbs import (
    binary_probit_gibbs, probit_mle, bayesian_residuals, posterior_summary
)
print('Import OK')

Import OK

1. Synthetic data: Gibbs posterior vs MLE¶

Design: $k=4$ (intercept + 3 covariates), $n=400$, $\beta^\star = [0.5,\; 1.2,\; -0.8,\; 0.5]$.

To create a challenging detection problem, 25 observations have their labels deliberately flipped ($y_i \leftarrow 1 - y_i$). These represent mislabelled data or structural anomalies — observations that contradict the model's predictions.

Goal of Section 1: recover $\beta$ despite label noise; compare Bayesian posterior with MLE. Goal of Section 2: detect the planted outliers using Bayesian residual diagnostics.

In [2]:
# ── Generate synthetic data ───────────────────────────────────────────────────
rng_data = np.random.default_rng(42)
n, k = 400, 4
beta_true = np.array([0.5, 1.2, -0.8, 0.5])

X = np.column_stack([np.ones(n), rng_data.standard_normal((n, k - 1))])
z_latent = X @ beta_true + rng_data.standard_normal(n)
y_clean  = (z_latent > 0).astype(int)

# Plant 25 label-flipped outliers
n_out    = 25
rng_out  = np.random.default_rng(99)
outlier_idx = rng_out.choice(n, size=n_out, replace=False)
y        = y_clean.copy()
y[outlier_idx] = 1 - y[outlier_idx]

is_outlier = np.zeros(n, dtype=bool)
is_outlier[outlier_idx] = True

print(f'n={n}, k={k}, beta_true={beta_true}')
print(f'Clean y=1 rate: {y_clean.mean():.3f}  |  Noisy y=1 rate: {y.mean():.3f}')
print(f'Outliers: {n_out} label-flipped observations')
print(f'Avg |z_latent| at outlier sites: {np.abs(z_latent[outlier_idx]).mean():.2f} sigma')

# Save design matrix + outcomes for bayesm notebook
df_save = pd.DataFrame(X, columns=['x0', 'x1', 'x2', 'x3'])
df_save['y'] = y
df_save['is_outlier'] = is_outlier.astype(int)
df_save.to_csv(r'c:\Users\user\project\binary_probit_data.csv',
               index=False)
print('Saved binary_probit_data.csv')

n=400, k=4, beta_true=[ 0.5  1.2 -0.8  0.5]
Clean y=1 rate: 0.593  |  Noisy y=1 rate: 0.590
Outliers: 25 label-flipped observations
Avg |z_latent| at outlier sites: 1.98 sigma
Saved binary_probit_data.csv
In [3]:
# ── Run Gibbs sampler and MLE ─────────────────────────────────────────────────
out = binary_probit_gibbs(y, X, R=15000, burn=3000, seed=1, nprint=0)
mle = probit_mle(y, X)

param_names = ['beta[0]', 'beta[1]', 'beta[2]', 'beta[3]']
summ = posterior_summary(out['beta'], param_names, true_vals=beta_true)

print('Gibbs posterior:')
print(summ)
print()
print('MLE estimates:')
for j, nm in enumerate(param_names):
    print(f'  {nm}: {mle["beta"][j]:7.4f}  (SE={mle["se"][j]:.4f})')
print(f'  Log-lik: {mle["loglik"]:.2f}')

Done in 0.0 min  |  kept draws: 12000
Gibbs posterior:
           mean      sd    q025    q975  true  covered
beta[0]  0.2851  0.0719  0.1461  0.4264   0.5    False
beta[1]  0.6278  0.0751  0.4820  0.7788   1.2    False
beta[2] -0.4720  0.0796 -0.6289 -0.3157  -0.8    False
beta[3]  0.3557  0.0758  0.2081  0.5056   0.5     True

MLE estimates:
  beta[0]:  0.2831  (SE=0.0710)
  beta[1]:  0.6219  (SE=0.0753)
  beta[2]: -0.4679  (SE=0.0792)
  beta[3]:  0.3534  (SE=0.0759)
  Log-lik: -210.33
In [4]:
# ── Plot: beta posteriors (Gibbs) with MLE overlaid ──────────────────────────
fig, axes = plt.subplots(1, 4, figsize=(16, 4))

for j, ax in enumerate(axes):
    ax.hist(out['beta'][:, j], bins=70, density=True,
            alpha=0.7, color='steelblue', label='Gibbs posterior')
    ax.axvline(beta_true[j], color='k', ls='--', lw=1.8,
               label=f'True {beta_true[j]}')
    ax.axvline(mle['beta'][j], color='red', ls='-', lw=1.5,
               label='MLE')
    ax.axvspan(mle['beta'][j] - 2 * mle['se'][j],
               mle['beta'][j] + 2 * mle['se'][j],
               alpha=0.15, color='red', label='MLE ±2 SE')
    ax.set_title(param_names[j])
    ax.legend(fontsize=7)

fig.suptitle(
    f'§1 Synthetic binary probit — Gibbs posterior (blue) vs MLE (red)\n'
    f'n={n}, k={k}, {n_out} label-flipped outliers  |  true beta={beta_true}',
    fontsize=10)
plt.tight_layout(); plt.show()
No description has been provided for this image

2. Bayesian residuals for outlier detection¶

Albert & Chib (1995): latent residual diagnostics¶

The data-augmentation Gibbs sampler draws the latent utilities $z_i$ at every iteration, making latent residuals $e_i = z_i - x_i'\beta$ directly available as posterior samples.

Under the correctly specified model, $e_i \sim N(0,1)$. For outliers (label-flipped observations):

  • Mislabelled $y_i=1$ (true $y=0$): the model expects $z_i < 0$, but the likelihood forces $z_i > 0$. The sampler squeezes $z_i$ just above 0, giving $\bar{e}_i = E[z_i - x_i'\beta \mid y] \gg 0$.
  • Mislabelled $y_i=0$ (true $y=1$): $z_i$ is squeezed just below 0, giving $\bar{e}_i \ll 0$.

The CPO (Conditional Predictive Ordinate) provides a proper model-based measure: $$\text{CPO}_i = p(y_i \mid y_{-i}) \approx \left(\frac{1}{R} \sum_r \frac{1}{p(y_i \mid \beta^{(r)})}\right)^{-1}$$

Low $\text{CPO}_i$ means the model — trained on all other observations — assigns low probability to $y_i$: a genuine outlier signal.

Classical alternative: Pearson residual $r_i = y_i - \hat{p}_i$ where $\hat{p}_i = \Phi(x_i'\hat{\beta}_\text{MLE})$. These are bounded in $(-1,1)$, ignore parameter uncertainty, and have no proper probabilistic interpretation.

In [5]:
# ── Compute Bayesian residuals ────────────────────────────────────────────────
resid = bayesian_residuals(y, X, out['beta'], out['z'])
resid['is_outlier'] = is_outlier

# Classical Pearson residual (MLE-based, for comparison)
p_hat_mle = ndtr(X @ mle['beta'])
resid['pearson'] = (y - p_hat_mle).round(4)

print('Bayesian residual summary:')
print(resid[['y', 'p_bar', 'CPO', 'e_bar', 'surprise', 'is_outlier']].describe().round(4))

print(f'\nLPML (log pseudo-marginal likelihood) = {resid["log_CPO"].sum():.2f}')
print(f'\nPlanted outliers in bottom-25 CPO:     ',
      resid.nsmallest(25, 'CPO')['is_outlier'].sum(), '/ 25')
print(f'Planted outliers in top-25 |Pearson|:  ',
      resid.reindex(resid['pearson'].abs().nlargest(25).index)['is_outlier'].sum(), '/ 25')
print(f'Planted outliers in top-25 |e_bar|:    ',
      resid.reindex(resid['e_bar'].abs().nlargest(25).index)['is_outlier'].sum(), '/ 25')

Bayesian residual summary:
              y     p_bar       CPO     e_bar  surprise
count  400.0000  400.0000  400.0000  400.0000  400.0000
mean     0.5900    0.5871    0.6460   -0.0000    0.3505
std      0.4924    0.2504    0.2209    0.7142    0.2188
min      0.0000    0.0210    0.0014   -3.1331    0.0020
25%      0.0000    0.3926    0.5023   -0.5414    0.1757
50%      1.0000    0.6137    0.6798    0.1664    0.3162
75%      1.0000    0.7983    0.8230    0.5137    0.4941
max      1.0000    0.9980    0.9980    2.0258    0.9969

LPML (log pseudo-marginal likelihood) = -214.75

Planted outliers in bottom-25 CPO:      13 / 25
Planted outliers in top-25 |Pearson|:   13 / 25
Planted outliers in top-25 |e_bar|:     13 / 25
In [6]:
# ── Plot 1: CPO ranked chart ──────────────────────────────────────────────────
order = resid['CPO'].argsort().values        # ascending CPO order
cols  = ['#d62728' if is_outlier[i] else 'steelblue' for i in order]

fig, axes = plt.subplots(1, 2, figsize=(16, 4))

ax = axes[0]
ax.bar(range(n), resid['CPO'].iloc[order], color=cols, width=1.0, alpha=0.85)
threshold_cpo = resid['CPO'].quantile(0.063)   # flag bottom ~25/400
ax.axhline(threshold_cpo, color='orange', ls='--', lw=1.5,
           label=f'Bottom 6% threshold ({threshold_cpo:.4f})')
ax.set_xlabel('Observation rank (ascending CPO)')
ax.set_ylabel('CPO')
ax.set_title('CPO ranked — red bars are planted outliers')
ax.legend(fontsize=8)

ax = axes[1]
ax.bar(range(60), resid['CPO'].iloc[order[:60]], color=cols[:60], width=0.8, alpha=0.85)
ax.axhline(threshold_cpo, color='orange', ls='--', lw=1.5)
for rank, obs_i in enumerate(order[:60]):
    if is_outlier[obs_i]:
        ax.annotate(str(obs_i), (rank, resid['CPO'].iloc[obs_i]),
                    textcoords='offset points', xytext=(0, 3),
                    fontsize=6, ha='center', color='darkred')
ax.set_xlabel('Rank'); ax.set_ylabel('CPO')
ax.set_title('Bottom 60 by CPO (zoomed, obs index labelled)')

fig.suptitle('§2 CPO-based outlier detection — red = planted label-flip', fontsize=11)
plt.tight_layout(); plt.show()
No description has been provided for this image
In [7]:
# ── Plot 2: latent residuals, p_bar, CPO vs Pearson ──────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(16, 4))

# Panel 1: e_bar sorted
e_order = resid['e_bar'].argsort().values
e_cols  = []
for i in e_order:
    if is_outlier[i]:
        e_cols.append('#d62728')
    elif y[i] == 1:
        e_cols.append('steelblue')
    else:
        e_cols.append('#aec7e8')

ax = axes[0]
ax.bar(range(n), resid['e_bar'].iloc[e_order], color=e_cols, width=1.0, alpha=0.85)
ax.axhline(0, color='k', lw=0.8, ls=':')
ax.set_xlabel('Observation rank (by ē_i)')
ax.set_ylabel('Posterior mean latent residual ē_i')
ax.set_title('Latent residuals ē_i = E[z_i − x_i\'β | y]\nred=outlier, blue=y=1, light=y=0')

# Panel 2: p_bar vs y
ax = axes[1]
jitter = np.random.default_rng(7).uniform(-0.04, 0.04, n)
ax.scatter(resid.loc[~is_outlier, 'p_bar'], y[~is_outlier] + jitter[~is_outlier],
           s=15, alpha=0.4, color='steelblue', label='Normal')
ax.scatter(resid.loc[is_outlier, 'p_bar'], y[is_outlier] + jitter[is_outlier],
           s=60, alpha=0.9, color='red', marker='x', lw=1.5, label='Planted outlier')
ax.axvline(0.5, color='k', ls=':', lw=0.8)
ax.set_xlabel('Posterior predictive probability p̄_i')
ax.set_ylabel('Observed y_i (jittered)')
ax.set_yticks([0, 1]); ax.set_yticklabels(['y=0', 'y=1'])
ax.set_title('p̄_i vs y_i — outliers appear on wrong side\n(high p̄ for y=0, low p̄ for y=1)')
ax.legend(fontsize=8)

# Panel 3: CPO vs |Pearson|
ax = axes[2]
ax.scatter(np.abs(resid.loc[~is_outlier, 'pearson']), resid.loc[~is_outlier, 'CPO'],
           s=15, alpha=0.4, color='steelblue', label='Normal')
ax.scatter(np.abs(resid.loc[is_outlier, 'pearson']), resid.loc[is_outlier, 'CPO'],
           s=60, alpha=0.9, color='red', marker='x', lw=1.5, label='Planted outlier')
ax.set_xlabel('|Pearson residual| (MLE-based)')
ax.set_ylabel('CPO (Bayesian)')
ax.set_title('CPO vs |Pearson residual|\nCPO separates outliers more cleanly')
ax.legend(fontsize=8)

fig.suptitle('§2 Bayesian residual diagnostics — Albert & Chib (1995)', fontsize=11)
plt.tight_layout(); plt.show()
No description has been provided for this image
In [8]:
# ── Detection performance: CPO vs Pearson vs |e_bar| ─────────────────────────
print(f'Detection summary — top-k flag vs {n_out} planted outliers')
print('-' * 65)
print(f'{"k flagged":>10s}  {"CPO":>14s}  {"Pearson":>14s}  {"|e_bar|":>14s}')
print(f'{"":>10s}  {"TP  Prec":>14s}  {"TP  Prec":>14s}  {"TP  Prec":>14s}')
print('-' * 65)

for k_flag in [n_out, 2*n_out, 3*n_out]:
    cpo_tp = sum(is_outlier[i] for i in resid['CPO'].nsmallest(k_flag).index)
    prs_tp = sum(is_outlier[i] for i in resid['pearson'].abs().nlargest(k_flag).index)
    ebar_tp= sum(is_outlier[i] for i in resid['e_bar'].abs().nlargest(k_flag).index)
    print(f'  k={k_flag:3d}       '
          f'  {cpo_tp:2d}/{n_out} ({cpo_tp/k_flag:.2f})   '
          f'  {prs_tp:2d}/{n_out} ({prs_tp/k_flag:.2f})   '
          f'  {ebar_tp:2d}/{n_out} ({ebar_tp/k_flag:.2f})')

Detection summary — top-k flag vs 25 planted outliers
-----------------------------------------------------------------
 k flagged             CPO         Pearson         |e_bar|
                  TP  Prec        TP  Prec        TP  Prec
-----------------------------------------------------------------
  k= 25         13/25 (0.52)     13/25 (0.52)     13/25 (0.52)
  k= 50         20/25 (0.40)     20/25 (0.40)     20/25 (0.40)
  k= 75         23/25 (0.31)     23/25 (0.31)     23/25 (0.31)

Interim summary — the sampler and the residual diagnostics (Sections 1–2)¶

Note on what Sections 1–2 actually showed. On this design the three diagnostics tied exactly (CPO, |Pearson| and |ē| each flagged 13/25, 20/25, 23/25 at k = 25/50/75). The reasons below are why CPO should dominate in principle; Section 3 constructs the scenario where it actually does.

Albert & Chib (1993) sampler¶

Two-block Gibbs via data augmentation: augmenting with latent utilities $z_i$ reduces the probit likelihood to a Normal regression, making $\beta$ conjugate given $z$. The posterior precision $\bar{B}^{-1} = X'X + B_0^{-1}$ is constant across iterations (σ²=1 fixed) — the Cholesky factorisation is computed once, making the sampler very fast.

Albert & Chib (1995) residuals¶

The latent draws $z_i^{(r)}$ provide a natural residual diagnostic with no classical analogue:

DiagnosticDefinitionOutlier signal
Latent residual $\bar{e}_i$$E[z_i - x_i'\beta \mid y]$$|\bar{e}_i|$ large — $z_i$ squeezed against boundary
CPO$p(y_i \mid y_{-i})$Low CPO — model predicts $y_i$ poorly without obs $i$
$\bar{p}_i$$E[\Phi(x_i'\beta) \mid y]$High for $y_i=0$ or low for $y_i=1$

CPO outperforms classical Pearson residuals because it:

  1. Accounts for parameter uncertainty (averages over posterior, not plug-in MLE)
  2. Uses a proper probability model ($p(y_i \mid y_{-i})$ is a probability, not a distance)
  3. Captures influence — how much observation $i$ would shift the posterior if removed

Gibbs vs MLE¶

For correctly specified probit with moderate $n$, MLE and Gibbs posterior means agree closely. The Bayesian advantage comes from: natural uncertainty quantification, direct access to $z_i$ draws for residual diagnostics, and robustness via prior when $n$ is small.

3. When CPO outperforms Pearson: influential masking with a strong prior¶

In Sections 1–2 all three diagnostics tied because:

  1. $n=400$ is large → posterior concentrates near the MLE → CPO ≈ plug-in
  2. Random label-flips were near the boundary → influential outliers with extreme $x$ were not planted

This section designs the scenario where CPO genuinely wins.

Mechanism¶

$$\text{CPO}_i = \left(\frac{1}{R}\sum_r \frac{1}{p(y_i \mid \beta^{(r)})}\right)^{-1}$$

For an influential outlier at $x_i=1.4$, $y_i=0$ (true $P(y=1)\approx0.999$):

  • The harmonic mean is dominated by draws where $p(y_i=0\mid\beta^{(r)})$ is near zero.
  • Such draws exist when $\beta_1^{(r)}$ is close to the true value (2.5): $p(y=0\mid\beta_1=2.5, x=1.4) = 1-\Phi(3.5)\approx 0.0002$.
  • If the prior keeps enough mass near $\beta_1=2.5$, those draws exist → CPO is very low.

Pearson uses only the MLE, which is completely contaminated by the outliers. With 5 influential outliers constituting 17% of $n=30$, the MLE shrinks $\hat{\beta}_1$ from 2.5 to ≈0.46. Under the contaminated MLE the outliers are assigned moderate predicted probabilities of only ≈0.56–0.60, so their |Pearson| residuals (0.56–0.60) fall below those of several perfectly clean observations (which reach ≈0.71). Pearson therefore ranks clean points as the worst offenders — the outliers are masked.

Design¶

  • $n=30$: $n_\text{clean}=25$, $n_\text{inf}=5$ outliers at $x_1 \in \{1.3, 1.4, 1.5\}$, $y=0$
  • $k=2$: intercept + one predictor, $\beta^\star=[0, 2.5]$
  • Prior: $\beta_1 \sim N(2.5,\; 0.25^2)$ — encodes domain knowledge that the slope is ≈2.5
  • Pearson computed from the unconstrained MLE (no prior)
  • Expected result: CPO detects 4–5/5 outliers; Pearson detects 0–2/5
In [9]:
# ── §3 data ───────────────────────────────────────────────────────────────────
rng3 = np.random.default_rng(7)   # seed=7 gives the clearest CPO vs Pearson contrast

k3         = 2
beta_true3 = np.array([0.0, 2.5])

# clean observations
n_clean3 = 25
x_clean3 = rng3.standard_normal(n_clean3)
Xc3      = np.column_stack([np.ones(n_clean3), x_clean3])
p_clean3 = ndtr(Xc3 @ beta_true3)
y_clean3 = (rng3.uniform(size=n_clean3) < p_clean3).astype(float)

# 5 influential outliers: x in [1.3, 1.5], y=0 (true P(y=1) ~ 0.996-0.999)
# These constitute 17% of the data and pull beta_1 from 2.5 toward 0.46 (MLE)
x_inf3 = np.array([1.3, 1.4, 1.5, 1.4, 1.3])
Xi3    = np.column_stack([np.ones(len(x_inf3)), x_inf3])
y_inf3 = np.zeros(len(x_inf3))    # FLIPPED: true P(y=1) ~ 0.996

X3      = np.vstack([Xc3, Xi3])
y3      = np.hstack([y_clean3, y_inf3])
n3      = len(y3)
is_inf3 = np.array([False]*n_clean3 + [True]*len(x_inf3))

print(f'n={n3}  (clean={n_clean3}, influential={is_inf3.sum()})')
print(f'True beta = {beta_true3}')
print(f'\nInfluential outlier x values:  {x_inf3}')
print(f'True P(y=1) at those x:        {ndtr(beta_true3[1]*x_inf3).round(4)}')
print(f'Observed y (all flipped to 0): {y_inf3.astype(int)}')
print(f'\nInfluential outliers = {is_inf3.sum()}/{n3} = {is_inf3.mean():.0%} of data')

n=30  (clean=25, influential=5)
True beta = [0.  2.5]

Influential outlier x values:  [1.3 1.4 1.5 1.4 1.3]
True P(y=1) at those x:        [0.9994 0.9998 0.9999 0.9998 0.9994]
Observed y (all flipped to 0): [0 0 0 0 0]

Influential outliers = 5/30 = 17% of data
In [10]:
# ── §3 Gibbs (informative prior) + MLE ───────────────────────────────────────
# Prior: beta_1 ~ N(2.5, 0.25^2)
# Encodes domain knowledge that the slope is around 2.5.
# MLE uses no prior — it is pulled by the influential outliers.
b0_3 = np.array([0.0, 2.5])
B0_3 = np.diag([100.0, 0.0625])   # flat for intercept, SD=0.25 for beta_1

out3 = binary_probit_gibbs(y3, X3, b0=b0_3, B0=B0_3, R=20000, burn=4000, seed=3, nprint=0)
mle3 = probit_mle(y3, X3)

print('True beta:    ', beta_true3)
print('Gibbs mean:   ', out3['beta'].mean(0).round(4), '  (prior-constrained)')
print('MLE:          ', mle3['beta'].round(4), '  (no prior, fully contaminated)')
print()
print(f'beta_1:  true=2.5  |  Gibbs posterior mean={out3["beta"][:,1].mean():.3f}  '
      f'|  MLE={mle3["beta"][1]:.3f}')
print(f'Fraction of Gibbs draws with beta_1 > 2.0: {(out3["beta"][:,1]>2.0).mean():.1%}')
print(f'  -> for those draws: p(y=0 | x=1.4) = {1-ndtr(2.0*1.4):.4f}  (drives CPO down)')

# Bayesian residuals
resid3 = bayesian_residuals(y3, X3, out3['beta'], out3['z'])
resid3['label'] = np.where(is_inf3, 'INF', 'clean')

# Pearson using the UNCONSTRAINED MLE
p_hat3_mle = ndtr(X3 @ mle3['beta'])
resid3['pearson'] = (y3 - p_hat3_mle).round(4)

print('\n--- 5 influential outliers ---')
print(resid3[is_inf3][['label','CPO','pearson','p_bar','e_bar']].to_string())

print('\n--- Top-8 clean obs by |Pearson| (some rank above the outliers!) ---')
clean_by_prs = resid3[~is_inf3].reindex(resid3[~is_inf3]['pearson'].abs().nlargest(8).index)
print(clean_by_prs[['label','CPO','pearson','p_bar','e_bar']].to_string())

print('\nDetection k=5:')
cpo5 = sum(is_inf3[i] for i in resid3['CPO'].nsmallest(5).index)
prs5 = sum(is_inf3[i] for i in resid3['pearson'].abs().nlargest(5).index)
ebar5 = sum(is_inf3[i] for i in resid3['e_bar'].abs().nlargest(5).index)
print(f'  CPO       bottom-5: {cpo5}/5 TP')
print(f'  |Pearson| top-5:    {prs5}/5 TP')
print(f'  |e_bar|   top-5:    {ebar5}/5 TP')

Done in 0.0 min  |  kept draws: 16000
True beta:     [0.  2.5]
Gibbs mean:    [-0.951   1.7835]   (prior-constrained)
MLE:           [-0.4536  0.4626]   (no prior, fully contaminated)

beta_1:  true=2.5  |  Gibbs posterior mean=1.783  |  MLE=0.463
Fraction of Gibbs draws with beta_1 > 2.0: 15.9%
  -> for those draws: p(y=0 | x=1.4) = 0.0026  (drives CPO down)

--- 5 influential outliers ---
   label       CPO  pearson   p_bar   e_bar
25   INF  0.067650  -0.5588  0.9035 -1.8310
26   INF  0.045599  -0.5770  0.9287 -1.9861
27   INF  0.029541  -0.5950  0.9484 -2.1366
28   INF  0.045599  -0.5770  0.9287 -1.9797
29   INF  0.067650  -0.5588  0.9035 -1.8242

--- Top-8 clean obs by |Pearson| (some rank above the outliers!) ---
    label       CPO  pearson   p_bar   e_bar
21  clean  0.066662   0.7131  0.0962  1.8360
14  clean  0.134746   0.6798  0.1686  1.5354
6   clean  0.174887   0.6649  0.2093  1.4069
12  clean  0.197670   0.6572  0.2321  1.3426
23  clean  0.294038   0.6286  0.3266  1.1177
1   clean  0.311719   0.6238  0.3437  1.0888
11  clean  0.350402   0.6135  0.3810  1.0194
10  clean  0.443542   0.5898  0.4703  0.8553

Detection k=5:
  CPO       bottom-5: 4/5 TP
  |Pearson| top-5:    0/5 TP
  |e_bar|   top-5:    4/5 TP
In [11]:
# ── §3 visualisation ─────────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

clean_idx = np.where(~is_inf3)[0]
inf_idx   = np.where(is_inf3)[0]

# ── Panel 1: Posterior of beta_1 — two competing forces
ax = axes[0]
ax.hist(out3['beta'][:, 1], bins=80, density=True,
        alpha=0.65, color='steelblue', label='Gibbs posterior β₁ (prior-constrained)')
ax.axvline(beta_true3[1], color='k', ls='--', lw=2,   label=f'True β₁ = {beta_true3[1]}')
ax.axvline(mle3['beta'][1], color='red', ls='-', lw=2, label=f'MLE β₁ = {mle3["beta"][1]:.2f}')
ax.axvspan(mle3['beta'][1] - 2*mle3['se'][1],
           mle3['beta'][1] + 2*mle3['se'][1],
           alpha=0.12, color='red')
frac20 = (out3['beta'][:, 1] > 2.0).mean()
ax.axvline(2.0, color='purple', ls=':', lw=1.5,
           label=f'β₁>2: {frac20:.1%} of draws\n→ p(y=0|x=1.4)={1-ndtr(2.0*1.4):.4f}')
ax.set_xlabel('β₁'); ax.set_ylabel('density')
ax.set_title('Posterior β₁ (blue) vs contaminated MLE (red)\n'
             'Prior keeps draws near true 2.5; MLE collapses to 0.46')
ax.legend(fontsize=7.5)

# ── Panel 2: CPO vs |Pearson| — the key comparison
ax = axes[1]
prs_clean = resid3.loc[clean_idx, 'pearson'].abs()
cpo_clean = resid3.loc[clean_idx, 'CPO']
prs_inf   = resid3.loc[inf_idx,   'pearson'].abs()
cpo_inf_v = resid3.loc[inf_idx,   'CPO']

ax.scatter(prs_clean, cpo_clean,
           s=35, alpha=0.7, color='steelblue', label='clean obs', zorder=2)
ax.scatter(prs_inf, cpo_inf_v,
           s=120, alpha=0.9, color='#d62728', marker='^',
           label='influential outlier (y=0, x≈1.4)', zorder=3)

# Shade the "missed by Pearson" region: moderate |Pearson| but low CPO
ax.axvline(resid3['pearson'].abs().nlargest(5).iloc[-1],
           color='orange', ls='--', lw=1.5, label='Pearson top-5 threshold')
ax.axhline(resid3['CPO'].nsmallest(5).iloc[-1],
           color='purple', ls='--', lw=1.5, label='CPO bottom-5 threshold')

ax.set_xlabel('|Pearson residual| — uses contaminated MLE')
ax.set_ylabel('CPO — uses prior-constrained posterior')
ax.set_title('CPO vs |Pearson|\nOutliers (▲): missed by Pearson, caught by CPO')
ax.legend(fontsize=7.5)

# ── Panel 3: ranked CPO and |Pearson| side by side
order_cpo = resid3['CPO'].argsort().values
order_prs = resid3['pearson'].abs().argsort().values[::-1]

ax = axes[2]
rank_cpo = {obs: r for r, obs in enumerate(order_cpo)}
rank_prs = {obs: r for r, obs in enumerate(order_prs)}

# Draw scatter: rank by CPO (x) vs rank by Pearson (y); color by type
x_cpo = [rank_cpo[i] for i in range(n3)]
y_prs = [rank_prs[i] for i in range(n3)]
colors = ['#d62728' if is_inf3[i] else 'steelblue' for i in range(n3)]
ax.scatter(x_cpo, y_prs, c=colors, s=40, alpha=0.8, zorder=2)

ax.axhline(4.5, color='purple', ls='--', lw=1.2, label='Pearson top-5 cutoff (rank 5)')
ax.axvline(4.5, color='orange', ls='--', lw=1.2, label='CPO bottom-5 cutoff (rank 5)')
ax.fill_between([0, 4.5], [0, 0], [4.5, 4.5], alpha=0.08, color='orange',
                label='Caught by both')
ax.fill_between([0, 4.5], [4.5, 4.5], [n3, n3], alpha=0.08, color='purple',
                label='Caught by CPO only (●)')

ax.set_xlabel('Rank by CPO (0 = most surprising)')
ax.set_ylabel('Rank by |Pearson| (0 = largest)')
ax.set_title('CPO rank vs |Pearson| rank\nRed=outlier; purple zone=caught by CPO only')
ax.legend(fontsize=7.5)

fig.suptitle(
    f'§3  n={n3}, 5 influential outliers at x∈[1.3,1.5], y=0  |  '
    f'Prior β₁~N(2.5, 0.25²)\n'
    f'MLE β₁={mle3["beta"][1]:.2f} (contaminated)  vs  Gibbs mean β₁={out3["beta"][:,1].mean():.2f} (prior-corrected)\n'
    f'CPO detection: {sum(is_inf3[i] for i in resid3["CPO"].nsmallest(5).index)}/5'
    f'   |Pearson| detection: {sum(is_inf3[i] for i in resid3["pearson"].abs().nlargest(5).index)}/5',
    fontsize=10)
plt.tight_layout(); plt.show()
No description has been provided for this image

3. Results — CPO catches exactly what Pearson masks¶

The contaminated MLE. The five influential outliers (17% of $n=30$, all at $x\approx1.3$–$1.5$ with $y=0$ where the true $P(y=1)\approx0.999$) drag the MLE slope from its true 2.5 down to 0.46 — the point estimate is effectively destroyed. The prior $\beta_1\sim N(2.5,\,0.25^2)$ holds the Gibbs posterior at 1.78, with ~16% of draws still above 2.0.

Left panel. The posterior of $\beta_1$ (blue) against the contaminated MLE (red). The prior keeps mass near the truth — and it is precisely those draws (where $p(y=0\mid x=1.4)\approx0.003$) that make the harmonic-mean CPO of the outliers collapse. This is the mechanism: CPO averages over the posterior, so a minority of draws that find the outlier very unlikely dominate the harmonic mean.

Centre / right panels — the masking. Under the contaminated MLE the outliers get predicted probabilities of only ≈0.56–0.60, so their $|$Pearson$|$ residuals (0.56–0.60) sit below those of several perfectly clean observations (up to 0.71). Pearson ranks clean points as the worst offenders.

The verdict (top-5 flagged out of 5 planted):

DiagnosticTrue positives
CPO (bottom 5)4 / 5
latent residual $|\bar e_i|$ (top 5)4 / 5
$|$Pearson$|$ (top 5)0 / 5

Pearson finds none of the five; CPO and the latent residual each find four.

Why Sections 1–2 could not show this. There, all three diagnostics tied exactly (13/25, 20/25, 23/25), because $n=400$ makes the posterior concentrate on the MLE (so CPO ≈ the plug-in) and the random label-flips sat near the decision boundary, where they have little leverage. The Bayesian diagnostic only separates when the outliers are influential — high-leverage enough to contaminate the very estimate that Pearson depends on. That is the honest scope of the claim: CPO is not universally better, it is better exactly when it matters.