Utilities

General-purpose helpers used across the package. Grouped by purpose:

• Data-generating processes — simulate_nonpaired_scores, simulate_paired_scores: generate paired and non-paired binary outcomes from configurable true effects. Used by the BFDA harness and the operating-characteristics notebooks.
• Frequentist baseline tests — fisher_exact_nonpaired_test, mcnemar_paired_test: classical tests that serve as the frequentist baseline in the operating-characteristics analyses (non-paired, paired).
• Bayes Factor Design Analysis (BFDA) — bfda_simulate, bfda_power_curve, find_n_for_power, plot_bfda_power, plot_bfda_sensitivity: Monte-Carlo sample-size planning under a Bayes-factor decision rule. See the sample-size planning guide.
• Decision helpers — bf10_to_ph0: convert a Bayes factor BF₁₀ into the posterior probability of the null under a given prior π_H₀.

utils

binarize_if_needed(y, threshold=0.5, *, name='y', verbose=False)

Return a 0/1 float64 array, binarising continuous inputs in [0, 1].

The package's models all operate on binary 0/1 data, but users often have continuous scores in [0, 1] (e.g. predicted probabilities, classifier confidences). This helper provides a single, consistent coercion across the non-paired and paired classes:

• If every entry of y is already exactly 0 or 1 (after casting to float), the array is returned unchanged.
• Otherwise, values strictly outside [0, 1] or any NaN raise :class:ValueError immediately — we never silently clip or truncate.
• Otherwise, y is binarised as (y >= threshold).astype(float).

Parameters:

Name	Type	Description	Default
y	ArrayLike	Input scores. Any array-like accepted by :func:numpy.asarray.	required
threshold	float	Cut-off for binarisation. Values ≥ threshold map to 1.0, the rest to 0.0. Defaults to 0.5.	0.5
name	str	Name of the variable in error / warning messages (e.g. "y_A_obs"). Purely cosmetic.	'y'
verbose	bool	If True, print a one-line notice when binarisation actually fires.	False

Returns:

Type	Description
ndarray	np.ndarray of float64 containing only 0.0 and 1.0.

Raises:

Type	Description
ValueError	If y contains NaN or any value outside [0, 1].

Source code in bayesprop/utils/utils.py

def binarize_if_needed(
    y: npt.ArrayLike,
    threshold: float = 0.5,
    *,
    name: str = "y",
    verbose: bool = False,
) -> np.ndarray:
    """Return a 0/1 ``float64`` array, binarising continuous inputs in ``[0, 1]``.

    The package's models all operate on binary 0/1 data, but users often
    have continuous scores in ``[0, 1]`` (e.g. predicted probabilities,
    classifier confidences). This helper provides a single, consistent
    coercion across the non-paired and paired classes:

    * If every entry of ``y`` is already exactly ``0`` or ``1`` (after
      casting to float), the array is returned unchanged.
    * Otherwise, values strictly outside ``[0, 1]`` or any ``NaN`` raise
      :class:`ValueError` immediately — we never silently clip or truncate.
    * Otherwise, ``y`` is binarised as ``(y >= threshold).astype(float)``.

    Args:
        y: Input scores. Any array-like accepted by :func:`numpy.asarray`.
        threshold: Cut-off for binarisation. Values ``≥ threshold`` map to
            ``1.0``, the rest to ``0.0``. Defaults to ``0.5``.
        name: Name of the variable in error / warning messages
            (e.g. ``"y_A_obs"``). Purely cosmetic.
        verbose: If ``True``, print a one-line notice when binarisation
            actually fires.

    Returns:
        ``np.ndarray`` of ``float64`` containing only ``0.0`` and ``1.0``.

    Raises:
        ValueError: If ``y`` contains ``NaN`` or any value outside
            ``[0, 1]``.
    """
    arr = np.asarray(y, dtype=float)
    if arr.size == 0:
        return arr

    # Fast path: already binary.
    is_binary = np.all((arr == 0.0) | (arr == 1.0))
    if is_binary:
        return arr

    # Validate range before binarising — silent truncation / clipping is
    # the failure mode we are protecting users against.
    if np.any(np.isnan(arr)):
        raise ValueError(f"{name} contains NaN values; cannot binarise.")
    if np.any((arr < 0.0) | (arr > 1.0)):
        raise ValueError(
            f"{name} contains values outside [0, 1] "
            f"(observed range: [{arr.min():.4g}, {arr.max():.4g}]); "
            f"refusing to binarise. Please rescale your scores to [0, 1] "
            f"or supply 0/1 inputs directly."
        )

    if verbose:
        print(
            f"Warning: non-binary {name} detected, binarising at threshold {threshold}"
        )
    return (arr >= threshold).astype(float)

simulate_nonpaired_scores(N=200, theta_A=0.75, theta_B=0.6, seed=0, rng=None)

Simulate independent binary outcomes for a non-paired A/B test.

Each group is sampled independently from a Bernoulli distribution with the specified success probability.

Parameters:

Name	Type	Description	Default
N	int	Number of observations per group.	200
theta_A	float	True success probability for model A.	0.75
theta_B	float	True success probability for model B.	0.6
seed	int	Random seed for reproducibility.	0
rng	Generator | None	Optional pre-seeded RNG; if provided, seed is ignored.	None

Returns:

Type	Description
NonPairedSimResult	class:NonPairedSimResult with fields y_A, y_B,
NonPairedSimResult	theta_A, theta_B, and true_params.

Source code in bayesprop/utils/utils.py

def simulate_nonpaired_scores(
    N: int = 200,
    theta_A: float = 0.75,
    theta_B: float = 0.60,
    seed: int = 0,
    rng: np.random.Generator | None = None,
) -> NonPairedSimResult:
    """Simulate independent binary outcomes for a non-paired A/B test.

    Each group is sampled independently from a Bernoulli distribution
    with the specified success probability.

    Args:
        N: Number of observations per group.
        theta_A: True success probability for model A.
        theta_B: True success probability for model B.
        seed: Random seed for reproducibility.
        rng: Optional pre-seeded RNG; if provided, *seed* is ignored.

    Returns:
        :class:`NonPairedSimResult` with fields ``y_A``, ``y_B``,
        ``theta_A``, ``theta_B``, and ``true_params``.
    """
    if rng is None:
        rng = np.random.default_rng(seed)
    y_A = rng.binomial(1, theta_A, size=N).astype(float)
    y_B = rng.binomial(1, theta_B, size=N).astype(float)
    return NonPairedSimResult(
        y_A=y_A,
        y_B=y_B,
        theta_A=theta_A,
        theta_B=theta_B,
        true_params=NonPairedTrueParams(N=N, theta_A=theta_A, theta_B=theta_B),
    )

simulate_paired_scores(N=200, theta_A=0.75, theta_B=0.6, sigma_theta=0.0, seed=0, rng=None)

Simulate paired binary outcomes from a logistic DGP.

Matches the paired model: y_A ~ Bern(σ(μ + δ_A)), y_B ~ Bern(σ(μ)), where μ = logit(theta_B) and δ_A = logit(theta_A) − logit(theta_B).

When sigma_theta > 0 each item i additionally receives a random effect ε_i ~ N(0, sigma_theta) so that θ_i = μ + ε_i (useful for more realistic BFDA simulations).

Parameters:

Name	Type	Description	Default
N	int	Number of paired observations.	200
theta_A	float	True success probability for model A.	0.75
theta_B	float	True success probability for model B.	0.6
sigma_theta	float	SD of optional per-item random effects (0 = fixed effects matching the model).	0.0
seed	int	Random seed for reproducibility.	0
rng	Generator | None	Optional pre-seeded RNG; if provided, seed is ignored.	None

Returns:

Type	Description
PairedSimResult	class:PairedSimResult with fields y_A, y_B,
PairedSimResult	theta_A, theta_B, p_A_true, p_B_true,
PairedSimResult	theta_true, and true_params.

Source code in bayesprop/utils/utils.py

def simulate_paired_scores(
    N: int = 200,
    theta_A: float = 0.75,
    theta_B: float = 0.60,
    sigma_theta: float = 0.0,
    seed: int = 0,
    rng: np.random.Generator | None = None,
) -> PairedSimResult:
    """Simulate paired binary outcomes from a logistic DGP.

    Matches the paired model: ``y_A ~ Bern(σ(μ + δ_A))``,
    ``y_B ~ Bern(σ(μ))``, where ``μ = logit(theta_B)`` and
    ``δ_A = logit(theta_A) − logit(theta_B)``.

    When ``sigma_theta > 0`` each item *i* additionally receives a
    random effect ``ε_i ~ N(0, sigma_theta)`` so that
    ``θ_i = μ + ε_i`` (useful for more realistic BFDA simulations).

    Args:
        N: Number of paired observations.
        theta_A: True success probability for model A.
        theta_B: True success probability for model B.
        sigma_theta: SD of optional per-item random effects
            (``0`` = fixed effects matching the model).
        seed: Random seed for reproducibility.
        rng: Optional pre-seeded RNG; if provided, *seed* is ignored.

    Returns:
        :class:`PairedSimResult` with fields ``y_A``, ``y_B``,
        ``theta_A``, ``theta_B``, ``p_A_true``, ``p_B_true``,
        ``theta_true``, and ``true_params``.
    """
    if rng is None:
        rng = np.random.default_rng(seed)

    mu = float(np.log(theta_B / (1.0 - theta_B)))
    delta_A = float(np.log(theta_A / (1.0 - theta_A))) - mu

    theta_true = mu + rng.normal(0.0, sigma_theta, size=N)
    p_A = _sigmoid(theta_true + delta_A)
    p_B = _sigmoid(theta_true)
    y_A = rng.binomial(1, p_A)
    y_B = rng.binomial(1, p_B)
    return PairedSimResult(
        y_A=y_A,
        y_B=y_B,
        theta_A=theta_A,
        theta_B=theta_B,
        p_A_true=p_A,
        p_B_true=p_B,
        theta_true=theta_true,
        true_params=PairedTrueParams(
            N=N,
            theta_A=theta_A,
            theta_B=theta_B,
            mu=mu,
            delta_A=delta_A,
            sigma_theta=sigma_theta,
        ),
    )

fisher_exact_nonpaired_test(y_A, y_B, *, alternative='two-sided')

Fisher's exact test for two independent proportions.

Frequentist baseline that mirrors the data contract of :class:bayesprop.resources.bayes_nonpaired.NonPairedBayesPropTest but returns a classical p-value instead of a Bayes factor. Useful as a calibration reference for operating-characteristic / Bayes-factor design analyses (e.g. overlaying a frequentist power curve onto the Bayesian three-way OC plot).

The test is exact under H₀: p_A = p_B (no large-sample approximation), so it remains valid at small n and at boundary success rates p near 0 or 1, where the Wald Z and Pearson chi-squared tests are unreliable.

Parameters:

Name	Type	Description	Default
y_A	ArrayLike	Binary outcomes (0/1) for group A. Non-binary inputs are rejected — binarise upstream if your scores are continuous.	required
y_B	ArrayLike	Binary outcomes (0/1) for group B.	required
alternative	Literal['two-sided', 'greater', 'less']	"two-sided" tests p_A ≠ p_B; "greater" tests p_A > p_B; "less" tests p_A < p_B.	'two-sided'

Returns:

Type	Description
FrequentistTestResult	class:FrequentistTestResult with the p-value, sample odds
FrequentistTestResult	ratio, and cell counts.

Raises:

Type	Description
ValueError	If either array contains values outside {0, 1}.

Source code in bayesprop/utils/utils.py

def fisher_exact_nonpaired_test(
    y_A: npt.ArrayLike,
    y_B: npt.ArrayLike,
    *,
    alternative: Literal["two-sided", "greater", "less"] = "two-sided",
) -> FrequentistTestResult:
    """Fisher's exact test for two independent proportions.

    Frequentist baseline that mirrors the data contract of
    :class:`bayesprop.resources.bayes_nonpaired.NonPairedBayesPropTest`
    but returns a classical p-value instead of a Bayes factor. Useful
    as a calibration reference for operating-characteristic /
    Bayes-factor design analyses (e.g. overlaying a frequentist power
    curve onto the Bayesian three-way OC plot).

    The test is exact under H₀: ``p_A = p_B`` (no large-sample
    approximation), so it remains valid at small ``n`` and at boundary
    success rates ``p`` near 0 or 1, where the Wald Z and Pearson
    chi-squared tests are unreliable.

    Args:
        y_A: Binary outcomes (0/1) for group A. Non-binary inputs are
            rejected — binarise upstream if your scores are continuous.
        y_B: Binary outcomes (0/1) for group B.
        alternative: ``"two-sided"`` tests ``p_A ≠ p_B``;
            ``"greater"`` tests ``p_A > p_B``; ``"less"`` tests
            ``p_A < p_B``.

    Returns:
        :class:`FrequentistTestResult` with the p-value, sample odds
        ratio, and cell counts.

    Raises:
        ValueError: If either array contains values outside ``{0, 1}``.
    """
    arr_A = np.asarray(y_A, dtype=float)
    arr_B = np.asarray(y_B, dtype=float)
    for name, arr in (("y_A", arr_A), ("y_B", arr_B)):
        if arr.size and not np.all((arr == 0.0) | (arr == 1.0)):
            raise ValueError(
                f"{name} must contain only 0/1 values; got non-binary entries"
            )

    n_A = int(arr_A.size)
    n_B = int(arr_B.size)
    sA = int(arr_A.sum())
    sB = int(arr_B.sum())

    # 2x2 contingency table:
    #               success   failure
    #   group A:      sA      n_A - sA
    #   group B:      sB      n_B - sB
    table = np.array([[sA, n_A - sA], [sB, n_B - sB]], dtype=int)
    odds_ratio_raw, p_value = fisher_exact(table, alternative=alternative)
    odds_ratio: float | None = (
        float(odds_ratio_raw) if np.isfinite(odds_ratio_raw) else None
    )

    return FrequentistTestResult(
        p_value=float(p_value),
        odds_ratio=odds_ratio,
        successes_A=sA,
        successes_B=sB,
        n_A=n_A,
        n_B=n_B,
        alternative=alternative,
        test="fisher_exact",
    )

mcnemar_paired_test(y_A, y_B, *, alternative='two-sided', exact=None)

McNemar's test for two paired (matched) proportions.

Frequentist baseline that mirrors the data contract of :class:bayesprop.resources.bayes_paired_laplace.PairedBayesPropTest but returns a classical p-value instead of a Bayes factor. Useful as a calibration reference for paired operating-characteristic / Bayes-factor design analyses.

Given two binary arrays y_A and y_B of equal length, the test conditions on the discordant pairs:

• b = #{i : y_A[i] = 1, y_B[i] = 0}
• c = #{i : y_A[i] = 0, y_B[i] = 1}

Under H_0: p_A = p_B and conditional on b + c, b follows Binomial(b + c, 0.5). The two-sided p-value is computed as 2 · min(P(B ≤ b), P(B ≥ b)) clipped to [0, 1] from an exact binomial calculation when b + c is small (default threshold 25, matching :func:scipy.stats.contingency.mcnemar's convention), and from the standard chi-squared approximation otherwise.

Parameters:

Name	Type	Description	Default
y_A	ArrayLike	Binary outcomes (0/1) for arm A. Non-binary inputs raise.	required
y_B	ArrayLike	Binary outcomes (0/1) for arm B. Must be the same length as y_A (pairs are aligned by index).	required
alternative	Literal['two-sided', 'greater', 'less']	"two-sided" tests p_A ≠ p_B; "greater" tests p_A > p_B (one-sided on b); "less" tests p_A < p_B.	'two-sided'
exact	bool | None	If True, force the exact binomial p-value; if False, force the chi-squared approximation; if None (default), the exact test is used when b + c <= 25.	None

Returns:

Type	Description
FrequentistTestResult	class:FrequentistTestResult with the p-value, the discordant
FrequentistTestResult	odds ratio b / c (or None if c == 0), and the
FrequentistTestResult	marginal counts.

Raises:

Type	Description
ValueError	If either array contains non-binary values or the two arrays have different lengths.

Source code in bayesprop/utils/utils.py

def mcnemar_paired_test(
    y_A: npt.ArrayLike,
    y_B: npt.ArrayLike,
    *,
    alternative: Literal["two-sided", "greater", "less"] = "two-sided",
    exact: bool | None = None,
) -> FrequentistTestResult:
    """McNemar's test for two paired (matched) proportions.

    Frequentist baseline that mirrors the data contract of
    :class:`bayesprop.resources.bayes_paired_laplace.PairedBayesPropTest`
    but returns a classical p-value instead of a Bayes factor. Useful
    as a calibration reference for paired operating-characteristic /
    Bayes-factor design analyses.

    Given two binary arrays ``y_A`` and ``y_B`` of equal length, the
    test conditions on the **discordant pairs**:

    * ``b = #{i : y_A[i] = 1, y_B[i] = 0}``
    * ``c = #{i : y_A[i] = 0, y_B[i] = 1}``

    Under ``H_0: p_A = p_B`` and conditional on ``b + c``, ``b`` follows
    ``Binomial(b + c, 0.5)``. The two-sided p-value is computed as
    ``2 · min(P(B ≤ b), P(B ≥ b))`` clipped to ``[0, 1]`` from an
    exact binomial calculation when ``b + c`` is small (default
    threshold 25, matching :func:`scipy.stats.contingency.mcnemar`'s
    convention), and from the standard chi-squared approximation
    otherwise.

    Args:
        y_A: Binary outcomes (0/1) for arm A. Non-binary inputs raise.
        y_B: Binary outcomes (0/1) for arm B. Must be the same length
            as ``y_A`` (pairs are aligned by index).
        alternative: ``"two-sided"`` tests ``p_A ≠ p_B``; ``"greater"``
            tests ``p_A > p_B`` (one-sided on ``b``); ``"less"`` tests
            ``p_A < p_B``.
        exact: If ``True``, force the exact binomial p-value; if
            ``False``, force the chi-squared approximation; if ``None``
            (default), the exact test is used when ``b + c <= 25``.

    Returns:
        :class:`FrequentistTestResult` with the p-value, the discordant
        odds ratio ``b / c`` (or ``None`` if ``c == 0``), and the
        marginal counts.

    Raises:
        ValueError: If either array contains non-binary values or the
            two arrays have different lengths.
    """
    arr_A = np.asarray(y_A, dtype=float)
    arr_B = np.asarray(y_B, dtype=float)
    if arr_A.shape != arr_B.shape:
        raise ValueError(
            f"y_A and y_B must have the same length for a paired test; "
            f"got shapes {arr_A.shape} and {arr_B.shape}"
        )
    for name, arr in (("y_A", arr_A), ("y_B", arr_B)):
        if arr.size and not np.all((arr == 0.0) | (arr == 1.0)):
            raise ValueError(
                f"{name} must contain only 0/1 values; got non-binary entries"
            )

    n = int(arr_A.size)
    # Discordant pair counts. McNemar conditions on these and ignores
    # the (a, d) cells (concordant pairs carry no information about the
    # marginal-rate difference).
    b = int(np.sum((arr_A == 1.0) & (arr_B == 0.0)))
    c = int(np.sum((arr_A == 0.0) & (arr_B == 1.0)))
    n_disc = b + c

    # Default: exact for small discordant counts, chi² otherwise.
    use_exact = exact if exact is not None else (n_disc <= 25)

    if n_disc == 0:
        # No information: by convention p = 1 (no evidence against H_0).
        p_value = 1.0
    elif use_exact:
        # binomtest gives the appropriate p-value for each alternative.
        # For McNemar, "greater" / "less" map to b being large / small
        # under Binomial(b+c, 0.5).
        bt = binomtest(
            b,
            n_disc,
            p=0.5,
            alternative=(
                "two-sided"
                if alternative == "two-sided"
                else ("greater" if alternative == "greater" else "less")
            ),
        )
        p_value = float(bt.pvalue)
    else:
        from scipy.stats import chi2

        chi2_stat = (b - c) ** 2 / n_disc
        p_two_sided = float(1.0 - chi2.cdf(chi2_stat, df=1))
        if alternative == "two-sided":
            p_value = p_two_sided
        else:
            # Split the two-sided p-value by direction of the observed
            # imbalance; this is the standard one-sided McNemar.
            half = 0.5 * p_two_sided
            if (alternative == "greater" and b > c) or (
                alternative == "less" and b < c
            ):
                p_value = half
            else:
                p_value = 1.0 - half

    odds_ratio: float | None = float(b) / c if c > 0 else None

    return FrequentistTestResult(
        p_value=p_value,
        odds_ratio=odds_ratio,
        successes_A=int(arr_A.sum()),
        successes_B=int(arr_B.sum()),
        n_A=n,
        n_B=n,
        alternative=alternative,
        test="mcnemar_exact" if use_exact else "mcnemar_chi2",
    )

bfda_simulate(data_generator, decision_fn, sample_sizes, n_sim=500, seed=42)

Generic BFDA engine -- works with any data-generating process and decision rule.

Parameters:

Name	Type	Description	Default
data_generator	Callable[[Generator, int], tuple[ndarray, ndarray]]	Callable(rng, n) -> (y_A, y_B). Generates one simulated dataset of size n per group using the provided RNG.	required
decision_fn	Callable[[ndarray, ndarray], bool]	Callable(y_A, y_B) -> bool. Returns True when the simulated dataset produces a "decisive" result (i.e. rejects H0).	required
sample_sizes	list[int]	List of per-group sample sizes to evaluate.	required
n_sim	int	Number of simulated datasets per sample size.	500
seed	int	Random seed for reproducibility.	42

Returns:

Type	Description
dict[int, float]	Dictionary mapping sample size -> P(decisive outcome).

Source code in bayesprop/utils/utils.py

def bfda_simulate(
    data_generator: Callable[[np.random.Generator, int], tuple[np.ndarray, np.ndarray]],
    decision_fn: Callable[[np.ndarray, np.ndarray], bool],
    sample_sizes: list[int],
    n_sim: int = 500,
    seed: int = 42,
) -> dict[int, float]:
    """Generic BFDA engine -- works with any data-generating process and decision rule.

    Args:
        data_generator: Callable(rng, n) -> (y_A, y_B). Generates one simulated
            dataset of size *n* per group using the provided RNG.
        decision_fn: Callable(y_A, y_B) -> bool. Returns ``True`` when the
            simulated dataset produces a "decisive" result (i.e. rejects H0).
        sample_sizes: List of per-group sample sizes to evaluate.
        n_sim: Number of simulated datasets per sample size.
        seed: Random seed for reproducibility.

    Returns:
        Dictionary mapping sample size -> P(decisive outcome).
    """
    rng = np.random.default_rng(seed)

    power: dict[int, float] = {}
    for n in sample_sizes:
        decisive_count = sum(decision_fn(*data_generator(rng, n)) for _ in range(n_sim))
        power[n] = decisive_count / n_sim

    return power

bf10_to_ph0(bf_10, prior_H0=0.5)

Convert BF_10 to posterior probability of H0.

Parameters:

Name	Type	Description	Default
bf_10	float	Bayes factor in favour of H1.	required
prior_H0	float	Prior probability of H0.	0.5

Returns:

Type	Description
float	P(H0 | data).

Source code in bayesprop/utils/utils.py

def bf10_to_ph0(bf_10: float, prior_H0: float = 0.5) -> float:
    """Convert BF_10 to posterior probability of H0.

    Args:
        bf_10: Bayes factor in favour of H1.
        prior_H0: Prior probability of H0.

    Returns:
        P(H0 | data).
    """
    bf_01 = 1.0 / max(bf_10, 1e-300)
    return (bf_01 * prior_H0) / (bf_01 * prior_H0 + (1 - prior_H0))

bfda_power_curve(theta_A_true, theta_B_true, sample_sizes, design='nonpaired', decision_rule='bayes_factor', *, bf_threshold=3.0, ph0_threshold=0.05, prior_H0=0.5, rope=(-0.02, 0.02), ci_mass=0.95, n_sim=500, seed=42, alpha0=1.0, beta0=1.0, sigma_theta=2.0, prior_sigma_delta=1.0, prior_sigma_mu=2.0, n_iter=1000, burn_in=300, n_chains=2)

Unified Bayes Factor Design Analysis for any design × decision-rule.

Simulates datasets under a known effect and estimates the probability that a given Bayesian decision rule will reject H0 as a function of sample size (i.e. Bayesian "power").

Supported combinations:

+----------------+----------------+--------------------+--------+ | design | decision_rule | key threshold | fast? | +================+================+====================+========+ | nonpaired | bayes_factor | bf_threshold | yes | | nonpaired | posterior_null | ph0_threshold | yes | | nonpaired | rope | rope | medium | | paired | bayes_factor | bf_threshold | slow | | paired | posterior_null | ph0_threshold | slow | | paired | rope | rope | slow | +----------------+----------------+--------------------+--------+

Parameters:

Name	Type	Description	Default
theta_A_true	float	Assumed true success rate for model A.	required
theta_B_true	float	Assumed true success rate for model B.	required
sample_sizes	list[int]	List of per-group sample sizes to evaluate.	required
design	str	"nonpaired" or "paired".	'nonpaired'
decision_rule	DecisionRuleType	"bayes_factor", "posterior_null", or "rope".	'bayes_factor'
bf_threshold	float	BF_10 threshold for decisive evidence (bayes_factor).	3.0
ph0_threshold	float	Reject H0 when P(H0|data) < this (posterior_null).	0.05
prior_H0	float	Prior probability of H0 (posterior_null).	0.5
rope	tuple[float, float]	(lower, upper) bounds of the ROPE (rope).	(-0.02, 0.02)
ci_mass	float	Credible interval mass for ROPE analysis (rope).	0.95
n_sim	int	Number of simulated datasets per sample size.	500
seed	int	Random seed for reproducibility.	42
alpha0	float	Prior Beta alpha parameter (non-paired only).	1.0
beta0	float	Prior Beta beta parameter (non-paired only).	1.0
sigma_theta	float	SD of the shared latent item effect (paired DGP).	2.0
prior_sigma_delta	float	SD of N(0, σ) prior on delta_A (paired only).	1.0
prior_sigma_mu	float	SD of N(0, σ) prior on mu (paired only).	2.0
n_iter	int	Total Gibbs iterations per chain (paired only).	1000
burn_in	int	Warm-up iterations per chain (paired only).	300
n_chains	int	Number of MCMC chains per dataset (paired only).	2

Returns:

Type	Description
dict[int, float]	Dictionary mapping sample size -> P(decisive outcome).

Source code in bayesprop/utils/utils.py

def bfda_power_curve(
    theta_A_true: float,
    theta_B_true: float,
    sample_sizes: list[int],
    design: str = "nonpaired",
    decision_rule: DecisionRuleType = "bayes_factor",
    *,
    # BF thresholds
    bf_threshold: float = 3.0,
    # P(H0) thresholds
    ph0_threshold: float = 0.05,
    prior_H0: float = 0.5,
    # ROPE thresholds
    rope: tuple[float, float] = (-0.02, 0.02),
    ci_mass: float = 0.95,
    # Simulation
    n_sim: int = 500,
    seed: int = 42,
    # Non-paired model priors
    alpha0: float = 1.0,
    beta0: float = 1.0,
    # Paired DGP
    sigma_theta: float = 2.0,
    # Paired model priors / MCMC
    prior_sigma_delta: float = 1.0,
    prior_sigma_mu: float = 2.0,
    n_iter: int = 1000,
    burn_in: int = 300,
    n_chains: int = 2,
) -> dict[int, float]:
    """Unified Bayes Factor Design Analysis for any design × decision-rule.

    Simulates datasets under a known effect and estimates the probability
    that a given Bayesian decision rule will reject H0 as a function of
    sample size (i.e. Bayesian "power").

    Supported combinations:

    +----------------+----------------+--------------------+--------+
    | design         | decision_rule  | key threshold      | fast?  |
    +================+================+====================+========+
    | ``nonpaired``  | bayes_factor   | ``bf_threshold``   | yes    |
    | ``nonpaired``  | posterior_null | ``ph0_threshold``  | yes    |
    | ``nonpaired``  | rope           | ``rope``           | medium |
    | ``paired``     | bayes_factor   | ``bf_threshold``   | slow   |
    | ``paired``     | posterior_null | ``ph0_threshold``  | slow   |
    | ``paired``     | rope           | ``rope``           | slow   |
    +----------------+----------------+--------------------+--------+

    Args:
        theta_A_true: Assumed true success rate for model A.
        theta_B_true: Assumed true success rate for model B.
        sample_sizes: List of per-group sample sizes to evaluate.
        design: ``"nonpaired"`` or ``"paired"``.
        decision_rule: ``"bayes_factor"``, ``"posterior_null"``, or ``"rope"``.
        bf_threshold: BF_10 threshold for decisive evidence (``bayes_factor``).
        ph0_threshold: Reject H0 when P(H0|data) < this (``posterior_null``).
        prior_H0: Prior probability of H0 (``posterior_null``).
        rope: (lower, upper) bounds of the ROPE (``rope``).
        ci_mass: Credible interval mass for ROPE analysis (``rope``).
        n_sim: Number of simulated datasets per sample size.
        seed: Random seed for reproducibility.
        alpha0: Prior Beta alpha parameter (non-paired only).
        beta0: Prior Beta beta parameter (non-paired only).
        sigma_theta: SD of the shared latent item effect (paired DGP).
        prior_sigma_delta: SD of N(0, σ) prior on delta_A (paired only).
        prior_sigma_mu: SD of N(0, σ) prior on mu (paired only).
        n_iter: Total Gibbs iterations per chain (paired only).
        burn_in: Warm-up iterations per chain (paired only).
        n_chains: Number of MCMC chains per dataset (paired only).

    Returns:
        Dictionary mapping sample size -> P(decisive outcome).
    """
    # Build data generator
    if design == "nonpaired":
        data_gen = _make_nonpaired_generator(theta_A_true, theta_B_true)
    elif design == "paired":
        data_gen = _make_paired_generator(
            theta_A_true, theta_B_true, sigma_theta=sigma_theta
        )
    else:
        raise ValueError(f"Unknown design: {design!r}. Use 'nonpaired' or 'paired'.")

    # Build decision function
    decide = _make_decision_fn(
        design=design,
        decision_rule=decision_rule,
        bf_threshold=bf_threshold,
        ph0_threshold=ph0_threshold,
        prior_H0=prior_H0,
        rope=rope,
        ci_mass=ci_mass,
        alpha0=alpha0,
        beta0=beta0,
        prior_sigma_delta=prior_sigma_delta,
        prior_sigma_mu=prior_sigma_mu,
        n_iter=n_iter,
        burn_in=burn_in,
        n_chains=n_chains,
        seed=seed,
    )

    return bfda_simulate(
        data_generator=data_gen,
        decision_fn=decide,
        sample_sizes=sample_sizes,
        n_sim=n_sim,
        seed=seed,
    )

find_n_for_power(power_curve, target_power=0.8)

Interpolate the sample size needed to achieve a target power level.

Parameters:

Name	Type	Description	Default
power_curve	dict[int, float]	Dictionary mapping sample size -> power (from BFDA).	required
target_power	float	Desired power level (default 0.80).	0.8

Returns:

Type	Description
float | None	Interpolated sample size, or None if target is outside the range.

Source code in bayesprop/utils/utils.py

def find_n_for_power(
    power_curve: dict[int, float],
    target_power: float = 0.80,
) -> float | None:
    """Interpolate the sample size needed to achieve a target power level.

    Args:
        power_curve: Dictionary mapping sample size -> power (from BFDA).
        target_power: Desired power level (default 0.80).

    Returns:
        Interpolated sample size, or ``None`` if target is outside the range.
    """
    ns = list(power_curve.keys())
    ps = list(power_curve.values())

    if min(ps) >= target_power:
        return float(ns[0])
    if max(ps) < target_power:
        return None

    f_interp = interp1d(ps, ns, kind="linear")
    return float(f_interp(target_power))

plot_bfda_power(power_curve, theta_A_true, theta_B_true, bf_threshold=3.0, target_power=0.8, title=None, ax=None)

Plot a BFDA power curve with 80%/95% reference lines.

Parameters:

Name	Type	Description	Default
power_curve	dict[int, float]	Dictionary mapping sample size -> power.	required
theta_A_true	float	Assumed true rate for model A (for title).	required
theta_B_true	float	Assumed true rate for model B (for title).	required
bf_threshold	float	BF_10 threshold used (for y-axis label).	3.0
target_power	float	Power level to highlight via interpolation.	0.8
title	str | None	Optional custom title.	None
ax	Axes | None	Optional matplotlib Axes to plot on.	None

Returns:

Type	Description
Figure	The matplotlib Figure.

Source code in bayesprop/utils/utils.py

def plot_bfda_power(
    power_curve: dict[int, float],
    theta_A_true: float,
    theta_B_true: float,
    bf_threshold: float = 3.0,
    target_power: float = 0.80,
    title: str | None = None,
    ax: plt.Axes | None = None,
) -> plt.Figure:
    """Plot a BFDA power curve with 80%/95% reference lines.

    Args:
        power_curve: Dictionary mapping sample size -> power.
        theta_A_true: Assumed true rate for model A (for title).
        theta_B_true: Assumed true rate for model B (for title).
        bf_threshold: BF_10 threshold used (for y-axis label).
        target_power: Power level to highlight via interpolation.
        title: Optional custom title.
        ax: Optional matplotlib Axes to plot on.

    Returns:
        The matplotlib Figure.
    """
    if ax is None:
        fig, ax = plt.subplots(figsize=(9, 5))
    else:
        fig = ax.get_figure()

    ns = list(power_curve.keys())
    ps = list(power_curve.values())

    ax.plot(ns, ps, "o-", color="#2196F3", linewidth=2, markersize=7)
    ax.axhline(
        0.80, color="gray", linestyle="--", linewidth=1, alpha=0.7, label="80% power"
    )
    ax.axhline(
        0.95, color="gray", linestyle=":", linewidth=1, alpha=0.7, label="95% power"
    )

    n_target = find_n_for_power(power_curve, target_power)
    if n_target is not None:
        ax.axvline(
            n_target,
            color="#E91E63",
            linestyle="-.",
            linewidth=1.5,
            alpha=0.7,
            label=f"n ~ {n_target:.0f} for {target_power:.0%} power",
        )

    delta = theta_A_true - theta_B_true
    if title is None:
        title = f"BFDA Power Curve -- Delta = {delta:.3f} (theta_A={theta_A_true:.2f}, theta_B={theta_B_true:.2f})"
    ax.set_xlabel("Sample size per group (n)")
    ax.set_ylabel(f"P(BF_10 > {bf_threshold:.0f})")
    ax.set_title(title, fontsize=12, fontweight="bold")
    ax.set_ylim(-0.02, 1.02)
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
    fig.tight_layout()

    return fig

plot_bfda_sensitivity(theta_A_true, theta_B_true, sample_sizes, thresholds=None, n_sim=500, seed=42, design='nonpaired', title=None, ax=None, **kwargs)

Plot BFDA power curves for multiple BF_10 thresholds.

Works for both paired and non-paired designs.

Parameters:

Name	Type	Description	Default
theta_A_true	float	Assumed true success rate for model A.	required
theta_B_true	float	Assumed true success rate for model B.	required
sample_sizes	list[int]	List of per-group sample sizes to evaluate.	required
thresholds	list[float] | None	BF_10 thresholds to compare (default: [3, 6, 10]).	None
n_sim	int	Number of simulated datasets per sample size.	500
seed	int	Random seed for reproducibility.	42
design	str	"nonpaired" or "paired".	'nonpaired'
title	str | None	Optional custom title.	None
ax	Axes | None	Optional matplotlib Axes to plot on.	None
**kwargs	Any	Extra arguments forwarded to :func:bfda_power_curve (e.g. alpha0, beta0, sigma_theta, prior_sigma_delta, n_iter, etc.).	{}

Returns:

Type	Description
Figure	The matplotlib Figure.

Source code in bayesprop/utils/utils.py

def plot_bfda_sensitivity(
    theta_A_true: float,
    theta_B_true: float,
    sample_sizes: list[int],
    thresholds: list[float] | None = None,
    n_sim: int = 500,
    seed: int = 42,
    design: str = "nonpaired",
    title: str | None = None,
    ax: plt.Axes | None = None,
    **kwargs: Any,
) -> plt.Figure:
    """Plot BFDA power curves for multiple BF_10 thresholds.

    Works for both paired and non-paired designs.

    Args:
        theta_A_true: Assumed true success rate for model A.
        theta_B_true: Assumed true success rate for model B.
        sample_sizes: List of per-group sample sizes to evaluate.
        thresholds: BF_10 thresholds to compare (default: [3, 6, 10]).
        n_sim: Number of simulated datasets per sample size.
        seed: Random seed for reproducibility.
        design: ``"nonpaired"`` or ``"paired"``.
        title: Optional custom title.
        ax: Optional matplotlib Axes to plot on.
        **kwargs: Extra arguments forwarded to :func:`bfda_power_curve`
            (e.g. ``alpha0``, ``beta0``, ``sigma_theta``,
            ``prior_sigma_delta``, ``n_iter``, etc.).

    Returns:
        The matplotlib Figure.
    """
    if thresholds is None:
        thresholds = [3.0, 6.0, 10.0]

    if ax is None:
        fig, ax = plt.subplots(figsize=(9, 5))
    else:
        fig = ax.get_figure()

    colors = ["#4CAF50", "#FF9800", "#E91E63", "#9C27B0", "#2196F3"]

    for thresh, col in zip(thresholds, colors, strict=False):
        curve = bfda_power_curve(
            theta_A_true=theta_A_true,
            theta_B_true=theta_B_true,
            sample_sizes=sample_sizes,
            design=design,
            decision_rule="bayes_factor",
            bf_threshold=thresh,
            n_sim=n_sim,
            seed=seed,
            **kwargs,
        )

        ax.plot(
            list(curve.keys()),
            list(curve.values()),
            "o-",
            color=col,
            linewidth=2,
            markersize=6,
            label=f"BF_10 > {thresh:.0f}",
        )

    ax.axhline(0.80, color="gray", linestyle="--", linewidth=1, alpha=0.6)
    ax.set_xlabel("Sample size per group (n)")
    ax.set_ylabel("P(BF_10 > threshold)")
    if title is None:
        title = f"BFDA Sensitivity -- {design.title()} Design"
    ax.set_title(title, fontsize=12, fontweight="bold")
    ax.set_ylim(-0.02, 1.02)
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
    fig.tight_layout()

    return fig