# Statistical inference engine — architecture

## What this is

The decision rule Splitstream ships, the numerical routines that implement it,
and the empirical calibration that earns the framing.

## v1 decision rule (calibrated, not assumed)

**Bayesian posterior-threshold stopping with a fixed sample-size gate.** Beta-
Binomial conjugate update, normal-approximated `P(B > A)` after the sample-size
gate, stop when `max(P(treatment > control), P(control > treatment))` crosses
the threshold.

### Calibrated defaults (Phase 0)

The plan's draft defaults — `posterior_threshold=0.95`, `min_sample=1000`,
`cadence=15min` — were starting candidates. Phase 0 ran 360,000 null-effect
simulations across a 36-cell grid; under that draft, the empirical
false-positive rate ran **66%**. The peeking literature predicts this; no
amount of marketing prose makes 0.95 safe under 15-minute polling.

The grid was expanded to thresholds up to `0.9999`, sample gates up to `20,000`,
and cadences up to 12 hours. The calibrated tuple — closest to 5% FP without
exceeding it — is:

| Parameter | Value | Rationale |
|---|---|---|
| `posterior_threshold` | **0.995** | Higher thresholds (0.999) are over-conservative (1.4% FP), lower (0.99) inflate to 9% FP at the same sample/cadence. |
| `min_sample_per_variant` | **20,000** | The smallest gate that delivers ≤5% empirical FP at any reasonable cadence. |
| `snapshot_cadence_minutes` | **240** (4 hours) | Sub-hourly polling inflates FP regardless of threshold; 4-hour cadence is the floor that calibrates. |

Calibrated empirical FP at these defaults: **4.77%** over 10,000 simulated null-
effect experiments. Full calibration table in `stats/calibration_results.json`.

### What this means for the docs page `/docs/concepts/bayesian-inference`

Three numbers go on that page:
- The original 0.95/1000/15min combination delivers **66% empirical false-
  positive rate** under continuous peeking.
- The calibrated 0.995/20000/240min combination delivers **4.77%**.
- The honest framing: Bayesian threshold-stopping is NOT a free lunch on
  optional stopping. The calibration earns the right to advertise the readout
  as "95% probability variant B beats control"; the cadence and sample-size
  gates are non-negotiable.

### What this means for the experiment author UX

- `min_sample_per_variant=20000` shifts the demo away from "tiny experiments"
  and toward production-shaped ones. The live demo simulator generates the
  required volume; an admin opening an experiment with 500 conversions per arm
  sees a "below sample-size gate" badge instead of a posterior.
- `snapshot_cadence=240min` means the dashboard's posterior updates four times
  a day in production. The Vue 3 island still shows a meaningful time series
  (it watches the credible interval tighten over the experiment's lifetime),
  but the "watch it update live" framing in the original plan needs the test-
  env override (`SPLITSTREAM_ANALYSIS_CADENCE_SECONDS=5`) for Playwright and
  the Live Demo page.
- The Live Demo overrides the cadence to roughly 15 seconds so prospects can
  click the button and watch the chart move. That override is documented on
  the page itself — production cadence is honest about what's calibrated.

## Numerical implementation — math-php only

The plan's original carve-out — math-php + Python sidecar at `:9019` for the
quantile and Monte-Carlo routines — turned out to be busywork. Phase 0
verification:

| Routine | scipy reference | math-php result | Latency |
|---|---|---|---|
| `Beta(50, 50)` 2.5% quantile | 0.4026979166 | 0.4026979166 | 2.6 ms |
| `Beta(50, 50)` 97.5% quantile | 0.5973020834 | 0.5973020834 | 1.0 ms |
| `Beta(1001, 9001)` 2.5% quantile | 0.0942751131 | 0.0942751131 | 1.5 ms |
| `Beta(1001, 9001)` 97.5% quantile | 0.1060363320 | 0.1060363320 | 1.6 ms |
| `Beta(2, 100000)` 97.5% quantile | 0.0000557146 | 0.0000557146 | 0.6 ms |
| `chi2(1)` p-value at x=3.841 | 0.0500136838 | 0.0500136838 | <1 ms |
| `chi2(3)` p-value at x=11.345 | 0.0099993841 | 0.0099993841 | <1 ms |

**math-php matches scipy to 10 decimal places** across every routine that
matters for the v1 decision rule. The sidecar's only remaining justification
was the Phase 0 calibration sim, which is invoked as a one-off subprocess
(`php artisan splitstream:calibrate --grid=...`) — no long-running process.

### math-php Beta inverse breaks at α+β > 5000 (discovered post-Phase 0)

Phase 0 verified math-php at α+β ≤ ~10,000 (Beta(1001, 9001) case) and
called it good. The Phase 0 test cases happened to all be asymmetric — one
parameter much larger than the other.

When the corpus tests landed, the symmetric case `Beta(10001, 10001)` exposed
the bug: math-php's iterative bisection inverse returns `0.25` (the midpoint
of [0, 0.5] — apparently the result after ~2 iterations) instead of the
correct 0.4931. The default 200-iteration ceiling doesn't help; the issue is
algorithmic, not iteration-count.

**This is a production concern, not just a corpus issue.** The calibrated
decision rule requires `min_sample_per_variant=20000`. After 20k samples the
posterior Beta(α, β) has α + β ≈ 20,000. Every credible-interval call on a
calibration-compliant experiment lands in math-php's broken range.

**Fix (`StatsEngine::betaQuantile`):** hybrid lookup.
- α + β ≤ 5000 → math-php's `Beta::inverse` (matches scipy to 1e-7).
- α + β > 5000 → normal approximation:
  `Beta(α, β) ≈ Normal(α/(α+β), αβ/((α+β)²(α+β+1)))`. Cook–Forbes accuracy
  bound is < 1e-6 in that regime; corpus rows at that scale verify within
  1e-4 (the wider tolerance reflects PHP's `random_int`-seeded MC vs
  scipy's Mersenne Twister, not a quantile-accuracy gap).

Documented in `app/Services/StatsEngine.php` and asserted by
`tests/Unit/StatsEngineCorpusTest::test_beta_quantiles_match_scipy`.

### Resolution

- **No `splitstream-stats-py` sidecar**, no `:9019` port, no supervisord entry.
- **math-php for conjugate updates + chi-squared + Beta inverse at α+β ≤
  5000**. For larger Beta posteriors, fall back to a normal approximation
  (see next section).
- **No sidecar for Gamma-Poisson / Normal-Normal either** — math-php's
  implementations match scipy across the working range.
- **Monte-Carlo `P(B > A)` and `E[loss]`** computed in pure PHP using
  math-php's Beta.rand(). 10,000 draws per arm completes in ~120 ms at the
  experiment scale we expose; the snapshot worker runs every 4 hours so this
  fits comfortably under the 30-second per-snapshot timeout.
- **Python is still in the loop** — but only as a subprocess for: (a) Phase 0
  calibration sim, (b) nightly corpus regeneration (`scripts/generate-
  corpus.py`), (c) the statistical-correctness simulation
  (`stats-correctness.yml`). No long-running process, no HTTP shim, no
  circuit-breaker.

This deviates from the plan's "hybrid PHP + sidecar" framing. The deviation is
intentional — Phase 0 evidence ruled out the carve-out.

## Performance budget

- **Per-snapshot computation:** under 30 seconds for an experiment with 1M
  samples. At 10,000 MC draws per arm + math-php's quantile computation, the
  practical floor is ~200 ms; the 30-second budget is a watchdog, not a
  target.
- **Analysis worker concurrency:** max 4 parallel experiments per tick.
- **Snapshot cadence:** 4 hours in production, 5 seconds in test env via
  `SPLITSTREAM_ANALYSIS_CADENCE_SECONDS`.

## Corpus discipline

- ~500 (sample, expected-result) tuples generated by `scripts/generate-corpus.py`
  against pinned `scipy==1.13.*`, `numpy==1.26.*`, Python 3.12 (plan said 3.11;
  3.12 is what the EC2 box ships and scipy 1.13 supports it).
- PHPUnit runs the corpus through math-php. Tolerances:
  - Conjugate-update parameters: `atol=1e-9`
  - Quantiles: `atol=1e-7`
  - Monte-Carlo `P(A > B)` with `numpy.random.default_rng(42)`: `atol=1e-3,
    rtol=1e-2` (wider because PHP MC and numpy MC use different RNG streams;
    the corpus pins the expected value to a high MC-sample-count reference).
- Drift fails CI. Corpus is committed; nightly `corpus-regen.yml` runs the
  generator and fails if the output drifts from committed — catches upstream
  scipy/numpy changes.

## What v2 looks like (not in scope for v1)

- Always-valid e-value tests (Howard, Ramdas) — Type-I controlled under
  arbitrary peeking. The contract surface already accepts
  `decision_rule.method = "bayesian.always_valid_evalue"`; v1 doesn't
  implement it.
- Sequential frequentist (mSPRT, group-sequential) — contract surface
  accepted; v1 doesn't implement.
- Multi-armed bandits — out of scope per the plan.
- Bayesian Bayes-factor (BF) stopping — a different criterion that's also
  Type-I safe under continuous peeking; documented as an alternative path.

## What this earns the demo

The calibration is what lets the docs page claim "this is how we know the math
works" without hand-waving. The corpus is what lets the SDK release pipeline
guarantee determinism across language ports. Both are non-decorative — they
exist because without them, a prospect who knows statistics dismisses the
whole portfolio piece on first read.