Operating Characteristics (Paired)

Monte-Carlo evaluation of the paired-Laplace Bayesian procedure's frequentist operating characteristics — three-way decision rates, CI coverage on Δ = p_A − p_B, matched-α McNemar baselines, and sequential stopping-time distributions. The companion notebook is src/notebooks/operating_characteristics_paired_laplace.ipynb; the narrative is in Frequentist Evaluation — Paired Laplace.

operation_characteristics_paired

Operating-characteristic / BFDA simulation harness for the paired-Laplace model.

This module mirrors :mod:bayesprop.utils.operation_characteristics but evaluates the paired Bayesian procedure built around :class:bayesprop.resources.bayes_paired_laplace.PairedBayesPropTest and its sequential variant :class:SequentialPairedBayesPropTest.

A paired binary design has the same three frequentist concerns as the non-paired one — Type-I rate, three-way decision rates, credible-interval coverage — plus a sequential variant. The analysis model is y_A ~ Bern(σ(μ + δ_A)), y_B ~ Bern(σ(μ)), with prior δ_A ~ N(0, σ_δ²). The Savage–Dickey BF tests H_0: δ_A = 0, which is equivalent to H_0: p_A = p_B.

For an apples-to-apples comparison with the non-paired analysis (and with the user's intuition), the public API takes the true marginal probabilities (p_A, p_B) on the probability scale and inverts the sigmoid internally to recover the analysis-model intercept μ = logit(p_B) and effect δ_A = logit(p_A) − logit(p_B). This parametrisation is the one expected by :func:simulate_paired_scores when δ_B = 0 and σ_θ = 0.

Three public entry points cover the workflow:

• :func:simulate_fixed_n_paired — one (p_A, p_B, n) cell of the OC grid, returning the Bayes three-way decision rates, the 95 % CI coverage of Δ = p_A − p_B, and the per-replicate McNemar exact p-values needed for a frequentist baseline.
• :func:grid_fixed_n_paired — sweeps :func:simulate_fixed_n_paired over an arbitrary list of (p_A, p_B) pairs.
• :func:simulate_sequential_paired — empirical stopping-time distribution of :class:SequentialPairedBayesPropTest at a single (p_A, p_B) point.

The matched-α helper :func:matched_calibration_alpha and the Wilson band helper :func:wilson_band are parametrisation-free and are re-exported from the non-paired module so callers don't have to import from two places.

matched_calibration_alpha(p_values, bayes_type1_rate, null_grid_index)

Empirical α that matches a frequentist test's Type-I rate to the Bayes rule's.

Given a matrix of Fisher (or other) p-values produced under a grid of true effects, the row indexed by null_grid_index corresponds to the null case Δ = 0. We pick the p-value cutoff α whose empirical rejection rate on that row equals bayes_type1_rate, by reading the corresponding quantile of the null p-value distribution. The resulting α is the fairest like-for-like calibration when overlaying the frequentist and Bayesian power curves.

Parameters:

Name	Type	Description	Default
p_values	ndarray	(n_grid, n_sim) matrix of p-values from a grid simulation. Typically the second return of :func:grid_fixed_n.	required
bayes_type1_rate	float	Empirical P(reject H₀ | Δ = 0) of the Bayes rule on that same grid (i.e. rates_df.iloc[null_grid_index]["reject"]).	required
null_grid_index	int	Row index into p_values corresponding to the null grid cell Δ = 0.	required

Returns:

Type	Description
float	Matched-calibration α, clipped to [0, 1]. Returns 0.0
float	if bayes_type1_rate <= 0 (the Bayes rule never rejects at
float	the null in this simulation), since no positive α can match.

Raises:

Type	Description
ValueError	If null_grid_index is out of bounds or p_values is not 2-D.

Source code in bayesprop/utils/operation_characteristics.py

def matched_calibration_alpha(
    p_values: np.ndarray,
    bayes_type1_rate: float,
    null_grid_index: int,
) -> float:
    """Empirical α that matches a frequentist test's Type-I rate to the Bayes rule's.

    Given a matrix of Fisher (or other) p-values produced under a grid
    of true effects, the row indexed by ``null_grid_index`` corresponds
    to the null case ``Δ = 0``. We pick the p-value cutoff ``α`` whose
    empirical rejection rate on that row equals
    ``bayes_type1_rate``, by reading the corresponding quantile of the
    null p-value distribution. The resulting ``α`` is the fairest
    like-for-like calibration when overlaying the frequentist and
    Bayesian power curves.

    Args:
        p_values: ``(n_grid, n_sim)`` matrix of p-values from a grid
            simulation. Typically the second return of
            :func:`grid_fixed_n`.
        bayes_type1_rate: Empirical ``P(reject H₀ | Δ = 0)`` of the
            Bayes rule on that same grid (i.e.
            ``rates_df.iloc[null_grid_index]["reject"]``).
        null_grid_index: Row index into ``p_values`` corresponding to
            the null grid cell ``Δ = 0``.

    Returns:
        Matched-calibration α, clipped to ``[0, 1]``. Returns ``0.0``
        if ``bayes_type1_rate <= 0`` (the Bayes rule never rejects at
        the null in this simulation), since no positive α can match.

    Raises:
        ValueError: If ``null_grid_index`` is out of bounds or
            ``p_values`` is not 2-D.
    """
    if p_values.ndim != 2:
        raise ValueError(
            f"p_values must be 2-D (n_grid, n_sim); got shape {p_values.shape}"
        )
    if not 0 <= null_grid_index < p_values.shape[0]:
        raise ValueError(
            f"null_grid_index={null_grid_index} out of range for "
            f"p_values of shape {p_values.shape}"
        )
    if bayes_type1_rate <= 0.0:
        return 0.0

    # The empirical α-quantile of the null p-value distribution is by
    # definition the cutoff whose tail mass equals bayes_type1_rate.
    alpha = float(np.quantile(p_values[null_grid_index], bayes_type1_rate))
    return float(np.clip(alpha, 0.0, 1.0))

wilson_band(rates, n_sim, *, confidence=0.95)

Wilson confidence band for a vector of Monte-Carlo binomial rates.

Each entry of rates is treated as an empirical binomial proportion k / n_sim with k = round(rate * n_sim) (the number of successes among the n_sim simulation replicates). The Wilson score interval is computed pointwise via :func:scipy.stats.binomtest's proportion_ci(method="wilson").

Unlike the Wald interval, Wilson stays inside [0, 1] and retains its nominal coverage near the boundaries (rate ≈ 0 or ≈ 1), which is exactly where OC curves spend most of their interesting structure. See Brown, Cai & DasGupta (2001).

Parameters:

Name	Type	Description	Default
rates	ndarray	1-D array of empirical rates in [0, 1] (e.g. one column of an OC DataFrame).	required
n_sim	int	Number of simulation replicates each rate was averaged over. Constant across the grid in this codebase.	required
confidence	float	Nominal coverage, e.g. 0.95 for a 95 % band. Must be in (0, 1).	0.95

Returns:

Type	Description
ndarray	Tuple (lower, upper) of 1-D np.ndarray matching the shape
ndarray	of rates. Both endpoints are clipped to [0, 1] by
tuple[ndarray, ndarray]	construction of the Wilson interval.

Raises:

Type	Description
ValueError	If n_sim < 1, confidence is not in (0, 1), or any rate lies outside [0, 1].

Source code in bayesprop/utils/operation_characteristics.py

def wilson_band(
    rates: np.ndarray,
    n_sim: int,
    *,
    confidence: float = 0.95,
) -> tuple[np.ndarray, np.ndarray]:
    """Wilson confidence band for a vector of Monte-Carlo binomial rates.

    Each entry of ``rates`` is treated as an empirical binomial proportion
    ``k / n_sim`` with ``k = round(rate * n_sim)`` (the number of successes
    among the ``n_sim`` simulation replicates). The Wilson score interval
    is computed pointwise via
    :func:`scipy.stats.binomtest`'s ``proportion_ci(method="wilson")``.

    Unlike the Wald interval, Wilson stays inside ``[0, 1]`` and retains
    its nominal coverage near the boundaries (``rate ≈ 0`` or ``≈ 1``),
    which is exactly where OC curves spend most of their interesting
    structure. See Brown, Cai & DasGupta (2001).

    Args:
        rates: 1-D array of empirical rates in ``[0, 1]`` (e.g. one
            column of an OC ``DataFrame``).
        n_sim: Number of simulation replicates each rate was averaged
            over. Constant across the grid in this codebase.
        confidence: Nominal coverage, e.g. ``0.95`` for a 95 % band.
            Must be in ``(0, 1)``.

    Returns:
        Tuple ``(lower, upper)`` of 1-D ``np.ndarray`` matching the shape
        of ``rates``. Both endpoints are clipped to ``[0, 1]`` by
        construction of the Wilson interval.

    Raises:
        ValueError: If ``n_sim < 1``, ``confidence`` is not in ``(0, 1)``,
            or any ``rate`` lies outside ``[0, 1]``.
    """
    rates = np.asarray(rates, dtype=float)
    if n_sim < 1:
        raise ValueError(f"n_sim must be >= 1; got {n_sim}")
    if not 0.0 < confidence < 1.0:
        raise ValueError(f"confidence must be in (0, 1); got {confidence}")
    if np.any((rates < 0.0) | (rates > 1.0)):
        raise ValueError("All rates must lie in [0, 1].")

    # binomtest works on integer counts, so we recover k from each rate.
    # round() is the right inverse here: rates produced by the harness
    # are themselves k / n_sim for integer k.
    counts = np.rint(rates * n_sim).astype(int)
    lower = np.empty_like(rates)
    upper = np.empty_like(rates)
    for i, k in enumerate(counts):
        ci = binomtest(int(k), n_sim).proportion_ci(
            confidence_level=confidence, method="wilson"
        )
        lower[i] = ci.low
        upper[i] = ci.high
    return lower, upper

simulate_fixed_n_paired(p_A, p_B, n, n_sim, rng, *, prior_sigma_delta=1.0, bf_upper=3.0, bf_lower=1.0 / 3.0, n_samples_mc=4000, track_ci=True)

Monte-Carlo operating characteristics at a single (p_A, p_B, n) cell.

For each of n_sim replicates we

1. generate one paired dataset of size n via :func:simulate_paired_scores with theta_A = p_A and theta_B = p_B (σ_θ = 0);
2. fit :class:PairedBayesPropTest with prior δ_A ~ N(0, prior_sigma_delta²);
3. compute the Savage–Dickey BF on δ_A = 0 and classify the result with :func:classify_bf (configurable bf_upper / bf_lower);
4. optionally check whether the 95 % credible interval on Δ = p_A − p_B covers the true effect (track_ci=True);
5. run :func:mcnemar_paired_test on the same simulated data so a frequentist baseline can be derived later.

Parameters:

Name	Type	Description	Default
p_A	float	True success probability for arm A.	required
p_B	float	True success probability for arm B (= the analysis-model intercept on the probability scale).	required
n	int	Number of paired observations.	required
n_sim	int	Number of Monte-Carlo replicates.	required
rng	Generator	Pre-seeded NumPy generator. Threaded through the entire loop so the harness is fully deterministic given the seed.	required
prior_sigma_delta	float	Standard deviation of the N(0, σ_δ²) prior on δ_A (logit scale).	1.0
bf_upper	float	Threshold above which the Bayes rule rejects H_0.	3.0
bf_lower	float	Threshold below which the Bayes rule accepts H_0. Must satisfy 0 < bf_lower < bf_upper.	1.0 / 3.0
n_samples_mc	int	Posterior Laplace draws per replicate. Affects the 95 % CI estimate; the BF is computed analytically.	4000
track_ci	bool	If False, skip the CI coverage check (ci_coverage will be NaN). Useful for pure null sweeps where coverage is not the quantity of interest.	True

Returns:

Type	Description
dict[str, float]	Tuple (summary, freq_pvalues). summary is a dict with
ndarray	keys "reject" / "accept" / "inconclusive" (Bayes
tuple[dict[str, float], ndarray]	decision rates, summing to 1), "ci_coverage" (empirical
tuple[dict[str, float], ndarray]	95 % CI coverage of Δ, or NaN if track_ci=False),
tuple[dict[str, float], ndarray]	and the echoes "p_A", "p_B", "delta", "n",
tuple[dict[str, float], ndarray]	"n_sim". freq_pvalues is a 1-D np.ndarray of
tuple[dict[str, float], ndarray]	length n_sim holding the per-replicate McNemar p-values.

Raises:

Type	Description
ValueError	If bf_lower >= bf_upper or bf_lower <= 0 (validated inside :func:classify_bf).

Source code in bayesprop/utils/operation_characteristics_paired.py

def simulate_fixed_n_paired(
    p_A: float,
    p_B: float,
    n: int,
    n_sim: int,
    rng: np.random.Generator,
    *,
    prior_sigma_delta: float = 1.0,
    bf_upper: float = 3.0,
    bf_lower: float = 1.0 / 3.0,
    n_samples_mc: int = 4000,
    track_ci: bool = True,
) -> tuple[dict[str, float], np.ndarray]:
    """Monte-Carlo operating characteristics at a single ``(p_A, p_B, n)`` cell.

    For each of ``n_sim`` replicates we

    1. generate one paired dataset of size ``n`` via
       :func:`simulate_paired_scores` with ``theta_A = p_A`` and
       ``theta_B = p_B`` (``σ_θ = 0``);
    2. fit :class:`PairedBayesPropTest` with prior
       ``δ_A ~ N(0, prior_sigma_delta²)``;
    3. compute the Savage–Dickey BF on ``δ_A = 0`` and classify the
       result with :func:`classify_bf` (configurable ``bf_upper`` /
       ``bf_lower``);
    4. optionally check whether the 95 % credible interval on
       ``Δ = p_A − p_B`` covers the true effect (``track_ci=True``);
    5. run :func:`mcnemar_paired_test` on the *same* simulated data
       so a frequentist baseline can be derived later.

    Args:
        p_A: True success probability for arm A.
        p_B: True success probability for arm B (= the analysis-model
            intercept on the probability scale).
        n: Number of paired observations.
        n_sim: Number of Monte-Carlo replicates.
        rng: Pre-seeded NumPy generator. Threaded through the entire
            loop so the harness is fully deterministic given the seed.
        prior_sigma_delta: Standard deviation of the ``N(0, σ_δ²)``
            prior on ``δ_A`` (logit scale).
        bf_upper: Threshold above which the Bayes rule rejects ``H_0``.
        bf_lower: Threshold below which the Bayes rule accepts ``H_0``.
            Must satisfy ``0 < bf_lower < bf_upper``.
        n_samples_mc: Posterior Laplace draws per replicate. Affects
            the 95 % CI estimate; the BF is computed analytically.
        track_ci: If False, skip the CI coverage check
            (``ci_coverage`` will be ``NaN``). Useful for pure null
            sweeps where coverage is not the quantity of interest.

    Returns:
        Tuple ``(summary, freq_pvalues)``. ``summary`` is a dict with
        keys ``"reject"`` / ``"accept"`` / ``"inconclusive"`` (Bayes
        decision rates, summing to 1), ``"ci_coverage"`` (empirical
        95 % CI coverage of ``Δ``, or ``NaN`` if ``track_ci=False``),
        and the echoes ``"p_A"``, ``"p_B"``, ``"delta"``, ``"n"``,
        ``"n_sim"``. ``freq_pvalues`` is a 1-D ``np.ndarray`` of
        length ``n_sim`` holding the per-replicate McNemar p-values.

    Raises:
        ValueError: If ``bf_lower >= bf_upper`` or ``bf_lower <= 0``
            (validated inside :func:`classify_bf`).
    """
    counts: dict[str, int] = {"reject": 0, "accept": 0, "inconclusive": 0}
    ci_covered = 0
    delta_true = p_A - p_B
    pvals = np.empty(n_sim, dtype=float)

    for i in range(n_sim):
        sim = simulate_paired_scores(
            N=n,
            theta_A=p_A,
            theta_B=p_B,
            sigma_theta=0.0,
            rng=rng,
        )

        bb = PairedBayesPropTest(
            prior_sigma_delta=prior_sigma_delta,
            n_samples=n_samples_mc,
            seed=int(rng.integers(0, 2**31 - 1)),
        )
        bb.fit(sim.y_A, sim.y_B)
        bf10 = bb.savage_dickey_test().BF_10
        counts[classify_bf(bf10, bf_upper, bf_lower)] += 1

        if track_ci:
            # PairedBayesPropTest.summary.ci_95 is already on the
            # probability scale (Δ = σ(μ+δ_A) − σ(μ)), so we can
            # compare directly to delta_true.
            ci = bb.summary.ci_95
            if ci.lower <= delta_true <= ci.upper:
                ci_covered += 1

        # Frequentist baseline on the same simulated data.
        freq = mcnemar_paired_test(sim.y_A, sim.y_B)
        pvals[i] = freq.p_value

    summary: dict[str, float] = {k: v / n_sim for k, v in counts.items()}
    summary["ci_coverage"] = ci_covered / n_sim if track_ci else float("nan")
    summary["p_A"] = p_A
    summary["p_B"] = p_B
    summary["delta"] = delta_true
    summary["n"] = float(n)
    summary["n_sim"] = float(n_sim)
    return summary, pvals

grid_fixed_n_paired(grid, n, n_sim, seed, *, prior_sigma_delta=1.0, bf_upper=3.0, bf_lower=1.0 / 3.0, n_samples_mc=4000, track_ci=True)

Sweep :func:simulate_fixed_n_paired over an arbitrary (p_A, p_B) grid.

A single deterministic np.random.Generator is created from seed and threaded through every grid cell, so a re-run with the same seed produces bit-identical results.

Parameters:

Name	Type	Description	Default
grid	list[tuple[float, float]]	List of (p_A, p_B) pairs to evaluate.	required
n	int	Number of paired observations (constant across grid cells).	required
n_sim	int	Replicates per grid cell.	required
seed	int	Top-level random seed.	required
prior_sigma_delta	float	See :func:simulate_fixed_n_paired.	1.0
bf_upper	float	See :func:simulate_fixed_n_paired.	3.0
bf_lower	float	See :func:simulate_fixed_n_paired.	1.0 / 3.0
n_samples_mc	int	See :func:simulate_fixed_n_paired.	4000
track_ci	bool	See :func:simulate_fixed_n_paired.	True

Returns:

Type	Description
DataFrame	Tuple (rates_df, p_values). rates_df is a
ndarray	pandas.DataFrame with one row per grid cell and columns
tuple[DataFrame, ndarray]	reject / accept / inconclusive / ci_coverage /
tuple[DataFrame, ndarray]	p_A / p_B / delta / n / n_sim. p_values
tuple[DataFrame, ndarray]	is a (len(grid), n_sim) np.ndarray of McNemar p-values
tuple[DataFrame, ndarray]	— needed to derive the matched calibration α via
tuple[DataFrame, ndarray]	func:matched_calibration_alpha.

Source code in bayesprop/utils/operation_characteristics_paired.py

def grid_fixed_n_paired(
    grid: list[tuple[float, float]],
    n: int,
    n_sim: int,
    seed: int,
    *,
    prior_sigma_delta: float = 1.0,
    bf_upper: float = 3.0,
    bf_lower: float = 1.0 / 3.0,
    n_samples_mc: int = 4000,
    track_ci: bool = True,
) -> tuple[pd.DataFrame, np.ndarray]:
    """Sweep :func:`simulate_fixed_n_paired` over an arbitrary ``(p_A, p_B)`` grid.

    A single deterministic ``np.random.Generator`` is created from
    ``seed`` and threaded through every grid cell, so a re-run with the
    same ``seed`` produces bit-identical results.

    Args:
        grid: List of ``(p_A, p_B)`` pairs to evaluate.
        n: Number of paired observations (constant across grid cells).
        n_sim: Replicates per grid cell.
        seed: Top-level random seed.
        prior_sigma_delta: See :func:`simulate_fixed_n_paired`.
        bf_upper: See :func:`simulate_fixed_n_paired`.
        bf_lower: See :func:`simulate_fixed_n_paired`.
        n_samples_mc: See :func:`simulate_fixed_n_paired`.
        track_ci: See :func:`simulate_fixed_n_paired`.

    Returns:
        Tuple ``(rates_df, p_values)``. ``rates_df`` is a
        ``pandas.DataFrame`` with one row per grid cell and columns
        ``reject`` / ``accept`` / ``inconclusive`` / ``ci_coverage`` /
        ``p_A`` / ``p_B`` / ``delta`` / ``n`` / ``n_sim``. ``p_values``
        is a ``(len(grid), n_sim)`` ``np.ndarray`` of McNemar p-values
        — needed to derive the matched calibration α via
        :func:`matched_calibration_alpha`.
    """
    rng = np.random.default_rng(seed)
    rows: list[dict[str, float]] = []
    pv_stack: list[np.ndarray] = []
    for p_A, p_B in grid:
        summary, pvals = simulate_fixed_n_paired(
            p_A,
            p_B,
            n,
            n_sim,
            rng,
            prior_sigma_delta=prior_sigma_delta,
            bf_upper=bf_upper,
            bf_lower=bf_lower,
            n_samples_mc=n_samples_mc,
            track_ci=track_ci,
        )
        rows.append(summary)
        pv_stack.append(pvals)
    return pd.DataFrame(rows), np.stack(pv_stack)

simulate_sequential_paired(p_A, p_B, n_sim, rng, *, prior_sigma_delta=1.0, bf_upper=3.0, bf_lower=1.0 / 3.0, n_min=50, n_max=600, batch_size=50, n_samples_mc=2000)

Empirical stopping-time distribution of the sequential paired-Laplace test.

For each of n_sim replicates we run a fresh :class:SequentialPairedBayesPropTest, stream batches of size batch_size of paired Bernoulli observations from :func:simulate_paired_scores, and record the per-arm sample size at which the procedure stops. Trials that hit n_max without the BF crossing either threshold are right-censored and classified as "inconclusive".

Parameters:

Name	Type	Description	Default
p_A	float	True success probability for arm A.	required
p_B	float	True success probability for arm B.	required
n_sim	int	Number of independent sequential trajectories.	required
rng	Generator	Pre-seeded NumPy generator.	required
prior_sigma_delta	float	Prior SD on δ_A (logit scale).	1.0
bf_upper	float	Stop for H_1 when BF_10 ≥ bf_upper.	3.0
bf_lower	float	Stop for H_0 when BF_10 ≤ bf_lower.	1.0 / 3.0
n_min	int	Minimum per-arm sample size before any BF-based stop is allowed.	50
n_max	int	Hard cap on per-arm sample size — the trajectory is terminated and the trial reported as censored.	600
batch_size	int	Per-arm batch size delivered to each update().	50
n_samples_mc	int	Posterior Laplace draws per look.	2000

Returns:

Type	Description
dict[str, float]	Dict with the per-trial stopping-time statistics. Keys
dict[str, float]	"p_A", "p_B", "delta" echo the inputs. Keys
dict[str, float]	"median_n", "q05", "q25", "q75", "q95"
dict[str, float]	are quantiles of the per-arm stopping sample size.
dict[str, float]	"frac_censored" is the fraction of trials that hit
dict[str, float]	n_max. "frac_reject" / "frac_accept" are the
dict[str, float]	fractions of trials whose final classification (via
dict[str, float]	func:classify_bf) is "reject" / "accept".

Source code in bayesprop/utils/operation_characteristics_paired.py

def simulate_sequential_paired(
    p_A: float,
    p_B: float,
    n_sim: int,
    rng: np.random.Generator,
    *,
    prior_sigma_delta: float = 1.0,
    bf_upper: float = 3.0,
    bf_lower: float = 1.0 / 3.0,
    n_min: int = 50,
    n_max: int = 600,
    batch_size: int = 50,
    n_samples_mc: int = 2000,
) -> dict[str, float]:
    """Empirical stopping-time distribution of the sequential paired-Laplace test.

    For each of ``n_sim`` replicates we run a fresh
    :class:`SequentialPairedBayesPropTest`, stream batches of size
    ``batch_size`` of paired Bernoulli observations from
    :func:`simulate_paired_scores`, and record the per-arm sample
    size at which the procedure stops. Trials that hit ``n_max``
    without the BF crossing either threshold are *right-censored* and
    classified as ``"inconclusive"``.

    Args:
        p_A: True success probability for arm A.
        p_B: True success probability for arm B.
        n_sim: Number of independent sequential trajectories.
        rng: Pre-seeded NumPy generator.
        prior_sigma_delta: Prior SD on ``δ_A`` (logit scale).
        bf_upper: Stop for ``H_1`` when ``BF_10 ≥ bf_upper``.
        bf_lower: Stop for ``H_0`` when ``BF_10 ≤ bf_lower``.
        n_min: Minimum per-arm sample size before any BF-based stop
            is allowed.
        n_max: Hard cap on per-arm sample size — the trajectory is
            terminated and the trial reported as censored.
        batch_size: Per-arm batch size delivered to each ``update()``.
        n_samples_mc: Posterior Laplace draws per look.

    Returns:
        Dict with the per-trial stopping-time statistics. Keys
        ``"p_A"``, ``"p_B"``, ``"delta"`` echo the inputs. Keys
        ``"median_n"``, ``"q05"``, ``"q25"``, ``"q75"``, ``"q95"``
        are quantiles of the per-arm stopping sample size.
        ``"frac_censored"`` is the fraction of trials that hit
        ``n_max``. ``"frac_reject"`` / ``"frac_accept"`` are the
        fractions of trials whose final classification (via
        :func:`classify_bf`) is ``"reject"`` / ``"accept"``.
    """
    max_looks = n_max // batch_size

    stop_n: list[int] = []
    censored: list[bool] = []
    decisions: list[BFDecision] = []

    for _ in range(n_sim):
        seq = SequentialPairedBayesPropTest(
            prior_sigma_delta=prior_sigma_delta,
            bf_upper=bf_upper,
            bf_lower=bf_lower,
            n_min=n_min,
            n_max=n_max,
            decision_rule="bayes_factor",
            n_samples=n_samples_mc,
        )
        last = None
        for _look in range(max_looks):
            batch = simulate_paired_scores(
                N=batch_size,
                theta_A=p_A,
                theta_B=p_B,
                sigma_theta=0.0,
                rng=rng,
            )
            last = seq.update(batch.y_A, batch.y_B)
            if seq.stopped:
                break
        assert last is not None, "max_looks=0 would skip the inner loop"

        stop_n.append(min(last.n_A, last.n_B))
        censored.append("n_max" in (last.stop_reason or "n_max reached"))

        bf10_final = (
            last.decision.bayes_factor.BF_10 if last.decision.bayes_factor else 1.0
        )
        decisions.append(classify_bf(bf10_final, bf_upper, bf_lower))

    arr = np.asarray(stop_n)
    decision_arr = np.asarray(decisions)
    return {
        "p_A": p_A,
        "p_B": p_B,
        "delta": p_A - p_B,
        "median_n": float(np.median(arr)),
        "q05": float(np.quantile(arr, 0.05)),
        "q25": float(np.quantile(arr, 0.25)),
        "q75": float(np.quantile(arr, 0.75)),
        "q95": float(np.quantile(arr, 0.95)),
        "frac_censored": float(np.mean(censored)),
        "frac_reject": float(np.mean(decision_arr == "reject")),
        "frac_accept": float(np.mean(decision_arr == "accept")),
    }