Skip to content

Paired Model — Bayesian Bootstrap

Overview

The paired Bayesian-bootstrap test (Rubin, 1981) is a nonparametric alternative to PairedBayesPropTest(method="laplace") and PairedBayesPropTest(method="pg"). Instead of specifying a parametric likelihood or a prior on a logit-scale effect, it places a flat Dirichlet "prior" over the empirical distribution of paired differences and reads the posterior on the average treatment effect off the resulting draws.

Use it when you want to:

  • sidestep prior elicitation entirely;
  • protect against model misspecification of the paired logistic likelihood;
  • generate a posterior on \(\Delta = p_A - p_B\) with no latent \(\delta_A\) on the logit scale.

The price you pay is that no Savage–Dickey BF is available — there is no parametric prior on \(\Delta\) to evaluate at the null. Decisions are routed through the ROPE / posterior-mass framework, and the decision surface is intentionally lean: three quantities are enough.

Generative model

For paired binary observations \((y_{A,i}, y_{B,i})\) form the per-pair differences

\[ D_i = y_{A,i} - y_{B,i} \in \{-1, 0, +1\}. \]

Each posterior draw of the average treatment effect is

\[ \Delta^{(s)} = \sum_{i=1}^n w_i^{(s)} D_i, \qquad \mathbf{w}^{(s)} \sim \text{Dirichlet}(\alpha, \dots, \alpha) \]

with \(\alpha = 1\) the standard noninformative choice. The class exposes dirichlet_alpha as a configuration knob: values \(< 1\) concentrate posterior mass on a small number of observations (sharper, more bootstrap-like); values \(> 1\) smooth toward the empirical mean.

Quick start

import numpy as np
from bayesprop.resources.bayes_paired import PairedBayesPropTest

# Paired binary outcomes (any 0/1 arrays of equal length).
y_A = np.array([1, 1, 0, 1, 1, 0, 1, 1, 1, 0])
y_B = np.array([0, 1, 0, 0, 1, 0, 0, 1, 0, 0])

model = PairedBayesPropTest(method="bootstrap", n_samples=20_000, seed=42).fit(y_A, y_B)

s = model.summary
print(f"Posterior mean Δ = {s.mean_delta:+.4f}")
print(f"95% CI = [{s.ci_95.lower:.4f}, {s.ci_95.upper:.4f}]")

Decision API

The class ships only the three quantities that are well-defined directly under the BB posterior — no Bayes-factor machinery, no synthetic prior on \(H_0\):

Quantity Where to read it
Posterior of null — \(P(\Delta \in \text{ROPE} \mid \text{data})\) model.rope_test().pct_in_rope
Posterior of superiority — \(P(p_A > p_B \mid \text{data})\) model.summary.p_A_greater_B
ROPE decision (reject / accept / undecided) model.rope_test().decision 
# (1) Posterior of null + full ROPE result
r = model.rope_test(rope=(-0.05, 0.05))
print(f"P(Δ ∈ ROPE) = {r.pct_in_rope:.3f}  → {r.decision}")
print(f"95% CI for Δ = [{r.ci_lower:.3f}, {r.ci_upper:.3f}]")

# (2) Posterior of superiority — straight off the fitted summary
print(f"P(p_A > p_B | data) = {model.summary.p_A_greater_B:.3f}")

# (3) Composite decision (bayes_factor and posterior_null are None by design)
d = model.decide()
assert d.bayes_factor is None
assert d.posterior_null is None
assert d.rule == "rope"
print(d.rope.decision)

What is intentionally not exposed:

  • savage_dickey_test() — no parametric prior on \(\Delta\) to evaluate at the null.
  • posterior_probability_H0() — under the BB this is just rope_test().pct_in_rope read off the posterior directly. Wrapping it would force the user to commit to a prior on \(H_0\) that has no role in the BB posterior itself, and any default flat-prior choice would be reparametrisation-non-invariant (Lindley–Jeffreys).

If you want a prior-dependent posterior probability of \(H_0\) that does react to evidence in a Bayes-factor sense, use a parametric paired model (Laplace or Pólya–Gamma) with a Savage–Dickey BF — see Paired Laplace and Paired Pólya–Gamma.

When to prefer this over the parametric paired classes

Question Use
Need a Savage–Dickey BF for a point null PairedBayesPropTest(method="laplace") or PairedBayesPropTest(method="pg")
Need sequential / early-stopping support SequentialPairedBayesPropTest
Worried about likelihood misspecification PairedBayesPropTest(method="bootstrap")
Sample size ≤ 30 and prior elicitation is acceptable PairedBayesPropTest(method="pg")
Sample size ≫ 100 and want a prior-free posterior on \(\Delta\) PairedBayesPropTest(method="bootstrap")
Want frequentist OC analysis with a McNemar baseline PairedBayesPropTest(method="laplace") (see Frequentist Evaluation — Paired Laplace)

Plotting

model.plot_posteriors()                    # θ_A and θ_B overlay
model.plot_posterior_delta()                # Δ = θ_A − θ_B with 95 % CI

plot_posteriors() shows the KDE overlay of the marginal θ_A and θ_B posteriors. plot_posterior_delta() shows the posterior of Δ = θ_A − θ_B on the probability scale with the 95 % CI band.

Performance notes

The implementation is fully vectorised — a single rng.dirichlet(α·1_n, size=S) produces all \(S\) weight vectors at once, followed by one matmul W @ D for the posterior draws. On the notebook's n=10, S=50 000 example this is ~200× faster than a per-draw Python loop. For large \(n\) the weight matrix is chunked internally to keep peak memory below ~400 MB.

References

  1. Rubin (1981). The Bayesian Bootstrap. The Annals of Statistics, 9(1), 130–134.
  2. Kruschke (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270–280. (ROPE-based decision making.)

Inputs and binarisation

PairedBayesPropTest(method="bootstrap") accepts both already-binary {0, 1} inputs and continuous scores in [0, 1]. Continuous inputs are auto-binarised at a configurable threshold (default 0.5):

model = PairedBayesPropTest(method="bootstrap", threshold=0.5, verbose=True).fit(scores_A, scores_B)

Values strictly outside [0, 1] or NaN raise ValueError instead of being silently truncated.

API

See API Reference — Paired Model (Bayesian Bootstrap) for full method-level documentation.