AB Testing is an Optimization Problem

Significance Testing is costing you revenue.

AB Testing Nightmare: 100 Close Calls

Your organization will run 100 A/B tests next year. In every single one, treatment outperforms control in the sample. But the p-value never quite crosses 0.05 — it lands at 0.06, 0.07, 0.08. Close, but not there.

Your vaunted stats expert calls it every time: “Not significant. Keep control.”

Do you abide? Or will you push back?

100 tests. 100 times, the data points toward the treatment. 100 times, you shipped nothing. The cost of that policy — across a year, across a product — how much revenue have you sacrificed in the name of statistical purity?

Why Are We AB Testing?

The purpose of an A/B test, in industry, is straightforward: find the variant that performs better and ship it. More clicks, more conversions, more revenue.

Meanwhile, Significance Testing was designed with science in mind — a field where false positives are costly and reproducibility depends on tight error control. α caps the false positive rate; β limits false negatives.

But neither tells you whether acting on the result is worth it, or the potential magnitude of an error.

Mindset Shift: Value-Focused Testing

This paper proposes a mindset shift. Rather than obsess over statistical significance in the search for low p-values, design experiments to maximize value.

Instead of asking “is this result statistically significant?”, ask “at what point does continuing this test cost more than it returns?”

The threshold that answers that question is the Optimal Stopping Rule.

Computing the Optimal Stopping Rule

Instead of fixing α, β, and a single — often unrealistic — Minimum Detectable Effect, we set three inputs:

• Max Test N — the absolute maximum number of users allowed into the test
• Post-Test N — the audience that sees the winner after the test ends
• Treatment Effects — drawn from a realistic distribution, not a single fixed point

With those inputs fixed, we simulate 100,000 experiments — each drawing a different treatment effect, running the test, and evaluating every possible stopping threshold. We pick the one that maximizes average long-run value across all simulations.

Why Use Simulations?

Simulations give us something no real experiment can: a ground truth to judge against. In practice, we never know the true effect — except in our simulations, where we encode it. Each draw is different, modeling the full range of experiments a team might actually run: large effects, small effects, null effects, and genuinely positive treatments that look like losers in the sample by chance. We measure which stopping rule performs best across all of it — without knowing the true effect in advance, just as in a real experiment.

Define Experiment Parameters

These parameters define your experiment — baseline conversion, expected treatment effects, test traffic, and post-test audience — and directly determine the optimal stopping threshold. Click any to learn more.

Conversion

Traffic

Split

Min / Max Days

Post-Test Traffic

The pre-experiment conversion rate. Sets the variance of your test statistic and determines how much room for improvement exists.

Test Traffic & Split

Total users in the experiment and their control/treatment split. The first 7 days are burn-in — no stopping decisions. The experiment runs at most 20 days.

Post-Test Traffic

The population exposed to the winner after the decision is made. More post-test traffic means more value at stake — potentially justifying a more aggressive stopping rule.

The realistic lift you expect comes from a distribution of possible effects — not a single number. Here, experiments have an equal chance of being good or bad.

Simulation: 1

Day: 1

Simulation: 1

Simulation: 1

How the Simulation Works

The same simulated data is evaluated under 50 different stopping rules — Pr(Best) thresholds from 0.50 to 0.99. For each:

1. Burn-in: No stopping before day min_test_days = 7
2. Daily check: If either variant’s Pr(Best) crosses the threshold, the experiment stops.
3. Rollout: The winning variant is deployed to all remaining test users and the post-test audience.
4. Score: Record the long-run value produced.

The threshold with the highest average score across all simulations is the Optimal Stopping Rule.

Metric	Description
Projected Mean	Average conversion rate across all users: both arms during the test, plus the winner deployed to the full post-test audience.
Marginal Events vs No Test	Incremental conversions above the baseline of never running an experiment.
Correct Decision Rate	How often the strategy deploys the correct variant. Covers the full distribution of effects — beneficial and harmful — not just the single MDE statistical power is designed for.
Avg Trials (Days)	Total user-experiment observations before stopping, with calendar duration in parentheses. Shorter experiments compound value by freeing up traffic faster.

Optimization Curve

We swept every Pr(Best) threshold from 0.50 to 0.99. For each, we applied that stopping rule to all 100,000 experiments — each drawing a different true effect, known only to us — and recorded the average long-run value produced across test traffic and the full post-test rollout. The curve plots the average across all simulations. The peak is the threshold that maximizes expected value.

Change the scenario — notice whether the peaks shift or hold.

Baseline Conversion 60%

Post-Test N 50,000 100,000 150,000 200,000 300,000 400,000 600,000 800,000 1,200,000 1,600,000

Expected Effects i ITC ~ N(100, 1)

Given our experiment design, the curve peaks at Pr(Best) = 0.81 — the threshold that balances speed and precision. Stop sooner and you trade accuracy for time. Wait longer and the extra certainty doesn't pay for itself.

Post-Test N is the audience that sees the winning variant after the experiment ends — the size of the real-world rollout. It matters because the value of an experiment isn't just what happens during the test: it's the cumulative impact of deploying the winner to every user who follows. A wrong call during a test with a small rollout is a minor cost; the same wrong call with a million-person rollout is a much larger one. As Post-Test N grows, the stakes of each stopping decision rise — and the optimal threshold rises with it, demanding stronger evidence before committing.

Change the scenario — notice whether the peaks shift or hold.

Baseline Conversion 60%

Expected Effects ITC ~ N(100, 1)

Metric Optimal Pr(Best) Avg Days Marginal Events

Optimal Stopping vs Significance: A Comparison

The same 100,000 simulated experiments — evaluated under Optimal Stopping and Fixed Sample Significance Testing. Neither strategy had access to the true effect we encoded. They only saw what any real analyst sees: the data as it accumulated.

• Optimal Stopping returns 2.0× as many marginal events as Significance Testing (α = 0.05) and takes half the time.
• Against Fixed Sample (α = 0.50), Optimal Stopping finishes 10 days sooner, producing 114 more marginal events.
  • Both strategies ship the same winner — so any difference in value is attributable entirely to early stopping.

Why It Works

The Optimal Stopping algorithm adapts to the signal it sees:

• Small effects — the value gap between variants is narrow, so the cost of a wrong call is low. The algorithm rolls out the winner quickly, sacrificing some accuracy for speed.
• Large positive effects — a clear winning condition emerges early. The algorithm stops, exploits the finding, and captures value before a Fixed Sample test would even reach significance.
• Large negative effects — harm is detectable fast. The algorithm cuts traffic to the losing variant before it does meaningful damage.

Conclusion

Stopping rules are arbitrary without context. Fixed Sample Significance Testing borrows a threshold designed for science — not for maximizing revenue. It runs longer than necessary, defaults to control when the data already points elsewhere, and treats a threshold with no connection to long-run value as the decision criterion.

The simulations show the cost. Compared to Fixed Sample Significance Testing (α = 0.05), Optimal Stopping ships the right variant more often, in fewer days, with more conversions over the full rollout window. The right stopping rule is derived from your experiment design — baseline conversion, audience size, effect distribution, and rollout size. It maximizes expected value, not statistical certainty.

Extras: The Optimization Toolkit

⚡ 100,000 simulated experiments × 20 days each × two recipes — 2.5s

To run your own optimization, swap in your experiment parameters — baseline conversion, Max Test N, Post-Test N, effect distribution — and the pipeline handles the rest.

The whole thing is fully vectorized. Every simulation, every day, every threshold is computed as a single set of array operations — no Python for loops anywhere. Under the hood, simulationengine uses Polars for columnar data transformations and NumPy for the numerical sweep.

Four steps: simulate 100\,000 experiments, compute daily running stats, score each potential stopping decision, then sweep every Pr(Best) threshold to find the one that maximizes long-run value.

import simulationengine as se
import numpy as np
import polars as pl

# ── Parameters ────────────────────────────────────────────────────────
CONTROL_MEAN       = 0.60      # baseline conversion rate
MAX_TEST_N         = 100_000   # total users in the experiment
SPLIT              = 0.50      # 50/50 control/treatment split
MIN_TEST_DAYS      = 7         # burn-in: no stopping before day 7
NUMBER_TEST_DAYS   = 20        # max experiment duration
POST_TEST_N        = 100_000   # post-test audience size
NUMBER_SIMULATIONS = 100_000
SEED               = 42
ITC_MEAN           = 100       # index-to-control mean (100 = no lift)
ITC_SE             = 1         # spread of treatment effects

# ── Step 1: Simulate experiments ──────────────────────────────────────
df = se.create_simulation_df(
    number_test_days=NUMBER_TEST_DAYS,
    number_simulations=NUMBER_SIMULATIONS,
    seed=SEED,
    condition_parameters=[{"mean": CONTROL_MEAN, "trials": int(MAX_TEST_N * SPLIT)}],
)
df = se.add_random_effect_condition(df, {
    "itc_mean": ITC_MEAN,
    "itc_se": ITC_SE,
    "trials": int(MAX_TEST_N * SPLIT),
})

# ── Step 2: Compute running stats ─────────────────────────────────────
df = se.add_cumulative_columns(df)
df = se.compute_cumulative_std_rates(df)
df = se.compute_delta_stats(df)
df = se.compute_test_stats(df)

# ── Step 3: Score each stopping decision ──────────────────────────────
df = se.calculate_remaining_test_potential(df)
df = se.calculate_launch_potential(df, post_test_n=POST_TEST_N)
df = se.calculate_total_potential(df)

# ── Step 4: Find the optimal Pr(Best) threshold ───────────────────────
p_thresholds = np.arange(0.5, 1.00, 0.01)

df_prepped = se.prepare_stopping_decision_columns(
    df.filter(pl.col("day") >= MIN_TEST_DAYS),
    metric_direction_improvement="increase",
).sort(["simulation", "day"])

results = se.evaluate_probability_thresholds(
    df_prepped,
    threshold_values=p_thresholds,
    metric_direction_improvement="increase",
    pre_sorted=True,
    select_best_at_end=True,
)

optimal_threshold = p_thresholds[results["mean_potential_value"].arg_max()]
print(f"Optimal Pr(Best) threshold: {optimal_threshold:.2f}")

create_simulation_df generates a DataFrame of daily binomial draws for one or more fixed-rate conditions. Each row is one day of one simulation.

import numpy as np
import polars as pl


def create_simulation_df(
    number_test_days,
    number_simulations,
    condition_parameters,
    seed=None,
):
    if seed is not None:
        np.random.seed(seed)

    # trials_per_day is the same for every day in the experiment
    n = np.array(
        [c["trials"] // number_test_days for c in condition_parameters]
    )[:, None]
    p = np.array([c["mean"] for c in condition_parameters])[:, None]

    rep = number_test_days * number_simulations
    n_mat = np.repeat(n, rep, axis=1)
    p_mat = np.repeat(p, rep, axis=1)

    # One binomial draw per (condition, day, simulation)
    samples = np.random.binomial(n=n_mat, p=p_mat)

    return pl.DataFrame(
        {
            "simulation": np.arange(rep) // number_test_days + 1,
            "day": np.arange(rep) % number_test_days + 1,
            **{
                f"condition_{i}_value": samples[i]
                for i in range(len(condition_parameters))
            },
            **{
                f"condition_{i}_trials": n_mat[i]
                for i in range(len(condition_parameters))
            },
            **{
                f"condition_{i}_true_mean": np.full(
                    rep, condition_parameters[i]["mean"]
                )
                for i in range(len(condition_parameters))
            },
        }
    )

add_random_effect_condition appends a treatment arm where the true effect varies across simulations but is fixed within one. The effect is drawn from an ITC (Index-to-Control) distribution: ITC = 100 → no lift; ITC ~ N(100, 1) → effects centered at parity with modest spread.

import numpy as np
import polars as pl


def add_random_effect_condition(df, random_effect_parameters, seed=None):
    if seed is not None:
        np.random.seed(seed)

    n_simulations = int(df["simulation"].max())
    n_days = int(df["day"].max())

    # One ITC draw per simulation — the unknown true effect for that experiment
    itc_samples = (
        np.random.normal(
            random_effect_parameters["itc_mean"],
            random_effect_parameters["itc_se"],
            n_simulations,
        )
        / 100
    )

    # treatment_mean = (itc / 100) x control_mean
    control_mean = float(df["condition_0_true_mean"][0])
    treatment_means = itc_samples * control_mean

    # Daily binomial draws; treatment_mean is constant across days within a simulation
    trials_per_day = random_effect_parameters["trials"] // n_days
    events = np.random.binomial(
        n=trials_per_day,
        p=treatment_means[:, None],
        size=(n_simulations, n_days),
    )

    # Reshape to tall format to match the existing DataFrame
    next_condition = (
        max(
            int(c.split("_")[1])
            for c in df.columns
            if c.startswith("condition_") and c.split("_")[1].isdigit()
        )
        + 1
    )
    return df.with_columns(
        [
            pl.Series(
                f"condition_{next_condition}_value",
                events.flatten(),
            ),
            pl.Series(
                f"condition_{next_condition}_trials",
                np.full(n_simulations * n_days, trials_per_day),
            ),
            pl.Series(
                f"condition_{next_condition}_true_mean",
                np.repeat(treatment_means, n_days),
            ),
        ]
    )

add_cumulative_columns converts daily event counts into running totals grouped by simulation. Each condition gets a cumulative event count, trial count, and conversion rate. Cross-condition totals and remaining trials are derived in a single vectorized pass — no loops.

import polars as pl


def add_cumulative_columns(df):
    value_cols  = [col for col in df.columns if col.endswith("_value")]
    trials_cols = [col for col in df.columns if col.endswith("_trials")]

    # Running totals per condition, grouped by simulation
    cumsum_exprs = [
        pl.col(col).cum_sum().over("simulation").alias(f"{col}_cumsum")
        for col in value_cols + trials_cols
    ]
    df = df.with_columns(cumsum_exprs)

    # Cross-condition totals
    value_cumsums  = [f"{c}_cumsum" for c in value_cols]
    trials_cumsums = [f"{c}_cumsum" for c in trials_cols]
    df = df.with_columns([
        pl.sum_horizontal(value_cumsums).alias("total_value_cumsum"),
        pl.sum_horizontal(trials_cumsums).alias("total_trials_cumsum"),
    ])

    # Trials remaining in the experiment window
    df = df.with_columns([
        (pl.col("total_trials_cumsum").max().over("simulation")
         - pl.col("total_trials_cumsum")).alias("total_remaining_test_trials")
    ])

    # Cumulative conversion rate per condition
    mean_exprs = [
        (pl.col(f"{c}_cumsum") / pl.col(f"{c.replace('_value','_trials')}_cumsum"))
        .alias(f"{c.replace('_value','')}_mean_cumsum")
        for c in value_cols
        if f"{c.replace('_value','_trials')}" in trials_cols
    ]
    return df.with_columns(mean_exprs)

compute_cumulative_std_rates applies the binomial standard deviation formula sqrt(p * (1 - p)) to each condition's running conversion rate. Note: this is the std of the proportion, not the standard error — divide by sqrt(n) to get SE.

import polars as pl


def compute_cumulative_std_rates(df):
    mean_cumsum_cols = [col for col in df.columns if col.endswith("_mean_cumsum")]

    std_exprs = [
        (pl.col(col) * (1 - pl.col(col)))
        .sqrt()
        .alias(col.replace("_mean_cumsum", "_std_cumsum"))
        for col in mean_cumsum_cols
    ]
    return df.with_columns(std_exprs)

compute_delta_stats computes two things for each treatment condition vs control: the standard error (SE = std / sqrt(n)) and the delta statistics (delta_mean = mean_i - mean_0, delta_se = sqrt(SE_i² + SE_0²)). These are the direct inputs to z-scores and p-values — this function does not compute those itself.

import polars as pl


def compute_delta_stats(df):
    mean_cols = [c for c in df.columns if c.endswith("_mean_cumsum")]

    # SE = std / sqrt(trials) for every condition
    se_exprs = [
        (pl.col(f"condition_{i}_std_cumsum") / pl.col(f"condition_{i}_trials_cumsum").sqrt())
        .alias(f"condition_{i}_se_cumsum")
        for i in [int(c.split("_")[1]) for c in mean_cols if c.startswith("condition_")]
        if f"condition_{i}_std_cumsum" in df.columns
    ]
    df = df.with_columns(se_exprs)

    # Delta mean and delta SE for each treatment vs condition_0 (control)
    treatment_nums = [
        int(c.split("_")[1]) for c in mean_cols
        if c.startswith("condition_") and int(c.split("_")[1]) > 0
    ]
    delta_exprs = []
    for i in treatment_nums:
        delta_exprs += [
            (pl.col(f"condition_{i}_mean_cumsum") - pl.col("condition_0_mean_cumsum"))
            .alias(f"condition_{i}_delta_mean"),
            (pl.col(f"condition_{i}_se_cumsum").pow(2) + pl.col("condition_0_se_cumsum").pow(2))
            .sqrt()
            .alias(f"condition_{i}_delta_se"),
        ]
    return df.with_columns(delta_exprs)

compute_test_stats converts delta statistics into probabilities and expected values using the normal CDF. For each treatment: Pr(treatment > control) = Φ(delta_mean / delta_se). Expected gain and loss use truncated normal formulas. These are the columns the stopping rule reads — it stops when Pr(Best) crosses the threshold α.

import numpy as np
import polars as pl
from scipy.stats import norm


def compute_test_stats(df):
    delta_mean_cols = [c for c in df.columns if c.endswith("_delta_mean")]
    condition_nums  = [int(c.split("_")[1]) for c in delta_mean_cols]

    pr_cols = []
    for i in condition_nums:
        z = df[f"condition_{i}_delta_mean"].to_numpy() / df[f"condition_{i}_delta_se"].to_numpy()
        pr = norm.cdf(z)
        df = df.with_columns([
            pl.lit(pr).alias(f"condition_{i}_pr_positive_diff"),
            pl.lit(1 - pr).alias(f"condition_{i}_pr_lowest_diff"),
        ])
        pr_cols.append(f"condition_{i}_pr_positive_diff")

        # Expected gain and loss (truncated normal)
        pdf_z = norm.pdf(z)
        df = df.with_columns([
            pl.lit(np.abs(pr   * (df[f"condition_{i}_delta_mean"].to_numpy() + df[f"condition_{i}_delta_se"].to_numpy() * pdf_z / np.maximum(pr,     1e-10)))).alias(f"condition_{i}_expected_gain"),
            pl.lit(np.abs((1-pr) * (-df[f"condition_{i}_delta_mean"].to_numpy() + df[f"condition_{i}_delta_se"].to_numpy() * pdf_z / np.maximum(1-pr, 1e-10)))).alias(f"condition_{i}_expected_loss"),
        ])

    # Control: Pr(best) = 1 - max(Pr(treatment > control))
    if pr_cols:
        df = df.with_columns([
            (1 - pl.max_horizontal(pr_cols)).alias("condition_0_pr_highest_diff"),
            pl.min_horizontal(pr_cols).alias("condition_0_pr_lowest_diff"),
        ])
    return df

For each condition, draws from Binomial(total_remaining_test_trials, true_mean) — one random realization of the value you'd capture by stopping now and deploying that condition to remaining test users.

import numpy as np
import polars as pl


def calculate_remaining_test_potential(df, seed=None):
    if not isinstance(df, pl.DataFrame):
        raise TypeError(f"df must be a polars DataFrame, got {type(df)}")
    if "total_remaining_test_trials" not in df.columns:
        raise ValueError("df must contain 'total_remaining_test_trials' column")
    if seed is not None:
        np.random.seed(seed)

    true_mean_cols = [col for col in df.columns if col.endswith("_true_mean")]
    if not true_mean_cols:
        raise ValueError("df must contain at least one column ending with '_true_mean'")

    potential_exprs = []
    for mean_col in true_mean_cols:
        condition_prefix = mean_col.replace("_true_mean", "")
        n = df["total_remaining_test_trials"].to_numpy().astype(int)
        p = np.clip(df[mean_col].to_numpy(), 0.0, 1.0)
        draws = np.random.binomial(n=n, p=p)
        potential_exprs.append(
            pl.lit(draws).alias(f"{condition_prefix}_test_potential_value")
        )

    return df.with_columns(potential_exprs)

Draws from Binomial(post_test_n, true_mean) for each condition — one random realization of the value you'd capture by launching that condition to the post-test population. This is where stopping decisions have the most leverage.

import numpy as np
import polars as pl


def calculate_launch_potential(df, post_test_n, seed=None):
    if not isinstance(df, pl.DataFrame):
        raise TypeError(f"df must be a polars DataFrame, got {type(df)}")
    if not isinstance(post_test_n, (int, np.integer)):
        raise TypeError(f"post_test_n must be an integer, got {type(post_test_n)}")
    if post_test_n < 0:
        raise ValueError(f"post_test_n must be non-negative, got {post_test_n}")
    if seed is not None:
        np.random.seed(seed)

    true_mean_cols = [col for col in df.columns if col.endswith("_true_mean")]
    if not true_mean_cols:
        raise ValueError("df must contain at least one column ending with '_true_mean'")

    potential_exprs = []
    for mean_col in true_mean_cols:
        condition_prefix = mean_col.replace("_true_mean", "")
        p = np.clip(df[mean_col].to_numpy(), 0.0, 1.0)
        draws = np.random.binomial(n=post_test_n, p=p)
        potential_exprs.append(
            pl.lit(draws).alias(f"{condition_prefix}_launch_potential_value")
        )

    return df.with_columns(potential_exprs)

total_potential = total_value_cumsum + test_potential + launch_potential. The cumulative observed value is shared across all conditions; only the forward-looking components differ. The condition with the highest total potential is the one to deploy.

import polars as pl


def calculate_total_potential(df):
    if not isinstance(df, pl.DataFrame):
        raise TypeError(f"df must be a polars DataFrame, got {type(df)}")
    if "total_value_cumsum" not in df.columns:
        raise ValueError("df must contain 'total_value_cumsum' column")

    test_potential_cols = [
        col for col in df.columns if col.endswith("_test_potential_value")
    ]
    if not test_potential_cols:
        raise ValueError("df must contain at least one column ending with '_test_potential_value'")

    total_exprs = []
    for test_col in test_potential_cols:
        condition_prefix = test_col.replace("_test_potential_value", "")
        launch_col = f"{condition_prefix}_launch_potential_value"
        if launch_col not in df.columns:
            raise ValueError(f"Missing required column: {launch_col}")
        total_exprs.append(
            (pl.col("total_value_cumsum") + pl.col(test_col) + pl.col(launch_col))
            .alias(f"{condition_prefix}_total_potential_value")
        )

    return df.with_columns(total_exprs)

The core of the optimization loop. Rather than iterating over thresholds in Python, it stacks all Pr(Best) columns into a 3D NumPy array (n_sims × n_days × n_conditions), finds the first crossing day per simulation per threshold in a single broadcast, then aggregates value, trial count, and decision accuracy across all thresholds simultaneously.

import numpy as np
import polars as pl


def evaluate_probability_thresholds(
    df_prepped,
    threshold_values,
    metric_direction_improvement="increase",
    pre_sorted=False,
    select_best_at_end=False,
):
    if not isinstance(df_prepped, pl.DataFrame):
        raise TypeError(f"df_prepped must be a polars DataFrame, got {type(df_prepped)}")
    if metric_direction_improvement not in ["increase", "decrease"]:
        raise ValueError(f"metric_direction_improvement must be 'increase' or 'decrease'")
    if not pre_sorted:
        df_prepped = df_prepped.sort(["simulation", "day"])

    n_sims = df_prepped["simulation"].n_unique()
    n_days = df_prepped.height // n_sims

    def _col2d(name):
        return df_prepped[name].to_numpy().reshape(n_sims, n_days)

    # Select probability columns based on metric direction
    if metric_direction_improvement == "increase":
        treatment_prob_cols = [
            c for c in df_prepped.columns
            if c.endswith("_pr_positive_diff") and not c.startswith("condition_0_")
        ]
        control_prob_col = "condition_0_pr_highest_diff"
    else:
        treatment_prob_cols = [
            c for c in df_prepped.columns
            if c.endswith("_pr_lowest_diff") and not c.startswith("condition_0_")
        ]
        control_prob_col = "condition_0_pr_lowest_diff"

    all_prob_cols = [control_prob_col] + treatment_prob_cols
    condition_nums = np.array([int(c.split("_")[1]) for c in all_prob_cols])

    # Stack: (n_sims, n_days, n_conditions)
    prob_stack = np.stack([_col2d(c) for c in all_prob_cols], axis=2)
    max_prob   = prob_stack.max(axis=2)            # (n_sims, n_days)
    best_cond  = condition_nums[prob_stack.argmax(axis=2)]

    # Cumulative max — first day cummax >= threshold is the stopping day
    cummax = np.maximum.accumulate(max_prob, axis=1)
    thresholds_arr = np.array(sorted(threshold_values))

    # First-crossing index, clamped to last day for never-crossed sims
    raw_stop   = (cummax[:, :, np.newaxis] < thresholds_arr).sum(axis=1)  # (n_sims, n_thresholds)
    any_crossed = raw_stop < n_days
    stop_idx   = np.minimum(raw_stop, n_days - 1)

    row = np.arange(n_sims)[:, np.newaxis]
    best_cond_at_stop = best_cond[row, stop_idx]
    winning_cond = np.where(any_crossed | select_best_at_end, best_cond_at_stop, 0)

    # Winning condition's total potential value at stopping day
    pot_cols = sorted(
        [c for c in df_prepped.columns
         if c.endswith("_total_potential_value") and c.startswith("condition_") and c.split("_")[1].isdigit()],
        key=lambda c: int(c.split("_")[1]),
    )
    pot_stack      = np.stack([_col2d(c) for c in pot_cols], axis=2)
    win_pot_at_stop = pot_stack[row, stop_idx, winning_cond]

    total_at_stop    = _col2d("total_trials_cumsum")[row, stop_idx]
    true_cond_at_stop = _col2d("true_best_condition")[row, stop_idx]

    return pl.DataFrame({
        "threshold":           thresholds_arr,
        "mean_potential_value": win_pot_at_stop.mean(axis=0),
        "mean_trials":          total_at_stop.mean(axis=0),
        "median_trials":        np.median(total_at_stop, axis=0),
        "p80_trials":           np.percentile(total_at_stop, 80, axis=0),
        "match_rate":           (winning_cond == true_cond_at_stop).mean(axis=0),
    })

Acknowledgements

1. The idea of designing experiments for value rather than significance is inspired by Test & Roll: Profit-Maximizing A/B Tests by Ron Berman and Elea McDonnell Feit. Their core insight: for tactical experiments, the cost of a wrong call on a small effect is negligible — the cost of oversizing the test to satisfy α is not. This paper extends that spirit by making the stopping rule, not the sample size, the optimization target.

2. My colleague Joseph Powers, PhD introduced me to the idea that A/B testing is an optimization problem and fundamentally changed how I approach experimentation. He showed me the power of simulation studies for building real intuition — not just running tests. A generous mentor and an invaluable teammate.