refactor(stats): extract statistical inference into a Services/Stats service

The command carried ~95 lines of private inference math (Wilson, Fisher,
Newcombe, pooled z, log-factorials), so it changed for two unrelated reasons
- presentation and statistical method - and the math was reachable only through
brittle toContain() assertions on rendered output. Move it to
App\Services\Stats\ProportionSignificance (matching the repo's 'compute in
Services/Stats, present in the command' convention) with a BinomialProportion
value object for the (successes, trials) clump, memoise the log-factorials, and
add unit tests that pin the exact interval and p-value numbers. Behaviour of the
command is unchanged.
This commit is contained in:
Víctor Falcón 2026-07-15 09:04:50 +02:00
parent d763a5b476
commit 013a237242
4 changed files with 265 additions and 136 deletions

View File

@ -4,7 +4,9 @@ namespace App\Console\Commands;
use App\Features\SubscriptionExperiment;
use App\Services\Discord\DiscordWebhook;
use App\Services\Stats\BinomialProportion;
use App\Services\Stats\ExperimentFunnelCollector;
use App\Services\Stats\ProportionSignificance;
use Carbon\CarbonImmutable;
use Illuminate\Console\Command;
@ -22,8 +24,10 @@ class SendExperimentFunnelReportCommand extends Command
SubscriptionExperiment::PAY_NOW => 'pay_now',
];
public function __construct(private ExperimentFunnelCollector $collector)
{
public function __construct(
private ExperimentFunnelCollector $collector,
private ProportionSignificance $significance,
) {
parent::__construct();
}
@ -78,7 +82,7 @@ class SendExperimentFunnelReportCommand extends Command
foreach (self::LABELS as $key => $label) {
$row = $report['variants'][$key];
$mature = $row['assignedMature'] > 0;
$money = $revenue && $mature;
$showMoney = $revenue && $mature;
$lines[] = sprintf(
'%-8s %5d %5d %5d %5d %5d %6s %7s %7s %7s %7s %7s',
@ -89,11 +93,11 @@ class SendExperimentFunnelReportCommand extends Command
$row['assignedMature'],
$row['convertedMature'],
$mature ? ((int) round($row['conversionRate'] * 100)).'%' : 'pend',
$money && $row['arpuCents'] !== null ? $this->money($row['arpuCents'], $currency) : '—',
$money ? $this->money($row['mrrCents'], $currency) : '—',
$showMoney && $row['arpuCents'] !== null ? $this->money($row['arpuCents'], $currency) : '—',
$showMoney ? $this->money($row['mrrCents'], $currency) : '—',
$mature ? $this->money($row['costCents'], $currency) : '—',
$mature ? $this->money($row['wastedCostCents'], $currency) : '—',
$money ? $this->money($row['contributionMarginCents'], $currency) : '—',
$showMoney ? $this->money($row['contributionMarginCents'], $currency) : '—',
);
}
@ -101,19 +105,18 @@ class SendExperimentFunnelReportCommand extends Command
}
/**
* Conversion-rate uncertainty per variant (95% Wilson interval) plus a
* two-proportion z-test between the two leaders, Bonferroni-corrected for
* the three pairwise comparisons, so "check significance before calling a
* winner" has the numbers behind it instead of just the warning.
* Per-variant conversion-rate uncertainty (95% Wilson interval) plus the
* leader-vs-runner-up verdict from {@see ProportionSignificance} a Fisher
* exact test and a Newcombe difference interval, Bonferroni-corrected so
* "check significance before calling a winner" has the numbers behind it.
*
* @param array{startedAt: ?CarbonImmutable, currency: string, revenueAvailable: bool, costPerConnectionCents: int, variants: array<string, array<string, mixed>>} $report
* @return list<string>
*/
private function significanceLines(array $report): array
{
$z = 1.96; // two-sided 95%
$lines = ['', 'Significance (95% Wilson CI on Conv%, n = MatU):'];
$scored = [];
$arms = [];
foreach (self::LABELS as $key => $label) {
$row = $report['variants'][$key];
@ -126,154 +129,43 @@ class SendExperimentFunnelReportCommand extends Command
continue;
}
[$low, $high] = $this->wilsonInterval($k, $n, $z);
[$low, $high] = $this->significance->wilsonInterval($k, $n);
$lines[] = sprintf(' %-8s %6s [%6s %6s] (n=%d)', $label, $this->percent($k / $n), $this->percent($low), $this->percent($high), $n);
$scored[] = ['label' => $label, 'k' => $k, 'n' => $n, 'rate' => $k / $n];
$arms[] = new BinomialProportion($label, $k, $n);
}
if (count($scored) < 2) {
if (count($arms) < 2) {
$lines[] = 'Not enough matured variants to compare yet.';
return $lines;
}
usort($scored, fn (array $a, array $b): int => $b['rate'] <=> $a['rate']);
[$leader, $runnerUp] = [$scored[0], $scored[1]];
// Decide with Fisher's exact test, not a normal-approximation z: at the
// small conversion counts this report has, the pooled z overstates the
// evidence (expected cell counts fall below the np>=5 the z-test needs).
// Fisher is exact at any n. Bonferroni-correct for the 3 pairwise arms.
$alpha = 0.05 / 3;
$delta = ($leader['rate'] - $runnerUp['rate']) * 100;
[$diffLow, $diffHigh] = $this->newcombeDiffInterval($leader['k'], $leader['n'], $runnerUp['k'], $runnerUp['n'], $z);
$fisherP = $this->fisherExactTwoSided(
$leader['k'], $leader['n'] - $leader['k'],
$runnerUp['k'], $runnerUp['n'] - $runnerUp['k'],
);
$significant = $fisherP < $alpha;
usort($arms, fn (BinomialProportion $a, BinomialProportion $b): int => $b->rate() <=> $a->rate());
[$leader, $runnerUp] = [$arms[0], $arms[1]];
$result = $this->significance->compare($leader, $runnerUp);
$lines[] = sprintf(
'Leader %s vs %s: Δ %+.1f pts (95%% CI %+.1f … %+.1f pts, Newcombe).',
$leader['label'], $runnerUp['label'], $delta, $diffLow * 100, $diffHigh * 100,
$leader->label, $runnerUp->label,
($leader->rate() - $runnerUp->rate()) * 100, $result['diffLow'] * 100, $result['diffHigh'] * 100,
);
$lines[] = sprintf(
'Fisher exact p=%.3f %s α=%.3f (Bonferroni×3) -> %s.%s',
$fisherP, $significant ? '<' : '≥', $alpha,
$significant ? 'significant' : 'not significant',
$significant ? '' : ' Keep running.',
$result['fisherP'], $result['significant'] ? '<' : '≥', $result['alpha'],
$result['significant'] ? 'significant' : 'not significant',
$result['significant'] ? '' : ' Keep running.',
);
$pooled = ($leader['k'] + $runnerUp['k']) / ($leader['n'] + $runnerUp['n']);
$minExpected = min($leader['n'], $runnerUp['n']) * min($pooled, 1 - $pooled);
if ($minExpected < 5.0) {
if ($result['minExpectedCount'] < 5.0) {
$lines[] = sprintf(
'(Small sample: min expected conversions %.1f < 5, so the normal-approx z=%.2f overstates — exact test used.)',
$minExpected, $this->twoProportionZ($leader['k'], $leader['n'], $runnerUp['k'], $runnerUp['n']),
$result['minExpectedCount'], $result['z'],
);
}
return $lines;
}
/**
* Wilson score interval for a binomial proportion accurate for small n
* and near 0/1, where the normal approximation misbehaves.
*
* @return array{0: float, 1: float} lower and upper bound, clamped to [0, 1]
*/
private function wilsonInterval(int $k, int $n, float $z): array
{
$p = $k / $n;
$z2 = $z * $z;
$denom = 1 + $z2 / $n;
$center = ($p + $z2 / (2 * $n)) / $denom;
$margin = ($z / $denom) * sqrt($p * (1 - $p) / $n + $z2 / (4 * $n * $n));
return [max(0.0, $center - $margin), min(1.0, $center + $margin)];
}
private function twoProportionZ(int $kA, int $nA, int $kB, int $nB): float
{
$pooled = ($kA + $kB) / ($nA + $nB);
$se = sqrt($pooled * (1 - $pooled) * (1 / $nA + 1 / $nB));
return $se > 0.0 ? (($kA / $nA) - ($kB / $nB)) / $se : 0.0;
}
/**
* Newcombe (Wilson-based) 95% interval for the difference of two proportions
* pA pB. The correct object when asking "is A better than B" overlapping
* marginal intervals do NOT imply the difference includes 0.
*
* @return array{0: float, 1: float}
*/
private function newcombeDiffInterval(int $kA, int $nA, int $kB, int $nB, float $z): array
{
$pA = $kA / $nA;
$pB = $kB / $nB;
[$lA, $uA] = $this->wilsonInterval($kA, $nA, $z);
[$lB, $uB] = $this->wilsonInterval($kB, $nB, $z);
$lower = ($pA - $pB) - sqrt(($pA - $lA) ** 2 + ($uB - $pB) ** 2);
$upper = ($pA - $pB) + sqrt(($uA - $pA) ** 2 + ($pB - $lB) ** 2);
return [$lower, $upper];
}
/**
* Two-sided Fisher exact test p-value for the 2x2 table [[a, b], [c, d]]
* (a/c = conversions, b/d = non-conversions). Sums the hypergeometric
* probabilities of every table, with the same margins, no more likely than
* the observed one. Exact at any sample size no normal approximation.
*/
private function fisherExactTwoSided(int $a, int $b, int $c, int $d): float
{
$rowA = $a + $b;
$rowB = $c + $d;
$col = $a + $c;
$total = $rowA + $rowB;
if ($rowA === 0 || $rowB === 0 || $col === 0 || $col === $total) {
return 1.0;
}
$logProbObserved = $this->hypergeometricLogProb($a, $rowA, $rowB, $col);
$p = 0.0;
for ($x = max(0, $col - $rowB); $x <= min($col, $rowA); $x++) {
$logProb = $this->hypergeometricLogProb($x, $rowA, $rowB, $col);
if ($logProb <= $logProbObserved + 1e-7) {
$p += exp($logProb);
}
}
return min(1.0, $p);
}
private function hypergeometricLogProb(int $x, int $rowA, int $rowB, int $col): float
{
return $this->logChoose($rowA, $x) + $this->logChoose($rowB, $col - $x) - $this->logChoose($rowA + $rowB, $col);
}
private function logChoose(int $n, int $k): float
{
if ($k < 0 || $k > $n) {
return -INF;
}
return $this->logFactorial($n) - $this->logFactorial($k) - $this->logFactorial($n - $k);
}
private function logFactorial(int $n): float
{
$sum = 0.0;
for ($i = 2; $i <= $n; $i++) {
$sum += log($i);
}
return $sum;
}
private function percent(float $rate): string
{
return number_format($rate * 100, 1).'%';

View File

@ -0,0 +1,22 @@
<?php
namespace App\Services\Stats;
/**
* One arm of a binomial experiment: a count of successes out of trials, with a
* label. Bundles the (successes, trials) pair that the significance methods
* would otherwise pass around as loose integers.
*/
final readonly class BinomialProportion
{
public function __construct(
public string $label,
public int $successes,
public int $trials,
) {}
public function rate(): float
{
return $this->trials > 0 ? (float) $this->successes / $this->trials : 0.0;
}
}

View File

@ -0,0 +1,149 @@
<?php
namespace App\Services\Stats;
/**
* Frequentist inference on binomial proportions for the experiment funnel:
* per-arm Wilson intervals, a Newcombe interval for the difference of two
* proportions, and a Fisher exact test for the leader-vs-runner-up comparison.
*
* Fisher (exact at any sample size) is the decision test because the report's
* conversion counts are too small for the two-proportion z normal approximation
* to be valid its expected cell counts fall well below the np>=5 rule.
*/
final class ProportionSignificance
{
/** Two-sided 95% standard-normal quantile. */
private const Z_95 = 1.96;
/** @var list<float> memoised log-factorials, index i = log(i!) */
private array $logFactorials = [0.0, 0.0];
/**
* Compare the two leading arms: Newcombe difference interval, Fisher exact
* p-value, and a family-wise (Bonferroni) corrected significance flag. `z`
* and `minExpectedCount` are returned for the small-sample caveat only.
*
* @return array{alpha: float, diffLow: float, diffHigh: float, fisherP: float, significant: bool, minExpectedCount: float, z: float}
*/
public function compare(BinomialProportion $leader, BinomialProportion $runnerUp, float $familyAlpha = 0.05, int $comparisons = 3): array
{
$alpha = $familyAlpha / $comparisons;
[$diffLow, $diffHigh] = $this->newcombeDiffInterval($leader, $runnerUp);
$fisherP = $this->fisherExactTwoSided(
$leader->successes, $leader->trials - $leader->successes,
$runnerUp->successes, $runnerUp->trials - $runnerUp->successes,
);
$pooled = ($leader->successes + $runnerUp->successes) / ($leader->trials + $runnerUp->trials);
$minExpectedCount = min($leader->trials, $runnerUp->trials) * min($pooled, 1 - $pooled);
return [
'alpha' => $alpha,
'diffLow' => $diffLow,
'diffHigh' => $diffHigh,
'fisherP' => $fisherP,
'significant' => $fisherP < $alpha,
'minExpectedCount' => $minExpectedCount,
'z' => $this->twoProportionZ($leader, $runnerUp),
];
}
/**
* Wilson score interval for a binomial proportion accurate for small n
* and near 0/1, where the normal approximation misbehaves.
*
* @return array{0: float, 1: float} lower and upper bound, clamped to [0, 1]
*/
public function wilsonInterval(int $successes, int $trials, float $z = self::Z_95): array
{
$p = $successes / $trials;
$z2 = $z * $z;
$denom = 1 + $z2 / $trials;
$center = ($p + $z2 / (2 * $trials)) / $denom;
$margin = ($z / $denom) * sqrt($p * (1 - $p) / $trials + $z2 / (4 * $trials * $trials));
return [max(0.0, $center - $margin), min(1.0, $center + $margin)];
}
/**
* Newcombe (Wilson-based) 95% interval for the difference pA pB. The
* correct object for "is A better than B": overlapping marginal intervals do
* NOT imply the difference includes 0.
*
* @return array{0: float, 1: float}
*/
public function newcombeDiffInterval(BinomialProportion $a, BinomialProportion $b, float $z = self::Z_95): array
{
$pA = $a->rate();
$pB = $b->rate();
[$lA, $uA] = $this->wilsonInterval($a->successes, $a->trials, $z);
[$lB, $uB] = $this->wilsonInterval($b->successes, $b->trials, $z);
$lower = ($pA - $pB) - sqrt(($pA - $lA) ** 2 + ($uB - $pB) ** 2);
$upper = ($pA - $pB) + sqrt(($uA - $pA) ** 2 + ($pB - $lB) ** 2);
return [$lower, $upper];
}
/** Pooled two-proportion z statistic — descriptive only, not the decision test. */
public function twoProportionZ(BinomialProportion $a, BinomialProportion $b): float
{
$pooled = ($a->successes + $b->successes) / ($a->trials + $b->trials);
$se = sqrt($pooled * (1 - $pooled) * (1 / $a->trials + 1 / $b->trials));
return $se > 0.0 ? ($a->rate() - $b->rate()) / $se : 0.0;
}
/**
* Two-sided Fisher exact p-value for the 2x2 table [[a, b], [c, d]]
* (a/c = successes, b/d = failures). Sums the hypergeometric probabilities
* of every same-margin table no more likely than the observed one. Exact at
* any sample size no normal approximation.
*/
public function fisherExactTwoSided(int $a, int $b, int $c, int $d): float
{
$rowA = $a + $b;
$rowB = $c + $d;
$col = $a + $c;
$total = $rowA + $rowB;
if ($rowA === 0 || $rowB === 0 || $col === 0 || $col === $total) {
return 1.0;
}
$logProbObserved = $this->hypergeometricLogProb($a, $rowA, $rowB, $col);
$p = 0.0;
for ($x = max(0, $col - $rowB); $x <= min($col, $rowA); $x++) {
$logProb = $this->hypergeometricLogProb($x, $rowA, $rowB, $col);
if ($logProb <= $logProbObserved + 1e-7) {
$p += exp($logProb);
}
}
return min(1.0, $p);
}
private function hypergeometricLogProb(int $x, int $rowA, int $rowB, int $col): float
{
return $this->logChoose($rowA, $x) + $this->logChoose($rowB, $col - $x) - $this->logChoose($rowA + $rowB, $col);
}
private function logChoose(int $n, int $k): float
{
if ($k < 0 || $k > $n) {
return -INF;
}
return $this->logFactorial($n) - $this->logFactorial($k) - $this->logFactorial($n - $k);
}
private function logFactorial(int $n): float
{
for ($i = count($this->logFactorials); $i <= $n; $i++) {
$this->logFactorials[$i] = $this->logFactorials[$i - 1] + log($i);
}
return $this->logFactorials[$n];
}
}

View File

@ -0,0 +1,66 @@
<?php
use App\Services\Stats\BinomialProportion;
use App\Services\Stats\ProportionSignificance;
beforeEach(function () {
$this->stats = new ProportionSignificance;
});
it('computes the Wilson score interval', function () {
[$low, $high] = $this->stats->wilsonInterval(6, 52);
expect($low)->toEqualWithDelta(0.0540, 0.0005)
->and($high)->toEqualWithDelta(0.2297, 0.0005);
[$low2, $high2] = $this->stats->wilsonInterval(5, 179);
expect($low2)->toEqualWithDelta(0.0120, 0.0005)
->and($high2)->toEqualWithDelta(0.0640, 0.0005);
});
it('keeps the Wilson interval inside [0, 1] at the k=0 and k=n boundaries', function () {
expect($this->stats->wilsonInterval(0, 30)[0])->toBe(0.0)
->and($this->stats->wilsonInterval(30, 30)[1])->toBe(1.0);
});
it('computes a two-sided Fisher exact p-value', function () {
// Real experiment table: reduced 6/52 vs pay_now 5/179.
expect($this->stats->fisherExactTwoSided(6, 46, 5, 174))->toEqualWithDelta(0.0182, 0.0005);
// Clean separation 5/5 vs 0/5 → 2 * C(5,5)/C(10,5) = 2/252.
expect($this->stats->fisherExactTwoSided(5, 0, 0, 5))->toEqualWithDelta(0.007936, 0.00001);
// Degenerate margin (no successes anywhere) → p = 1.
expect($this->stats->fisherExactTwoSided(0, 10, 0, 10))->toBe(1.0);
});
it('computes the Newcombe difference interval and the descriptive z', function () {
$reduced = new BinomialProportion('reduced', 6, 52);
$payNow = new BinomialProportion('pay_now', 5, 179);
[$low, $high] = $this->stats->newcombeDiffInterval($reduced, $payNow);
expect($low)->toEqualWithDelta(0.0164, 0.0005)
->and($high)->toEqualWithDelta(0.2029, 0.0005)
->and($this->stats->twoProportionZ($reduced, $payNow))->toEqualWithDelta(2.607, 0.01);
});
it('calls the borderline real comparison NOT significant under Bonferroni', function () {
$result = $this->stats->compare(
new BinomialProportion('reduced', 6, 52),
new BinomialProportion('pay_now', 5, 179),
);
// Fisher p ≈ 0.018 exceeds the corrected bar 0.05/3 ≈ 0.0167 → not a winner yet.
expect($result['significant'])->toBeFalse()
->and($result['fisherP'])->toEqualWithDelta(0.0182, 0.0005)
->and($result['alpha'])->toEqualWithDelta(0.0167, 0.0005)
->and($result['minExpectedCount'])->toEqualWithDelta(2.48, 0.05);
});
it('calls a clean separation significant', function () {
$result = $this->stats->compare(
new BinomialProportion('a', 5, 5),
new BinomialProportion('b', 0, 5),
);
// Fisher p ≈ 0.008 clears the corrected bar.
expect($result['significant'])->toBeTrue()
->and($result['fisherP'])->toEqualWithDelta(0.007936, 0.00001);
});