Paired Model — Laplace Approximation

Paired logistic model with MAP + Hessian (Laplace) posterior inference.

Supports two modes:

• Fixed priors (default) — \(\mu \sim \mathcal{N}(0, \sigma_\mu)\), \(\delta_A \sim \mathcal{N}(0, \sigma_\delta)\) with user-chosen scales. 2-D Newton solver.
• Hierarchical (learned scales) — places Inverse-Gamma hyperpriors on \(\sigma_\mu^2\) and \(\sigma_\delta^2\) so the prior widths are learned from data. 4-D Newton solver; the marginal prior on \(\delta_A\) becomes Student-\(t\), making the Savage–Dickey Bayes factor robust to prior misspecification. Activated by passing hyperprior_mu and hyperprior_delta to the constructor.

The sequential variant SequentialPairedBayesPropTest is documented on the Sequential designs page.

bayes_paired_laplace

Pooled Bernoulli logistic regression for paired A/B model comparison (Laplace).

This module provides :class:PairedBayesPropTest, a self-contained class for fitting a pooled Bayesian Bernoulli logistic model to paired binary scores via Laplace approximation (MAP + analytical Hessian), performing hypothesis testing via the Savage-Dickey density ratio, running posterior-predictive diagnostics, and generating publication-ready plots.

For exact MCMC inference via Pólya-Gamma data augmentation, see :mod:bayesprop.resources.bayes_paired_pg.

Typical workflow::

from bayesprop.resources.bayes_paired import PairedBayesPropTest

model = PairedBayesPropTest(seed=42).fit(y_A, y_B)
model.print_summary()
model.plot_posterior_delta()
model.plot_savage_dickey()

For multi-metric comparisons::

results = {"Relevancy": model_rel, "Faithfulness": model_faith}
PairedBayesPropTest.plot_forest(results, label_A="v2", label_B="v1")

PairedBayesPropTestLaplace = _PairedLaplace module-attribute

Explicit alias for the Laplace backend class.

_PairedLaplace(prior_sigma_delta=1.0, seed=0, n_samples=8000, decision_rule='all', rope_epsilon=0.02, threshold=0.5, verbose=False, hyperprior_mu=None, hyperprior_delta=None)

Bases: BaseBayesPropTest

Pooled Bernoulli logistic model for paired A/B comparison (Laplace backend).

Uses Laplace approximation (MAP + Hessian) instead of full MCMC for fast, analytic posterior inference on binarized scores.

Generative model (fixed priors, default)::

μ      ~ N(0, 2)              (overall intercept)
δ_A    ~ N(0, σ_δ)            (model-A advantage)
y_A,i  ~ Bernoulli(σ(μ + δ_A))
y_B,i  ~ Bernoulli(σ(μ))

Hierarchical variant (when hyperprior_mu and hyperprior_delta are supplied)::

σ²_μ   ~ IG(a_μ, b_μ)
σ²_δ   ~ IG(a_δ, b_δ)
μ      ~ N(0, σ²_μ)
δ_A    ~ N(0, σ²_δ)
y_A,i  ~ Bernoulli(σ(μ + δ_A))
y_B,i  ~ Bernoulli(σ(μ))

In the hierarchical case the prior scales are learned from the data via joint MAP optimisation over (μ, δ_A, log σ_μ, log σ_δ); the Laplace covariance on (μ, δ_A) is the marginal of the 4-D inverse Hessian and therefore already accounts for hyperparameter uncertainty.

Only 2 regression parameters — no confounding between item effects and model effect when >90% of pairs are concordant.

The prior width on delta_A is fixed (not learned) so the Savage-Dickey density ratio remains exactly consistent. In the hierarchical variant the marginal prior on δ_A is a scaled Student-t (integrated over the IG hyperprior).

Attributes:

Name	Type	Description
laplace	dict[str, Any] | None	Dict with MAP estimate, covariance, Hessian, and posterior samples (None before :meth:fit).
summary	dict[str, Any] | None	Dict with mean_delta, ci_95, P(A > B), and delta_A_posterior_mean on the probability scale.
trace_summary	DataFrame | None	pandas.DataFrame with posterior summary for delta_A and mu.
delta_A_samples	ndarray | None	1-D array of posterior draws for delta_A (logit scale), shape (n_samples,).
y_A_obs	ndarray | None	Observed binary scores for model A (set by :meth:fit).
y_B_obs	ndarray | None	Observed binary scores for model B (set by :meth:fit).

Initialise model configuration.

Parameters:

Name	Type	Description	Default
prior_sigma_delta	float	Standard deviation of the N(0, σ) prior on delta_A (logit scale). When hyperpriors are active this serves as the Newton warm-start for σ_δ.	1.0
seed	int	Random seed for reproducibility.	0
n_samples	int	Number of draws from the Laplace posterior.	8000
decision_rule	DecisionRuleType	Default decision framework — one of "bayes_factor", "posterior_null", "rope", or "all".	'all'
rope_epsilon	float	Half-width of the ROPE interval (default 0.02 = 2 pp).	0.02
threshold	float	Cutoff used to binarise continuous inputs in [0, 1] passed to :meth:fit. Already-binary inputs are left untouched. Defaults to 0.5.	0.5
verbose	bool	If True, emit a one-line notice whenever continuous inputs are binarised.	False
hyperprior_mu	tuple[float, float] | None	(a, b) shape and scale of an Inverse-Gamma hyperprior on σ²_μ. When set together with hyperprior_delta the model becomes hierarchical and both prior scales are learned from the data. None (default) keeps σ_μ fixed at 2.0.	None
hyperprior_delta	tuple[float, float] | None	(a, b) shape and scale of an Inverse-Gamma hyperprior on σ²_δ. None (default) keeps σ_δ fixed at prior_sigma_delta.	None

Source code in bayesprop/resources/bayes_paired_laplace.py

def __init__(
    self,
    prior_sigma_delta: float = 1.0,
    seed: int = 0,
    n_samples: int = 8000,
    decision_rule: DecisionRuleType = "all",
    rope_epsilon: float = 0.02,
    threshold: float = 0.5,
    verbose: bool = False,
    hyperprior_mu: tuple[float, float] | None = None,
    hyperprior_delta: tuple[float, float] | None = None,
) -> None:
    """Initialise model configuration.

    Args:
        prior_sigma_delta: Standard deviation of the N(0, σ) prior on
            ``delta_A`` (logit scale).  When hyperpriors are active
            this serves as the Newton warm-start for σ_δ.
        seed: Random seed for reproducibility.
        n_samples: Number of draws from the Laplace posterior.
        decision_rule: Default decision framework — one of
            ``"bayes_factor"``, ``"posterior_null"``, ``"rope"``, or ``"all"``.
        rope_epsilon: Half-width of the ROPE interval (default 0.02 = 2 pp).
        threshold: Cutoff used to binarise continuous inputs in
            ``[0, 1]`` passed to :meth:`fit`. Already-binary inputs are
            left untouched. Defaults to ``0.5``.
        verbose: If ``True``, emit a one-line notice whenever
            continuous inputs are binarised.
        hyperprior_mu: ``(a, b)`` shape and scale of an
            Inverse-Gamma hyperprior on σ²_μ.  When set together
            with *hyperprior_delta* the model becomes hierarchical
            and both prior scales are learned from the data.
            ``None`` (default) keeps σ_μ fixed at 2.0.
        hyperprior_delta: ``(a, b)`` shape and scale of an
            Inverse-Gamma hyperprior on σ²_δ.  ``None`` (default)
            keeps σ_δ fixed at *prior_sigma_delta*.
    """
    if (hyperprior_mu is None) != (hyperprior_delta is None):
        raise ValueError(
            "hyperprior_mu and hyperprior_delta must both be set or both be None."
        )

    self.prior_sigma_delta: float = prior_sigma_delta
    self.seed: int = seed
    self.n_samples: int = n_samples
    self.decision_rule: DecisionRuleType = decision_rule
    self.rope_epsilon: float = rope_epsilon
    self.threshold: float = threshold
    self.verbose: bool = verbose
    self.hyperprior_mu: tuple[float, float] | None = hyperprior_mu
    self.hyperprior_delta: tuple[float, float] | None = hyperprior_delta

    # --- Populated by .fit() ---
    self.laplace: dict[str, Any] | None = None
    self.summary: dict[str, Any] | None = None
    self.trace_summary: pd.DataFrame | None = None
    self.delta_A_samples: np.ndarray | None = None
    self.delta_samples: np.ndarray | None = None
    self.y_A_obs: np.ndarray | None = None
    self.y_B_obs: np.ndarray | None = None

__repr__()

Return an informative string representation.

Source code in bayesprop/resources/bayes_paired_laplace.py

def __repr__(self) -> str:
    """Return an informative string representation."""
    cls = type(self).__name__
    header = (
        f"{cls}(n_samples={self.n_samples}, "
        f"prior_\u03c3_\u03b4={self.prior_sigma_delta}, seed={self.seed})"
    )
    if self.summary is None:
        return header
    s = self.summary
    return (
        f"{header}\n"
        f"  \u03b8_A = {s.theta_A_mean:.4f},  \u03b8_B = {s.theta_B_mean:.4f}\n"
        f"  Mean \u0394 = {s.mean_delta:+.4f},  "
        f"95% CI = [{s.ci_95.lower:.4f}, {s.ci_95.upper:.4f}]\n"
        f"  P(A > B) = {s.p_A_greater_B:.4f}"
    )

fit(y_A_obs, y_B_obs)

Fit the pooled Bernoulli model via Laplace approximation.

Reduces (y_A_obs, y_B_obs) to the four sufficient statistics (n_A, k_A, n_B, k_B) and delegates to :func:_paired_laplace_from_counts (fixed priors) or :func:_hierarchical_laplace_from_counts (learned priors).

Fixed-prior mode (default): solves for the 2-D MAP of (μ, δ_A) via damped Newton with closed-form gradient and Hessian, returning the Laplace covariance Σ = H⁻¹.

Hierarchical mode (when hyperprior_mu and hyperprior_delta are set): places Inverse-Gamma hyperpriors on the prior variances and jointly optimises over (μ, δ_A, log σ_μ, log σ_δ) via 4-D Newton. The returned 2×2 marginal Laplace covariance on (μ, δ_A) is the top-left block of the 4×4 H⁻¹ and therefore already incorporates hyperparameter uncertainty.

Log-posterior (up to constant, fixed-prior case)::

log p(μ, δ|y) = Σᵢ [y_Aᵢ log σ(μ+δ) + (1-y_Aᵢ) log(1-σ(μ+δ))]
              + Σᵢ [y_Bᵢ log σ(μ)   + (1-y_Bᵢ) log(1-σ(μ))]
              - μ²/(2σ_μ²) - δ²/(2σ_δ²)

Gradient::

∂/∂μ = (k_A - n_A·p_A) + (k_B - n_B·p_B) - μ/σ_μ²
∂/∂δ = (k_A - n_A·p_A)                    - δ/σ_δ²

Hessian of the negative log-posterior (observed information)::

H[0,0] = n_A·w_A + n_B·w_B + 1/σ_μ²
H[1,1] = n_A·w_A           + 1/σ_δ²
H[0,1] = H[1,0] = n_A·w_A

where w_A = p_A(1-p_A), w_B = p_B(1-p_B), evaluated at MAP.

Parameters:

Name	Type	Description	Default
y_A_obs	ndarray	Observed scores for model A — either binary {0, 1} or continuous in [0, 1]. Continuous inputs are binarised at self.threshold (default 0.5); values outside [0, 1] raise :class:ValueError.	required
y_B_obs	ndarray	Observed scores for model B — same conventions.	required

Returns:

Type	Description
_PairedLaplace	self (for method chaining).

Source code in bayesprop/resources/bayes_paired_laplace.py

def fit(self, y_A_obs: np.ndarray, y_B_obs: np.ndarray) -> _PairedLaplace:
    """Fit the pooled Bernoulli model via Laplace approximation.

    Reduces ``(y_A_obs, y_B_obs)`` to the four sufficient statistics
    ``(n_A, k_A, n_B, k_B)`` and delegates to
    :func:`_paired_laplace_from_counts` (fixed priors) or
    :func:`_hierarchical_laplace_from_counts` (learned priors).

    **Fixed-prior mode** (default): solves for the 2-D MAP of
    ``(μ, δ_A)`` via damped Newton with closed-form gradient
    and Hessian, returning the Laplace covariance ``Σ = H⁻¹``.

    **Hierarchical mode** (when *hyperprior_mu* and
    *hyperprior_delta* are set): places Inverse-Gamma hyperpriors
    on the prior variances and jointly optimises over
    ``(μ, δ_A, log σ_μ, log σ_δ)`` via 4-D Newton.  The
    returned 2×2 marginal Laplace covariance on ``(μ, δ_A)`` is
    the top-left block of the 4×4 ``H⁻¹`` and therefore already
    incorporates hyperparameter uncertainty.

    Log-posterior (up to constant, fixed-prior case)::

        log p(μ, δ|y) = Σᵢ [y_Aᵢ log σ(μ+δ) + (1-y_Aᵢ) log(1-σ(μ+δ))]
                      + Σᵢ [y_Bᵢ log σ(μ)   + (1-y_Bᵢ) log(1-σ(μ))]
                      - μ²/(2σ_μ²) - δ²/(2σ_δ²)

    Gradient::

        ∂/∂μ = (k_A - n_A·p_A) + (k_B - n_B·p_B) - μ/σ_μ²
        ∂/∂δ = (k_A - n_A·p_A)                    - δ/σ_δ²

    Hessian of the *negative* log-posterior (observed information)::

        H[0,0] = n_A·w_A + n_B·w_B + 1/σ_μ²
        H[1,1] = n_A·w_A           + 1/σ_δ²
        H[0,1] = H[1,0] = n_A·w_A

    where w_A = p_A(1-p_A), w_B = p_B(1-p_B), evaluated at MAP.

    Args:
        y_A_obs: Observed scores for model A — either binary ``{0, 1}``
            or continuous in ``[0, 1]``. Continuous inputs are
            binarised at ``self.threshold`` (default ``0.5``); values
            outside ``[0, 1]`` raise :class:`ValueError`.
        y_B_obs: Observed scores for model B — same conventions.

    Returns:
        self (for method chaining).
    """
    y_A_bin = binarize_if_needed(
        y_A_obs, self.threshold, name="y_A_obs", verbose=self.verbose
    )
    y_B_bin = binarize_if_needed(
        y_B_obs, self.threshold, name="y_B_obs", verbose=self.verbose
    )
    self.y_A_obs = y_A_bin.astype(int)
    self.y_B_obs = y_B_bin.astype(int)

    n_A = int(len(self.y_A_obs))
    k_A = int(self.y_A_obs.sum())
    n_B = int(len(self.y_B_obs))
    k_B = int(self.y_B_obs.sum())

    # Closed-form MAP + Hessian directly from sufficient statistics.
    if self.hyperprior_mu is not None:
        theta_map, cov, H, sigma_mu_map, sigma_delta_map = (
            _hierarchical_laplace_from_counts(
                n_A=n_A,
                k_A=k_A,
                n_B=n_B,
                k_B=k_B,
                hp_mu=self.hyperprior_mu,
                hp_delta=self.hyperprior_delta,
                x0=(
                    0.0,
                    0.0,
                    np.log(2.0),
                    np.log(self.prior_sigma_delta),
                ),
            )
        )
    else:
        theta_map, cov, H = _paired_laplace_from_counts(
            n_A=n_A,
            k_A=k_A,
            n_B=n_B,
            k_B=k_B,
            prior_sigma_delta=self.prior_sigma_delta,
        )
        sigma_mu_map = None
        sigma_delta_map = None
    mu_map, delta_map = float(theta_map[0]), float(theta_map[1])

    # Sample from Gaussian posterior
    rng = np.random.default_rng(self.seed)
    samples = rng.multivariate_normal(theta_map, cov, size=self.n_samples)
    mu_s = samples[:, 0]
    delta_A_s = samples[:, 1]

    pA_s = sigmoid(mu_s + delta_A_s)
    pB_s = sigmoid(mu_s)
    Delta_s = pA_s - pB_s

    self.delta_A_samples = delta_A_s
    self.delta_samples = Delta_s

    # Batch all quantile calls per array — np.quantile shares a single
    # internal partition() across requested probabilities, so this cuts
    # six independent O(n log n) sorts down to three.
    delta_lo, delta_hi = np.quantile(Delta_s, [0.025, 0.975])
    dA_lo, dA_hi = np.quantile(delta_A_s, [0.03, 0.97])
    mu_lo, mu_hi = np.quantile(mu_s, [0.03, 0.97])

    self.summary = PairedSummary(
        mean_delta=float(Delta_s.mean()),
        ci_95=CredibleInterval(lower=float(delta_lo), upper=float(delta_hi)),
        **{"P(A > B)": float((Delta_s > 0).mean())},
        theta_A_mean=float(pA_s.mean()),
        theta_B_mean=float(pB_s.mean()),
        delta_A_posterior_mean=float(delta_A_s.mean()),
    )

    self.trace_summary = pd.DataFrame(
        {
            "mean": [delta_A_s.mean(), mu_s.mean()],
            "sd": [delta_A_s.std(), mu_s.std()],
            "hdi_3%": [dA_lo, mu_lo],
            "hdi_97%": [dA_hi, mu_hi],
            "MAP": [delta_map, mu_map],
        },
        index=["delta_A", "mu"],
    )

    self.laplace = {
        "map": theta_map,
        "cov": cov,
        "H": H,
        "mu_samples": mu_s,
        "delta_A_samples": delta_A_s,
        "n_A": n_A,
        "k_A": k_A,
        "n_B": n_B,
        "k_B": k_B,
        "prior_sigma_delta": self.prior_sigma_delta,
        "hierarchical": self.hyperprior_mu is not None,
        "sigma_mu_map": sigma_mu_map,
        "sigma_delta_map": sigma_delta_map,
    }

    return self

_check_fitted()

Raise RuntimeError if the model has not been fitted yet.

Source code in bayesprop/resources/bayes_paired_laplace.py

def _check_fitted(self) -> None:
    """Raise RuntimeError if the model has not been fitted yet."""
    if self.laplace is None:
        raise RuntimeError("Model has not been fitted yet. Call .fit() first.")

savage_dickey_test(null_value=0.0)

Savage-Dickey density-ratio Bayes factor for H0: delta_A = null_value.

Uses the analytical marginal of the Laplace posterior: the joint posterior on (μ, δ_A) is by construction Gaussian with mean theta_map and covariance cov, so the marginal on δ_A is exactly N(δ_A_MAP, cov[1, 1]). We evaluate the density at null_value from this closed form rather than re-estimating it via Gaussian-KDE on samples drawn from the same Gaussian — the closed form is faster, exact, and free of KDE bandwidth-quantisation noise (the latter being especially relevant for the BF tail behaviour at small |null_value − δ_A_MAP| / √cov[1,1]).

The prior on δ_A is N(0, prior_sigma_delta²), so the prior density at null_value is also computed in closed form rather than dispatching through scipy.stats.norm.pdf.

Parameters:

Name	Type	Description	Default
null_value	float	The point null hypothesis value for delta_A.	0.0

Returns:

Type	Description
SavageDickeyResult	class:SavageDickeyResult with BF_01, BF_10, densities,
SavageDickeyResult	interpretation, and decision.

Source code in bayesprop/resources/bayes_paired_laplace.py

def savage_dickey_test(self, null_value: float = 0.0) -> SavageDickeyResult:
    """Savage-Dickey density-ratio Bayes factor for H0: delta_A = *null_value*.

    Uses the *analytical* marginal of the Laplace posterior: the
    joint posterior on ``(μ, δ_A)`` is by construction Gaussian
    with mean ``theta_map`` and covariance ``cov``, so the marginal
    on ``δ_A`` is exactly ``N(δ_A_MAP, cov[1, 1])``. We evaluate
    the density at ``null_value`` from this closed form rather
    than re-estimating it via Gaussian-KDE on samples drawn from
    the same Gaussian — the closed form is faster, exact, and
    free of KDE bandwidth-quantisation noise (the latter being
    especially relevant for the BF tail behaviour at small
    ``|null_value − δ_A_MAP| / √cov[1,1]``).

    The prior on ``δ_A`` is ``N(0, prior_sigma_delta²)``, so the
    prior density at ``null_value`` is also computed in closed
    form rather than dispatching through ``scipy.stats.norm.pdf``.

    Args:
        null_value: The point null hypothesis value for delta_A.

    Returns:
        :class:`SavageDickeyResult` with BF_01, BF_10, densities,
        interpretation, and decision.
    """
    self._check_fitted()

    # Marginal Laplace posterior on δ_A is exactly Gaussian.
    # read off the mean and variance from the joint MAP and Hessian.
    delta_A_map = float(self.laplace["map"][1])
    sigma_post = float(np.sqrt(self.laplace["cov"][1, 1]))
    z = (null_value - delta_A_map) / sigma_post
    posterior_at_null = float(
        np.exp(-0.5 * z * z) / (sigma_post * np.sqrt(2.0 * np.pi))
    )
    if self.laplace["hierarchical"]:
        # Marginal prior on δ_A integrating out σ²_δ ~ IG(a, b)
        # is a scaled Student-t with ν = 2a d.f. and scale √(b/a).
        a_d, b_d = self.hyperprior_delta
        prior_at_null = float(
            student_t.pdf(
                null_value, df=2.0 * a_d, loc=0.0, scale=np.sqrt(b_d / a_d)
            )
        )
    else:
        prior_at_null = float(
            np.exp(-0.5 * (null_value / self.prior_sigma_delta) ** 2)
            / (self.prior_sigma_delta * np.sqrt(2.0 * np.pi))
        )

    BF_01 = posterior_at_null / prior_at_null
    BF_10 = 1.0 / BF_01 if BF_01 > 0 else float("inf")

    if BF_10 > 100:
        interpretation = "Decisive evidence against H0"
    elif BF_10 > 30:
        interpretation = "Very strong evidence against H0"
    elif BF_10 > 10:
        interpretation = "Strong evidence against H0"
    elif BF_10 > 3:
        interpretation = "Moderate evidence against H0"
    elif BF_10 > 1:
        interpretation = "Anecdotal evidence against H0"
    elif BF_10 == 1:
        interpretation = "No evidence either way"
    elif BF_01 > 100:
        interpretation = "Decisive evidence for H0"
    elif BF_01 > 30:
        interpretation = "Very strong evidence for H0"
    elif BF_01 > 10:
        interpretation = "Strong evidence for H0"
    elif BF_01 > 3:
        interpretation = "Moderate evidence for H0"
    else:
        interpretation = "Anecdotal evidence for H0"

    decision = "Reject H0" if BF_10 > 3 else "Fail to reject H0"

    return SavageDickeyResult(
        BF_01=BF_01,
        BF_10=BF_10,
        posterior_density_at_0=posterior_at_null,
        prior_density_at_0=prior_at_null,
        interpretation=interpretation,
        decision=decision,
    )

posterior_probability_H0(BF_01, prior_H0=0.5) staticmethod

Convert BF_01 to posterior probability of H0 (spike-and-slab).

Parameters:

Name	Type	Description	Default
BF_01	float	Bayes factor in favour of H0.	required
prior_H0	float	Prior probability of H0 (default 0.5).	0.5

Returns:

Type	Description
PosteriorProbH0Result	class:PosteriorProbH0Result with posterior and prior odds
PosteriorProbH0Result	and model probabilities.

Source code in bayesprop/resources/bayes_paired_laplace.py

@staticmethod
def posterior_probability_H0(
    BF_01: float, prior_H0: float = 0.5
) -> PosteriorProbH0Result:
    """Convert BF_01 to posterior probability of H0 (spike-and-slab).

    Args:
        BF_01: Bayes factor in favour of H0.
        prior_H0: Prior probability of H0 (default 0.5).

    Returns:
        :class:`PosteriorProbH0Result` with posterior and prior odds
        and model probabilities.
    """
    prior_odds = prior_H0 / (1 - prior_H0)
    posterior_odds = BF_01 * prior_odds
    P_H0 = posterior_odds / (1 + posterior_odds)
    P_H1 = 1 - P_H0

    if P_H1 > 0.95:
        decision = "Reject H0"
    elif P_H0 > 0.95:
        decision = "Fail to reject H0"
    else:
        decision = "Undecided"

    return PosteriorProbH0Result(
        **{"P(H0|data)": P_H0, "P(H1|data)": P_H1},
        prior_odds=prior_odds,
        posterior_odds=posterior_odds,
        decision=decision,
    )

rope_test(rope=None, ci_mass=0.95)

ROPE analysis on the posterior of Δ = p_A − p_B (probability scale).

Parameters:

Name	Type	Description	Default
rope	tuple[float, float] | None	(lower, upper) ROPE bounds. Defaults to (-self.rope_epsilon, +self.rope_epsilon).	None
ci_mass	float	Credible interval mass (default 95%).	0.95

Returns:

Type	Description
ROPEResult	class:ROPEResult with CI, ROPE overlap fractions, and
ROPEResult	decision.

Source code in bayesprop/resources/bayes_paired_laplace.py

def rope_test(
    self,
    rope: tuple[float, float] | None = None,
    ci_mass: float = 0.95,
) -> ROPEResult:
    """ROPE analysis on the posterior of Δ = p_A − p_B (probability scale).

    Args:
        rope: (lower, upper) ROPE bounds. Defaults to
            ``(-self.rope_epsilon, +self.rope_epsilon)``.
        ci_mass: Credible interval mass (default 95%).

    Returns:
        :class:`ROPEResult` with CI, ROPE overlap fractions, and
        decision.
    """
    self._check_fitted()
    if rope is None:
        rope = (-self.rope_epsilon, self.rope_epsilon)
    return ROPEResult.from_samples(self.delta_samples, rope=rope, ci_mass=ci_mass)

decide(rule=None)

Run the chosen decision framework(s) and return a composite result.

Parameters:

Name	Type	Description	Default
rule	DecisionRuleType | None	Override the default decision_rule. One of "bayes_factor", "posterior_null", "rope", or "all".	None

Returns:

Type	Description
HypothesisDecision	class:HypothesisDecision with the requested sub-results
HypothesisDecision	populated.

Source code in bayesprop/resources/bayes_paired_laplace.py

def decide(self, rule: DecisionRuleType | None = None) -> HypothesisDecision:
    """Run the chosen decision framework(s) and return a composite result.

    Args:
        rule: Override the default ``decision_rule``. One of
            ``"bayes_factor"``, ``"posterior_null"``, ``"rope"``,
            or ``"all"``.

    Returns:
        :class:`HypothesisDecision` with the requested sub-results
        populated.
    """
    self._check_fitted()
    rule = rule or self.decision_rule

    bf: SavageDickeyResult | None = None
    pn: PosteriorProbH0Result | None = None
    rp: ROPEResult | None = None

    if rule in ("bayes_factor", "posterior_null", "all"):
        bf = self.savage_dickey_test()
    if rule in ("posterior_null", "all"):
        assert bf is not None  # noqa: S101
        pn = self.posterior_probability_H0(bf.BF_01)
    if rule in ("rope", "all"):
        rp = self.rope_test()

    return HypothesisDecision(
        bayes_factor=bf, posterior_null=pn, rope=rp, rule=rule
    )

ppc_pvalues(seed=None)

Posterior predictive p-values for summary statistics.

Returns:

Type	Description
dict[str, PPCStatistic]	Dict mapping statistic name to :class:PPCStatistic.

Source code in bayesprop/resources/bayes_paired_laplace.py

def ppc_pvalues(self, seed: int | None = None) -> dict[str, PPCStatistic]:
    """Posterior predictive p-values for summary statistics.

    Returns:
        Dict mapping statistic name to :class:`PPCStatistic`.
    """
    self._check_fitted()

    rng = np.random.default_rng(seed if seed is not None else self.seed)

    mu_s = self.laplace["mu_samples"]
    delta_s = self.laplace["delta_A_samples"]
    n = len(self.y_A_obs)

    p_A_s = sigmoid(mu_s + delta_s)
    p_B_s = sigmoid(mu_s)
    y_A_rep = (rng.random((len(p_A_s), n)) < p_A_s[:, None]).astype(int)
    y_B_rep = (rng.random((len(p_B_s), n)) < p_B_s[:, None]).astype(int)

    checks = {
        "mean(y_A)": (self.y_A_obs.mean(), y_A_rep.mean(axis=1)),
        "mean(y_B)": (self.y_B_obs.mean(), y_B_rep.mean(axis=1)),
        "mean(y_A-y_B)": (
            (self.y_A_obs - self.y_B_obs).mean(),
            (y_A_rep - y_B_rep).mean(axis=1),
        ),
        "std(y_A-y_B)": (
            (self.y_A_obs - self.y_B_obs).std(),
            (y_A_rep - y_B_rep).std(axis=1),
        ),
        "n_disagree": (
            float(np.sum(self.y_A_obs != self.y_B_obs)),
            np.sum(y_A_rep != y_B_rep, axis=1).astype(float),
        ),
    }

    results: dict[str, PPCStatistic] = {}
    for stat_name, (obs_val, rep_vals) in checks.items():
        p_val = min(
            2
            * min(
                float((rep_vals >= obs_val).mean()),
                float((rep_vals <= obs_val).mean()),
            ),
            1.0,
        )
        results[stat_name] = PPCStatistic(
            observed=float(obs_val),
            p_value=p_val,
            status="OK" if p_val > 0.05 else "WARN",
        )
    return results

plot_laplace_posterior(**kwargs)

Two-panel posterior plot: overlaid p_A / p_B and Δ = p_A − p_B.

The implied success probabilities p_A = σ(μ + δ_A) and p_B = σ(μ) are computed from the Laplace posterior samples and displayed as overlaid KDE densities in the left panel. The right panel shows the difference Δ = p_A − p_B.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Accepts figsize (default (14, 5)) and title (default "Laplace Posterior (Pooled Binomial)").	{}

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_laplace_posterior(self, **kwargs: Any) -> None:
    """Two-panel posterior plot: overlaid p_A / p_B and Δ = p_A − p_B.

    The implied success probabilities ``p_A = σ(μ + δ_A)`` and
    ``p_B = σ(μ)`` are computed from the Laplace posterior samples
    and displayed as overlaid KDE densities in the left panel.
    The right panel shows the difference Δ = p_A − p_B.

    Args:
        **kwargs: Accepts ``figsize`` (default ``(14, 5)``) and
            ``title`` (default ``"Laplace Posterior (Pooled Binomial)"``).
    """
    import matplotlib.pyplot as plt

    self._check_fitted()
    assert self.laplace is not None
    mu_s = self.laplace["mu_samples"]
    delta_s = self.laplace["delta_A_samples"]

    p_A_s = sigmoid(mu_s + delta_s)
    p_B_s = sigmoid(mu_s)
    Delta_s = p_A_s - p_B_s

    figsize = kwargs.pop("figsize", (14, 5))
    fig, axes = plt.subplots(1, 2, figsize=figsize)

    # Panel 1: p_A and p_B overlaid
    ax = axes[0]
    kde_A = gaussian_kde(p_A_s)
    kde_B = gaussian_kde(p_B_s)
    lo = min(p_A_s.min(), p_B_s.min())
    hi = max(p_A_s.max(), p_B_s.max())
    x = np.linspace(max(0, lo - 0.05), min(1, hi + 0.05), 500)

    pdf_A = kde_A(x)
    pdf_B = kde_B(x)
    ax.plot(
        x,
        pdf_A,
        color="#2196F3",
        linewidth=2,
        label=f"p_A = σ(μ+δ_A)  mean={p_A_s.mean():.3f}",
    )
    ax.fill_between(x, pdf_A, alpha=0.15, color="#2196F3")
    ax.plot(
        x,
        pdf_B,
        color="#4CAF50",
        linewidth=2,
        label=f"p_B = σ(μ)  mean={p_B_s.mean():.3f}",
    )
    ax.fill_between(x, pdf_B, alpha=0.15, color="#4CAF50")

    ax.axvline(
        p_A_s.mean(), color="#2196F3", linestyle="--", linewidth=1, alpha=0.6
    )
    ax.axvline(
        p_B_s.mean(), color="#4CAF50", linestyle="--", linewidth=1, alpha=0.6
    )
    ax.set_xlabel("Success probability")
    ax.set_ylabel("Density")
    ax.set_title("Implied Probability Posteriors", fontsize=11, fontweight="bold")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)

    # Panel 2: Δ = p_A − p_B
    ax = axes[1]
    ax.hist(
        Delta_s,
        bins=60,
        density=True,
        alpha=0.6,
        color="#9C27B0",
        edgecolor="white",
    )
    ax.axvline(0, color="gray", linestyle="--", linewidth=1, alpha=0.6)
    ax.axvline(
        Delta_s.mean(),
        color="#9C27B0",
        linewidth=1.5,
        label=f"Mean = {Delta_s.mean():.4f}",
    )
    ax.set_xlabel("Δ = p_A − p_B")
    ax.set_ylabel("Density")
    ax.set_title("Difference Posterior", fontsize=11, fontweight="bold")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)

    fig.suptitle(
        kwargs.pop(
            "suptitle", kwargs.pop("title", "Laplace Posterior (Pooled Binomial)")
        ),
        fontsize=13,
        fontweight="bold",
        y=1.02,
    )
    plt.tight_layout()
    plt.show()

plot_posteriors(**kwargs)

Overlaid KDE posteriors of θ_A and θ_B (probability scale).

Plots the implied success probabilities θ_A = σ(μ + δ_A) and θ_B = σ(μ) derived from the Laplace posterior samples as overlaid KDE densities.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Accepts figsize (default (7, 5)) and title (default "Posterior: θ_A and θ_B").	{}

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_posteriors(self, **kwargs: Any) -> None:
    """Overlaid KDE posteriors of θ_A and θ_B (probability scale).

    Plots the implied success probabilities ``θ_A = σ(μ + δ_A)``
    and ``θ_B = σ(μ)`` derived from the Laplace posterior samples
    as overlaid KDE densities.

    Args:
        **kwargs: Accepts ``figsize`` (default ``(7, 5)``) and
            ``title`` (default ``"Posterior: θ_A and θ_B"``).
    """
    import matplotlib.pyplot as plt

    self._check_fitted()
    assert self.laplace is not None
    mu_s = self.laplace["mu_samples"]
    delta_s = self.laplace["delta_A_samples"]

    p_A_s = sigmoid(mu_s + delta_s)
    p_B_s = sigmoid(mu_s)

    figsize = kwargs.pop("figsize", (7, 5))
    fig, ax = plt.subplots(figsize=figsize)

    kde_A = gaussian_kde(p_A_s)
    kde_B = gaussian_kde(p_B_s)
    lo = min(p_A_s.min(), p_B_s.min())
    hi = max(p_A_s.max(), p_B_s.max())
    x = np.linspace(max(0, lo - 0.05), min(1, hi + 0.05), 500)

    pdf_A = kde_A(x)
    pdf_B = kde_B(x)
    ax.plot(
        x,
        pdf_A,
        color="#2196F3",
        linewidth=2,
        label=f"θ_A  mean={p_A_s.mean():.3f}",
    )
    ax.fill_between(x, pdf_A, alpha=0.15, color="#2196F3")
    ax.plot(
        x,
        pdf_B,
        color="#4CAF50",
        linewidth=2,
        label=f"θ_B  mean={p_B_s.mean():.3f}",
    )
    ax.fill_between(x, pdf_B, alpha=0.15, color="#4CAF50")

    ax.axvline(
        p_A_s.mean(), color="#2196F3", linestyle="--", linewidth=1, alpha=0.6
    )
    ax.axvline(
        p_B_s.mean(), color="#4CAF50", linestyle="--", linewidth=1, alpha=0.6
    )
    ax.set_xlabel("Success probability")
    ax.set_ylabel("Density")
    ax.set_title(
        kwargs.pop("title", "Posterior: θ_A and θ_B"),
        fontsize=12,
        fontweight="bold",
    )
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
    plt.tight_layout()
    plt.show()

plot_posterior_delta(color='#9C27B0', **kwargs)

KDE posterior density of Δ = θ_A − θ_B (probability scale) with 95% CI.

Parameters:

Name	Type	Description	Default
color	str	Colour for the density curve and fill.	'#9C27B0'
**kwargs	Any	Accepts figsize (default (7, 5)), title (default "Posterior: Δ = θ_A − θ_B"), xlabel, ylabel.	{}

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_posterior_delta(self, color: str = "#9C27B0", **kwargs: Any) -> None:
    """KDE posterior density of Δ = θ_A − θ_B (probability scale) with 95% CI.

    Args:
        color: Colour for the density curve and fill.
        **kwargs: Accepts ``figsize`` (default ``(7, 5)``),
            ``title`` (default ``"Posterior: Δ = θ_A − θ_B"``),
            ``xlabel``, ``ylabel``.
    """
    import matplotlib.pyplot as plt

    self._check_fitted()
    assert self.laplace is not None
    mu_s = self.laplace["mu_samples"]
    delta_s = self.laplace["delta_A_samples"]

    p_A_s = sigmoid(mu_s + delta_s)
    p_B_s = sigmoid(mu_s)
    Delta_s = p_A_s - p_B_s

    ci_low, ci_high = np.quantile(Delta_s, [0.025, 0.975])
    mean_val = float(Delta_s.mean())

    kde = gaussian_kde(Delta_s)
    x_grid = np.linspace(Delta_s.min() - 0.05, Delta_s.max() + 0.05, 500)
    density = kde(x_grid)

    figsize = kwargs.pop("figsize", (7, 5))
    fig, ax = plt.subplots(figsize=figsize)
    ax.plot(x_grid, density, color=color, linewidth=2)
    ax.fill_between(x_grid, density, alpha=0.15, color=color)
    mask = (x_grid >= ci_low) & (x_grid <= ci_high)
    ax.fill_between(
        x_grid[mask], density[mask], alpha=0.35, color=color, label="95% CI"
    )
    ax.axvline(
        mean_val,
        color=color,
        linestyle="-",
        linewidth=1.5,
        alpha=0.8,
        label=f"Mean = {mean_val:.4f}",
    )
    ax.axvline(
        0,
        color="gray",
        linestyle="--",
        linewidth=1,
        alpha=0.6,
        label="Δ = 0 (no difference)",
    )
    ax.set_xlabel(kwargs.pop("xlabel", "Δ = θ_A − θ_B"), fontsize=11)
    ax.set_ylabel(kwargs.pop("ylabel", "Density"), fontsize=11)
    ax.set_title(
        kwargs.pop("title", "Posterior: Δ = θ_A − θ_B"),
        fontsize=12,
        fontweight="bold",
    )
    ax.legend(fontsize=9, loc="upper right")
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()
    plt.show()

plot_savage_dickey(color='#2196F3', **kwargs)

Posterior vs prior density with Savage-Dickey BF annotation.

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_savage_dickey(self, color: str = "#2196F3", **kwargs) -> None:
    """Posterior vs prior density with Savage-Dickey BF annotation."""
    import matplotlib.pyplot as plt

    bf = self.savage_dickey_test()
    samples = self.delta_A_samples

    kde = gaussian_kde(samples)
    x_grid = np.linspace(samples.min() - 0.5, samples.max() + 0.5, 500)
    density = kde(x_grid)
    prior_density = norm.pdf(x_grid, 0, self.prior_sigma_delta)

    figsize = kwargs.pop("figsize", (7, 5))
    fig, ax = plt.subplots(figsize=figsize)
    ax.plot(x_grid, density, color=color, linewidth=2, label="Posterior")
    ax.fill_between(x_grid, density, alpha=0.15, color=color)
    ax.plot(
        x_grid,
        prior_density,
        color="gray",
        linewidth=1.5,
        linestyle="--",
        alpha=0.7,
        label=f"Prior N(0,{self.prior_sigma_delta})",
    )
    ax.plot(
        0,
        bf.posterior_density_at_0,
        "o",
        color="red",
        markersize=10,
        zorder=5,
        label=f"Post. at \u03b4=0: {bf.posterior_density_at_0:.2e}",
    )
    ax.plot(
        0,
        bf.prior_density_at_0,
        "s",
        color="gray",
        markersize=8,
        zorder=5,
        label=f"Prior at \u03b4=0: {bf.prior_density_at_0:.3f}",
    )

    bf10_label = _format_bf(bf.BF_10)
    log10_bf = np.log10(bf.BF_10)
    ax.text(
        0.02,
        0.97,
        f"$BF_{{10}}$ = {bf10_label}\n$\\log_{{10}}BF_{{10}}$ = {log10_bf:.1f}\n{bf.decision}",
        fontsize=10,
        fontweight="bold",
        color="darkred",
        transform=ax.transAxes,
        verticalalignment="top",
        bbox=dict(
            boxstyle="round,pad=0.3",
            facecolor="lightyellow",
            edgecolor="darkred",
            alpha=0.9,
        ),
    )
    ax.set_xlabel(kwargs.pop("xlabel", "\u03b4_A (logit scale)"), fontsize=11)
    ax.set_ylabel(kwargs.pop("ylabel", "Density"), fontsize=11)
    ax.set_title(
        kwargs.pop("title", "Savage-Dickey Test (Binomial)"),
        fontsize=12,
        fontweight="bold",
    )
    ax.legend(fontsize=9, loc="upper right")
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()
    plt.show()

plot_ppc(seed=None, **kwargs)

Three-column PPC plot: P(perfect) A, P(perfect) B, rate difference.

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_ppc(self, seed: int | None = None, **kwargs) -> None:
    """Three-column PPC plot: P(perfect) A, P(perfect) B, rate difference."""
    import matplotlib.pyplot as plt

    self._check_fitted()
    rng = np.random.default_rng(seed if seed is not None else self.seed)

    mu_s = self.laplace["mu_samples"]
    delta_s = self.laplace["delta_A_samples"]
    n = len(self.y_A_obs)

    p_A_s = sigmoid(mu_s + delta_s)
    p_B_s = sigmoid(mu_s)
    y_A_rep = (rng.random((len(p_A_s), n)) < p_A_s[:, None]).astype(int)
    y_B_rep = (rng.random((len(p_B_s), n)) < p_B_s[:, None]).astype(int)

    figsize = kwargs.pop("figsize", (18, 5))
    fig, axes = plt.subplots(1, 3, figsize=figsize)

    # P(perfect) Group A
    ax = axes[0]
    frac_A_rep = y_A_rep.mean(axis=1)
    frac_A_obs = self.y_A_obs.mean()
    ax.hist(
        frac_A_rep,
        bins=40,
        density=True,
        color="#2196F3",
        alpha=0.6,
        edgecolor="white",
        label="Replicated",
    )
    ax.axvline(
        frac_A_obs,
        color="#E53935",
        linewidth=2.5,
        label=f"Observed = {frac_A_obs:.3f}",
        zorder=10,
    )
    ax.set_xlabel("Fraction perfect (y=1)")
    ax.set_ylabel("Density")
    ax.set_title("PPC: P(perfect) Group A", fontsize=11, fontweight="bold")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)

    # P(perfect) Group B
    ax = axes[1]
    frac_B_rep = y_B_rep.mean(axis=1)
    frac_B_obs = self.y_B_obs.mean()
    ax.hist(
        frac_B_rep,
        bins=40,
        density=True,
        color="#4CAF50",
        alpha=0.6,
        edgecolor="white",
        label="Replicated",
    )
    ax.axvline(
        frac_B_obs,
        color="#E53935",
        linewidth=2.5,
        label=f"Observed = {frac_B_obs:.3f}",
        zorder=10,
    )
    ax.set_xlabel("Fraction perfect (y=1)")
    ax.set_ylabel("Density")
    ax.set_title("PPC: P(perfect) Group B", fontsize=11, fontweight="bold")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)

    # Rate difference
    ax = axes[2]
    diff_rep = frac_A_rep - frac_B_rep
    diff_obs = frac_A_obs - frac_B_obs
    ax.hist(
        diff_rep,
        bins=40,
        density=True,
        color="#9C27B0",
        alpha=0.6,
        edgecolor="white",
        label="Replicated",
    )
    ax.axvline(
        diff_obs,
        color="#E53935",
        linewidth=2.5,
        label=f"Observed = {diff_obs:.3f}",
        zorder=10,
    )
    ax.axvline(0, color="gray", linestyle="--", linewidth=1, alpha=0.6)
    ax.set_xlabel("P(perfect)_A \u2212 P(perfect)_B")
    ax.set_ylabel("Density")
    ax.set_title(
        "PPC: Rate Difference (A \u2212 B)", fontsize=11, fontweight="bold"
    )
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)

    fig.suptitle(
        kwargs.pop(
            "suptitle",
            kwargs.pop("title", "Posterior Predictive Checks (Laplace Binomial)"),
        ),
        fontsize=14,
        fontweight="bold",
        y=1.02,
    )
    plt.tight_layout()
    plt.show()

plot_sensitivity(prior_H0=0.5, **kwargs)

Two-panel sensitivity: P(H0|data) vs prior P(H0), and slab-width sweep.

Source code in bayesprop/resources/bayes_paired_laplace.py

def plot_sensitivity(self, prior_H0: float = 0.5, **kwargs) -> None:
    """Two-panel sensitivity: P(H0|data) vs prior P(H0), and slab-width sweep."""
    import matplotlib.pyplot as plt

    bf = self.savage_dickey_test()

    figsize = kwargs.pop("figsize", (14, 5))
    fig, axes = plt.subplots(1, 2, figsize=figsize)

    # Left: P(H0|data) vs prior P(H0)
    ax = axes[0]
    prior_grid = np.linspace(0.01, 0.99, 200)
    p_h0_grid = [
        self.posterior_probability_H0(bf.BF_01, p).p_H0 for p in prior_grid
    ]
    bf10 = bf.BF_10
    bf_label = (
        f"log\u2081\u2080BF\u2081\u2080={np.log10(bf10):.0f}"
        if bf10 > 1e4
        else f"BF\u2081\u2080={bf10:.1f}"
    )
    ax.plot(prior_grid, p_h0_grid, linewidth=2, label=bf_label)
    ax.axhline(
        0.05, color="red", linestyle="--", alpha=0.5, label="P(H\u2080)=0.05"
    )
    ax.axvline(prior_H0, color="gray", linestyle=":", alpha=0.5)
    ax.set_xlabel("Prior P(H\u2080)")
    ax.set_ylabel("Posterior P(H\u2080 | data)")
    ax.set_title(
        "Sensitivity: P(H\u2080 | data) vs Prior P(H\u2080)", fontweight="bold"
    )
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)

    # Right: BF_10 vs slab width sigma_s
    ax2 = axes[1]
    sigma_grid = np.linspace(0.25, 5.0, 100)
    # Analytical posterior density at 0 (Laplace marginal on δ_A is
    # exactly Gaussian), matching savage_dickey_test's closed form.
    sigma_post = float(np.sqrt(self.laplace["cov"][1, 1]))
    delta_A_map = float(self.laplace["map"][1])
    post_at_0 = float(
        np.exp(-0.5 * (delta_A_map / sigma_post) ** 2)
        / (sigma_post * np.sqrt(2.0 * np.pi))
    )
    bf10_vals = [norm.pdf(0, 0, s) / post_at_0 for s in sigma_grid]
    ax2.plot(sigma_grid, bf10_vals, linewidth=2)
    ax2.axhline(
        3, color="red", linestyle="--", alpha=0.5, label="BF\u2081\u2080 = 3"
    )
    ax2.axhline(
        1, color="gray", linestyle=":", alpha=0.5, label="BF\u2081\u2080 = 1"
    )
    ax2.axvline(
        self.prior_sigma_delta,
        color="gray",
        linestyle="--",
        alpha=0.3,
        label=f"\u03c3_s = {self.prior_sigma_delta} (used)",
    )
    ax2.set_xlabel("Slab width \u03c3_s")
    ax2.set_ylabel("BF\u2081\u2080")
    ax2.set_title("Sensitivity: BF\u2081\u2080 vs Slab Width", fontweight="bold")
    ax2.set_yscale("log")
    ax2.legend(fontsize=9)
    ax2.grid(alpha=0.3)

    fig.suptitle(
        kwargs.pop(
            "suptitle",
            kwargs.pop("title", "Jeffreys-Lindley Sensitivity (Binomial)"),
        ),
        fontsize=13,
        fontweight="bold",
        y=1.04,
    )
    plt.tight_layout()
    plt.show()

print_summary()

Print posterior summary, Savage-Dickey test, and PPC p-values.

Source code in bayesprop/resources/bayes_paired_laplace.py

def print_summary(self) -> None:
    """Print posterior summary, Savage-Dickey test, and PPC p-values."""
    self._check_fitted()

    mu_map, delta_map = self.laplace["map"]
    cov = self.laplace["cov"]
    n_A = self.laplace["n_A"]
    k_A = self.laplace["k_A"]
    n_B = self.laplace["n_B"]
    k_B = self.laplace["k_B"]

    # Laplace posterior info
    s = self.summary
    verdict = (
        "A wins"
        if s.p_A_greater_B > 0.95
        else ("Tied" if s.p_A_greater_B > 0.5 else "B wins")
    )
    print("Laplace posterior summary (Paired)")
    print("=" * 60)
    print(f"  \u03b8_A  mean={s.theta_A_mean:.4f}  " f"[n_A={n_A}, k_A={k_A}]")
    print(f"  \u03b8_B  mean={s.theta_B_mean:.4f}  " f"[n_B={n_B}, k_B={k_B}]")
    print(f"  Mean \u0394 (\u03b8_A \u2212 \u03b8_B):  {s.mean_delta:.4f}")
    print(f"  95% CI:              [{s.ci_95.lower:.4f}, {s.ci_95.upper:.4f}]")
    print(f"  P(A > B):            {s.p_A_greater_B:.4f}")
    print(f"  Verdict:             {verdict}")
    print()
    print("Logit-scale parameters")
    print("-" * 60)
    print(f"  MAP: \u03bc={mu_map:.4f}, \u03b4_A={delta_map:.4f}")
    print(
        f"  Posterior sd: \u03bc={np.sqrt(cov[0, 0]):.4f}, \u03b4_A={np.sqrt(cov[1, 1]):.4f}"
    )
    print(f"  Correlation: {cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1]):.3f}")
    if self.laplace["hierarchical"]:
        print()
        print("Learned prior scales (MAP)")
        print("-" * 60)
        print(f"  \u03c3_\u03bc = {self.laplace['sigma_mu_map']:.4f}")
        print(f"  \u03c3_\u03b4 = {self.laplace['sigma_delta_map']:.4f}")

    # Savage-Dickey
    bf = self.savage_dickey_test()
    print()
    print("Savage-Dickey Bayes Factor: H0 (\u0394 = 0) vs H1 (\u0394 \u2260 0)")
    print("=" * 60)
    print(f"  Prior  density at \u03b4=0: {bf.prior_density_at_0:.6f}")
    print(f"  Post.  density at \u03b4=0: {bf.posterior_density_at_0:.2e}")
    print(f"  BF_01 (for H0):        {_format_bf(bf.BF_01)}")
    print(f"  BF_10 (against H0):    {_format_bf(bf.BF_10)}")
    print(f"  log\u2081\u2080(BF_10):          {np.log10(bf.BF_10):.1f}")
    print(f"  \u2192 {bf.interpretation}")
    print(f"  \u2192 Decision: {bf.decision}")

    # Posterior probability of H0
    post = self.posterior_probability_H0(bf.BF_01)
    print()
    print("Posterior model probabilities (prior P(H0) = 0.5)")
    print("=" * 60)
    print(f"  P(H0|data): {post.p_H0:.2e}")
    print(f"  P(H1|data): {post.p_H1:.6f}")

    # PPC p-values
    ppc = self.ppc_pvalues()
    print()
    print("Posterior Predictive p-values")
    print("=" * 60)
    print(f"  {'Statistic':<20} {'Observed':>10} {'p-value':>10} {'Status':>8}")
    print("  " + "-" * 50)
    for stat, vals in ppc.items():
        print(
            f"  {stat:<20} {vals.observed:>10.4f} {vals.p_value:>10.3f} {vals.status:>8}"
        )

    # Trace summary
    print()
    print("Laplace trace diagnostics")
    print("=" * 60)
    print(self.trace_summary.to_string())

plot_forest(results, label_A='Group A', label_B='Group B', **kwargs) staticmethod

Forest plot + P(A>B) bar chart for multiple metrics.

Source code in bayesprop/resources/bayes_paired_laplace.py

@staticmethod
def plot_forest(
    results: dict[str, "_PairedLaplace"],
    label_A: str = "Group A",
    label_B: str = "Group B",
    **kwargs,
) -> None:
    """Forest plot + P(A>B) bar chart for multiple metrics."""
    import matplotlib.patches as mpatches
    import matplotlib.pyplot as plt

    metrics = list(results.keys())
    means = [results[m].summary.mean_delta for m in metrics]
    ci_lows = [results[m].summary.ci_95.lower for m in metrics]
    ci_highs = [results[m].summary.ci_95.upper for m in metrics]
    probs = [results[m].summary.p_A_greater_B for m in metrics]

    colors = [
        "#2196F3" if p > 0.95 else "#FF9800" if p > 0.5 else "#F44336"
        for p in probs
    ]
    y_pos = np.arange(len(metrics))

    figsize = kwargs.pop("figsize", (14, max(4, 2 * len(metrics))))
    fig, axes = plt.subplots(1, 2, figsize=figsize)

    ax = axes[0]
    for i, (m, ci_l, ci_h, col) in enumerate(
        zip(means, ci_lows, ci_highs, colors, strict=False)
    ):
        ax.plot(
            [ci_l, ci_h], [i, i], color=col, linewidth=2.5, solid_capstyle="round"
        )
        ax.plot(m, i, "o", color=col, markersize=8, zorder=5)
    ax.axvline(0, color="gray", linestyle="--", linewidth=1, alpha=0.7)
    ax.set_yticks(y_pos)
    ax.set_yticklabels(metrics, fontsize=11)
    ax.set_xlabel(
        f"Mean \u0394 P(perfect)\n\u2190 {label_B} better | {label_A} better \u2192"
    )
    ax.set_title("Posterior Mean Difference with 95% CI", fontweight="bold")
    ax.invert_yaxis()
    ax.grid(axis="x", alpha=0.3)

    ax2 = axes[1]
    ax2.barh(y_pos, probs, color=colors, height=0.5, alpha=0.85)
    ax2.axvline(0.5, color="gray", linestyle="--", linewidth=1, alpha=0.7)
    ax2.set_yticks(y_pos)
    ax2.set_yticklabels(metrics, fontsize=11)
    ax2.set_xlabel(f"P({label_A} > {label_B})")
    ax2.set_title("Posterior Probability of Superiority", fontweight="bold")
    ax2.set_xlim(0, 1.05)
    ax2.invert_yaxis()
    ax2.grid(axis="x", alpha=0.3)
    for i, p in enumerate(probs):
        ax2.text(
            p + 0.02, i, f"{p:.2f}", va="center", fontsize=10, fontweight="bold"
        )

    legend_elements = [
        mpatches.Patch(color="#2196F3", label="Strong (P > 0.95)"),
        mpatches.Patch(color="#FF9800", label="Moderate (0.5 < P \u2264 0.95)"),
        mpatches.Patch(color="#F44336", label="Reversed (P \u2264 0.5)"),
    ]
    fig.legend(
        handles=legend_elements,
        loc="lower center",
        ncol=3,
        fontsize=9,
        bbox_to_anchor=(0.5, -0.02),
    )
    fig.suptitle(
        kwargs.pop(
            "suptitle",
            kwargs.pop(
                "title", f"{label_A} vs {label_B} \u2014 Pooled Binomial Comparison"
            ),
        ),
        fontsize=14,
        fontweight="bold",
        y=1.02,
    )
    plt.tight_layout()
    plt.show()

print_comparison_table(results) staticmethod

Print a formatted comparison table across metrics.

Source code in bayesprop/resources/bayes_paired_laplace.py

@staticmethod
def print_comparison_table(results: dict[str, "_PairedLaplace"]) -> None:
    """Print a formatted comparison table across metrics."""
    print("=" * 80)
    print(
        f"{'Metric':<25} {'Mean \u0394':>8} {'95% CI':>20} {'P(A>B)':>8} {'Verdict':>12}"
    )
    print("=" * 80)
    for m, model in results.items():
        s = model.summary
        verdict = (
            "A wins"
            if s.p_A_greater_B > 0.95
            else ("Tied" if s.p_A_greater_B > 0.5 else "B wins")
        )
        print(
            f"{m:<25} {s.mean_delta:>8.4f} "
            f"[{s.ci_95.lower:>7.4f}, {s.ci_95.upper:>7.4f}] "
            f"{s.p_A_greater_B:>8.4f} {verdict:>12}"
        )
    print("=" * 80)

sigmoid(x)

Element-wise sigmoid (logistic) function.

Source code in bayesprop/resources/bayes_paired_laplace.py

def sigmoid(x: npt.ArrayLike) -> np.ndarray:
    """Element-wise sigmoid (logistic) function."""
    return 1.0 / (1.0 + np.exp(-x))

_format_bf(value)

Format a Bayes Factor for human-readable display.

Source code in bayesprop/resources/bayes_paired_laplace.py

def _format_bf(value: float) -> str:
    """Format a Bayes Factor for human-readable display."""
    if value > 1e4:
        return f"10^{np.log10(value):.0f}"
    elif value < 1e-4 and value > 0:
        return f"10^{np.log10(value):.0f}"
    else:
        return f"{value:.2f}"

_paired_laplace_from_counts(n_A, k_A, n_B, k_B, prior_sigma_delta, prior_sigma_mu=2.0, x0=(0.0, 0.0), tol=1e-08, max_iter=50)

Compute the Laplace posterior of (mu, delta_A) directly from counts.

Solves for the MAP via damped Newton iterations using the closed-form gradient and Hessian of the pooled Bernoulli logistic log-posterior (no raw data is materialised, no external optimizer invoked). The objective depends on the data only through the four sufficient statistics (n_A, k_A, n_B, k_B).

Newton converges quadratically; warm-starting x0 from the previous MAP (as the sequential test does) typically requires only 1-3 iterations per update.

Parameters:

Name	Type	Description	Default
n_A	int	Cumulative sample size for arm A.	required
k_A	int	Cumulative successes for arm A.	required
n_B	int	Cumulative sample size for arm B.	required
k_B	int	Cumulative successes for arm B.	required
prior_sigma_delta	float	Std of the N(0, sigma) prior on delta_A.	required
prior_sigma_mu	float	Std of the N(0, sigma) prior on mu.	2.0
x0	tuple[float, float]	Warm-start for Newton as (mu0, delta0).	(0.0, 0.0)
tol	float	Convergence tolerance on the gradient infinity-norm.	1e-08
max_iter	int	Maximum number of Newton iterations.	50

Returns:

Type	Description
ndarray	Tuple (theta_map, cov, H) where theta_map is the MAP
ndarray	of (mu, delta_A), cov is the 2x2 Laplace covariance
ndarray	(closed-form inverse of the observed information), and H is
tuple[ndarray, ndarray, ndarray]	the 2x2 Hessian of the negative log-posterior at the MAP.

Source code in bayesprop/resources/bayes_paired_laplace.py

def _paired_laplace_from_counts(
    n_A: int,
    k_A: int,
    n_B: int,
    k_B: int,
    prior_sigma_delta: float,
    prior_sigma_mu: float = 2.0,
    x0: tuple[float, float] = (0.0, 0.0),
    tol: float = 1e-8,
    max_iter: int = 50,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Compute the Laplace posterior of (mu, delta_A) directly from counts.

    Solves for the MAP via damped Newton iterations using the closed-form
    gradient and Hessian of the pooled Bernoulli logistic log-posterior
    (no raw data is materialised, no external optimizer invoked). The
    objective depends on the data only through the four sufficient
    statistics ``(n_A, k_A, n_B, k_B)``.

    Newton converges quadratically; warm-starting ``x0`` from the
    previous MAP (as the sequential test does) typically requires only
    1-3 iterations per update.

    Args:
        n_A: Cumulative sample size for arm A.
        k_A: Cumulative successes for arm A.
        n_B: Cumulative sample size for arm B.
        k_B: Cumulative successes for arm B.
        prior_sigma_delta: Std of the N(0, sigma) prior on delta_A.
        prior_sigma_mu: Std of the N(0, sigma) prior on mu.
        x0: Warm-start for Newton as (mu0, delta0).
        tol: Convergence tolerance on the gradient infinity-norm.
        max_iter: Maximum number of Newton iterations.

    Returns:
        Tuple ``(theta_map, cov, H)`` where ``theta_map`` is the MAP
        of ``(mu, delta_A)``, ``cov`` is the 2x2 Laplace covariance
        (closed-form inverse of the observed information), and ``H`` is
        the 2x2 Hessian of the negative log-posterior at the MAP.
    """
    inv_var_mu = 1.0 / (prior_sigma_mu**2)
    inv_var_delta = 1.0 / (prior_sigma_delta**2)

    def neg_log_post(mu_: float, delta_: float) -> float:
        """Closed-form negative log-posterior up to an additive constant."""
        zA = mu_ + delta_
        zB = mu_
        # softplus via np.logaddexp(0, x) for numerical stability.
        nll = (
            k_A * np.logaddexp(0.0, -zA)
            + (n_A - k_A) * np.logaddexp(0.0, zA)
            + k_B * np.logaddexp(0.0, -zB)
            + (n_B - k_B) * np.logaddexp(0.0, zB)
        )
        nlp = 0.5 * inv_var_mu * mu_ * mu_ + 0.5 * inv_var_delta * delta_ * delta_
        return float(nll + nlp)

    def hessian_terms(p_A: float, p_B: float) -> tuple[float, float, float, float]:
        """Closed-form Hessian of the negative log-posterior at (mu, delta).

        Returns the unique entries (a, b, c) of the symmetric 2x2 matrix
        ``H = [[a, c], [c, b]]`` (always positive definite here) along
        with its determinant ``det = a*b - c*c``.
        """
        w_A = p_A * (1.0 - p_A)
        w_B = p_B * (1.0 - p_B)
        a_ = n_A * w_A + n_B * w_B + inv_var_mu
        b_ = n_A * w_A + inv_var_delta
        c_ = n_A * w_A
        return a_, b_, c_, a_ * b_ - c_ * c_

    mu, delta = float(x0[0]), float(x0[1])
    a = b = c = det = 0.0  # populated each iteration; reused for cov
    for _ in range(max_iter):
        # Gradient of negative log-posterior.
        p_A = 1.0 / (1.0 + np.exp(-(mu + delta)))
        p_B = 1.0 / (1.0 + np.exp(-mu))
        g_mu = -((k_A - n_A * p_A) + (k_B - n_B * p_B) - inv_var_mu * mu)
        g_delta = -((k_A - n_A * p_A) - inv_var_delta * delta)

        a, b, c, det = hessian_terms(p_A, p_B)

        if max(abs(g_mu), abs(g_delta)) < tol:
            break

        # Newton step: dx = -H^{-1} g, using closed-form 2x2 inverse.
        d_mu = -(b * g_mu - c * g_delta) / det
        d_delta = -(-c * g_mu + a * g_delta) / det

        # Backtracking line search (Armijo with c1=1e-4) to guarantee descent
        # even from a poor warm-start. The full Newton step is accepted on the
        # first try in the typical warm-started case, so this adds ~1 objective
        # evaluation per iteration in steady state.
        old_obj = neg_log_post(mu, delta)
        directional = g_mu * d_mu + g_delta * d_delta  # = -g^T H^{-1} g <= 0
        step = 1.0
        for _ in range(20):
            new_mu = mu + step * d_mu
            new_delta = delta + step * d_delta
            if neg_log_post(new_mu, new_delta) <= old_obj + 1e-4 * step * directional:
                break
            step *= 0.5
        mu, delta = new_mu, new_delta
    else:
        grad_norm = max(abs(g_mu), abs(g_delta))
        warnings.warn(
            f"_paired_laplace_from_counts did not converge in {max_iter} "
            f"iterations (||∇||_∞ = {grad_norm:.2e}, tol = {tol:.1e}). "
            "Consider increasing max_iter or relaxing tol.",
            stacklevel=2,
        )

    H = np.array([[a, c], [c, b]])
    cov = np.array([[b, -c], [-c, a]]) / det

    return np.array([mu, delta]), cov, H

_hierarchical_laplace_from_counts(n_A, k_A, n_B, k_B, hp_mu, hp_delta, x0=None, tol=1e-08, max_iter=100)

Compute the hierarchical Laplace posterior of (μ, δ_A) with learned scales.

Places Inverse-Gamma hyperpriors on the prior variances::

σ²_μ  ~ IG(a_μ, b_μ)
σ²_δ  ~ IG(a_δ, b_δ)
μ     ~ N(0, σ²_μ)
δ_A   ~ N(0, σ²_δ)

Optimisation is performed in the reparameterised space (μ, δ_A, ψ_μ, ψ_δ) where ψ = log σ to enforce positivity. Newton iterations use the closed-form 4×4 gradient and Hessian with Armijo backtracking.

The marginal posterior on (μ, δ_A) is extracted from the top-left 2×2 block of the 4×4 Laplace covariance, which already incorporates uncertainty from the learned scales.

Parameters:

Name	Type	Description	Default
n_A	int	Cumulative sample size for arm A.	required
k_A	int	Cumulative successes for arm A.	required
n_B	int	Cumulative sample size for arm B.	required
k_B	int	Cumulative successes for arm B.	required
hp_mu	tuple[float, float]	(a, b) shape and scale of the IG hyperprior on σ²_μ.	required
hp_delta	tuple[float, float]	(a, b) shape and scale of the IG hyperprior on σ²_δ.	required
x0	tuple[float, float, float, float] | None	Warm-start as (μ, δ, log σ_μ, log σ_δ). Defaults to (0, 0, log 2, 0).	None
tol	float	Convergence tolerance on the gradient infinity-norm.	1e-08
max_iter	int	Maximum number of Newton iterations.	100

Returns:

Type	Description
ndarray	Tuple (theta_map, cov, H, sigma_mu_map, sigma_delta_map)
ndarray	where theta_map is the 2-D MAP [μ, δ_A], cov is
ndarray	the marginal 2×2 Laplace covariance, H is the
float	corresponding 2×2 marginal precision, and
float	sigma_mu_map / sigma_delta_map are the MAP estimates
tuple[ndarray, ndarray, ndarray, float, float]	of the learned prior standard deviations.

Source code in bayesprop/resources/bayes_paired_laplace.py

def _hierarchical_laplace_from_counts(
    n_A: int,
    k_A: int,
    n_B: int,
    k_B: int,
    hp_mu: tuple[float, float],
    hp_delta: tuple[float, float],
    x0: tuple[float, float, float, float] | None = None,
    tol: float = 1e-8,
    max_iter: int = 100,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float]:
    """Compute the hierarchical Laplace posterior of (μ, δ_A) with learned scales.

    Places Inverse-Gamma hyperpriors on the prior variances::

        σ²_μ  ~ IG(a_μ, b_μ)
        σ²_δ  ~ IG(a_δ, b_δ)
        μ     ~ N(0, σ²_μ)
        δ_A   ~ N(0, σ²_δ)

    Optimisation is performed in the reparameterised space
    ``(μ, δ_A, ψ_μ, ψ_δ)`` where ``ψ = log σ`` to enforce
    positivity.  Newton iterations use the closed-form 4×4 gradient
    and Hessian with Armijo backtracking.

    The marginal posterior on ``(μ, δ_A)`` is extracted from the
    top-left 2×2 block of the 4×4 Laplace covariance, which already
    incorporates uncertainty from the learned scales.

    Args:
        n_A: Cumulative sample size for arm A.
        k_A: Cumulative successes for arm A.
        n_B: Cumulative sample size for arm B.
        k_B: Cumulative successes for arm B.
        hp_mu: ``(a, b)`` shape and scale of the IG hyperprior on σ²_μ.
        hp_delta: ``(a, b)`` shape and scale of the IG hyperprior on σ²_δ.
        x0: Warm-start as ``(μ, δ, log σ_μ, log σ_δ)``.  Defaults to
            ``(0, 0, log 2, 0)``.
        tol: Convergence tolerance on the gradient infinity-norm.
        max_iter: Maximum number of Newton iterations.

    Returns:
        Tuple ``(theta_map, cov, H, sigma_mu_map, sigma_delta_map)``
        where ``theta_map`` is the 2-D MAP ``[μ, δ_A]``, ``cov`` is
        the marginal 2×2 Laplace covariance, ``H`` is the
        corresponding 2×2 marginal precision, and
        ``sigma_mu_map`` / ``sigma_delta_map`` are the MAP estimates
        of the learned prior standard deviations.
    """
    a_mu, b_mu = hp_mu
    a_delta, b_delta = hp_delta

    if x0 is None:
        x0 = (0.0, 0.0, np.log(2.0), 0.0)

    def neg_log_post(
        mu_: float, delta_: float, psi_mu_: float, psi_delta_: float
    ) -> float:
        """Closed-form negative log-posterior (up to additive constant).

        Evaluated in the reparameterised space ``(μ, δ, ψ_μ, ψ_δ)``
        where ``ψ = log σ``.  The Jacobian of the log-transform is
        absorbed into the prior term.

        Args:
            mu_: Intercept on the logit scale.
            delta_: Treatment effect on the logit scale.
            psi_mu_: Log prior standard deviation for μ.
            psi_delta_: Log prior standard deviation for δ_A.

        Returns:
            Scalar value of the negative log-posterior.
        """
        zA = mu_ + delta_
        zB = mu_
        tau_mu_ = np.exp(-2.0 * psi_mu_)
        tau_delta_ = np.exp(-2.0 * psi_delta_)
        nll = (
            k_A * np.logaddexp(0.0, -zA)  # np.log(1 + np.exp(z))
            + (n_A - k_A) * np.logaddexp(0.0, zA)
            + k_B * np.logaddexp(0.0, -zB)
            + (n_B - k_B) * np.logaddexp(0.0, zB)
        )
        nlp = (
            (2.0 * a_mu + 1.0) * psi_mu_
            + (mu_**2 / 2.0 + b_mu) * tau_mu_
            + (2.0 * a_delta + 1.0) * psi_delta_
            + (delta_**2 / 2.0 + b_delta) * tau_delta_
        )
        return float(nll + nlp)

    mu, delta = float(x0[0]), float(x0[1])
    psi_mu, psi_delta = float(x0[2]), float(x0[3])

    for _ in range(max_iter):
        p_A = 1.0 / (1.0 + np.exp(-(mu + delta)))
        p_B = 1.0 / (1.0 + np.exp(-mu))
        w_A = p_A * (1.0 - p_A)
        w_B = p_B * (1.0 - p_B)
        tau_mu = np.exp(-2.0 * psi_mu)
        tau_delta = np.exp(-2.0 * psi_delta)

        r_A = k_A - n_A * p_A
        r_B = k_B - n_B * p_B

        # 4-D gradient of the negative log-posterior.
        g = np.array(
            [
                -(r_A + r_B) + mu * tau_mu,
                -r_A + delta * tau_delta,
                (2.0 * a_mu + 1.0) - (mu**2 + 2.0 * b_mu) * tau_mu,
                (2.0 * a_delta + 1.0) - (delta**2 + 2.0 * b_delta) * tau_delta,
            ]
        )

        if np.max(np.abs(g)) < tol:
            break

        # 4×4 Hessian of the negative log-posterior.
        H4 = np.zeros((4, 4))
        H4[0, 0] = n_A * w_A + n_B * w_B + tau_mu
        H4[0, 1] = H4[1, 0] = n_A * w_A
        H4[1, 1] = n_A * w_A + tau_delta
        H4[0, 2] = H4[2, 0] = -2.0 * mu * tau_mu
        H4[1, 3] = H4[3, 1] = -2.0 * delta * tau_delta
        H4[2, 2] = 2.0 * (mu**2 + 2.0 * b_mu) * tau_mu
        H4[3, 3] = 2.0 * (delta**2 + 2.0 * b_delta) * tau_delta

        # Newton step via dense solve.
        try:
            step_dir = np.linalg.solve(H4, -g)
        except np.linalg.LinAlgError:
            H4 += 1e-6 * np.eye(4)
            step_dir = np.linalg.solve(H4, -g)

        # Armijo backtracking line search.
        old_obj = neg_log_post(mu, delta, psi_mu, psi_delta)
        directional = float(g @ step_dir)
        alpha = 1.0
        for _ in range(20):
            new_mu = mu + alpha * step_dir[0]
            new_delta = delta + alpha * step_dir[1]
            new_psi_mu = psi_mu + alpha * step_dir[2]
            new_psi_delta = psi_delta + alpha * step_dir[3]
            if (
                neg_log_post(new_mu, new_delta, new_psi_mu, new_psi_delta)
                <= old_obj + 1e-4 * alpha * directional
            ):
                break
            alpha *= 0.5
        mu, delta = new_mu, new_delta
        psi_mu, psi_delta = new_psi_mu, new_psi_delta
    else:
        grad_norm = float(np.max(np.abs(g)))
        warnings.warn(
            f"_hierarchical_laplace_from_counts did not converge in "
            f"{max_iter} iterations (||∇||_∞ = {grad_norm:.2e}, "
            f"tol = {tol:.1e}). Consider increasing max_iter or "
            "relaxing tol.",
            stacklevel=2,
        )

    # Recompute Hessian at converged point for the covariance.
    p_A = 1.0 / (1.0 + np.exp(-(mu + delta)))
    p_B = 1.0 / (1.0 + np.exp(-mu))
    w_A = p_A * (1.0 - p_A)
    w_B = p_B * (1.0 - p_B)
    tau_mu = np.exp(-2.0 * psi_mu)
    tau_delta = np.exp(-2.0 * psi_delta)

    H4 = np.zeros((4, 4))
    H4[0, 0] = n_A * w_A + n_B * w_B + tau_mu
    H4[0, 1] = H4[1, 0] = n_A * w_A
    H4[1, 1] = n_A * w_A + tau_delta
    H4[0, 2] = H4[2, 0] = -2.0 * mu * tau_mu
    H4[1, 3] = H4[3, 1] = -2.0 * delta * tau_delta
    H4[2, 2] = 2.0 * (mu**2 + 2.0 * b_mu) * tau_mu
    H4[3, 3] = 2.0 * (delta**2 + 2.0 * b_delta) * tau_delta

    cov_full = np.linalg.inv(H4)
    cov_2d = cov_full[:2, :2]
    H_2d = np.linalg.inv(cov_2d)

    sigma_mu_map = float(np.exp(psi_mu))
    sigma_delta_map = float(np.exp(psi_delta))

    return np.array([mu, delta]), cov_2d, H_2d, sigma_mu_map, sigma_delta_map