Statistical Inference Playbook

1/**2 * Example 48: Statistical Inference Playbook3 *4 * Covers the v1.0.0 inference layer that was not represented in the earlier5 * stats examples: confidence intervals, bootstrap uncertainty, Gaussian KDE,6 * multiple-comparison correction, and power analysis.7 */89import { mkdir } from "node:fs/promises";10import { tensor } from "deepbox/ndarray";11import { axhline, figure, groupedBar, kdeplot, legend, saveFig } from "deepbox/plot";12import {13  benjaminiHochberg,14  bonferroni,15  bootstrap,16  cohenD,17  gaussian_kde,18  meanConfidenceInterval,19  meanConfidenceIntervalZ,20  meanDiffConfidenceInterval,21  proportionConfidenceInterval,22  tTestPower,23  ttest_ind,24} from "deepbox/stats";2526const OUTPUT_DIR = "docs/examples/48-statistical-inference-playbook/output";2728console.log("=".repeat(72));29console.log("Example 48: Statistical Inference Playbook");30console.log("=".repeat(72));3132await mkdir(OUTPUT_DIR, { recursive: true });3334// A compact A/B rollout dataset for a checkout experiment.35const controlRevenuePerSession = [86, 91, 88, 94, 90, 96, 84, 89, 92, 87, 95, 90];36const treatmentRevenuePerSession = [94, 101, 99, 104, 100, 107, 92, 98, 102, 97, 105, 100];3738const controlLatencyMs = [1210, 1185, 1224, 1178, 1202, 1194, 1218, 1189];39const treatmentLatencyMs = [1124, 1098, 1116, 1108, 1121, 1095, 1113, 1104];4041const controlCsat = [4.1, 4.2, 4.0, 4.3, 4.1, 4.2, 4.0, 4.1];42const treatmentCsat = [4.3, 4.4, 4.2, 4.5, 4.4, 4.3, 4.2, 4.4];4344const controlResolutionRate = [0.74, 0.71, 0.76, 0.73, 0.75, 0.72, 0.74, 0.73];45const treatmentResolutionRate = [0.81, 0.79, 0.83, 0.82, 0.8, 0.81, 0.84, 0.8];4647const controlConversions = { successes: 158, total: 200 };48const treatmentConversions = { successes: 186, total: 205 };4950// ============================================================================51// Part 1: Mean confidence intervals and uplift intervals52// ============================================================================53console.log("\n📊 Part 1: Confidence Intervals");54console.log("-".repeat(72));5556const controlRevenueCi = meanConfidenceInterval(controlRevenuePerSession, 0.95);57const treatmentRevenueCi = meanConfidenceInterval(treatmentRevenuePerSession, 0.95);58const revenueUpliftCi = meanDiffConfidenceInterval(59  treatmentRevenuePerSession,60  controlRevenuePerSession,61  0.9562);6364console.log(65  `Control mean revenue/session:   ${controlRevenueCi.mean.toFixed(2)}  | 95% CI [${controlRevenueCi.lower.toFixed(2)}, ${controlRevenueCi.upper.toFixed(2)}]`66);67console.log(68  `Treatment mean revenue/session: ${treatmentRevenueCi.mean.toFixed(2)} | 95% CI [${treatmentRevenueCi.lower.toFixed(2)}, ${treatmentRevenueCi.upper.toFixed(2)}]`69);70console.log(71  `Treatment - control uplift:     ${revenueUpliftCi.mean.toFixed(2)}  | 95% CI [${revenueUpliftCi.lower.toFixed(2)}, ${revenueUpliftCi.upper.toFixed(2)}]`72);7374const treatmentLatencyCiZ = meanConfidenceIntervalZ(treatmentLatencyMs, 36, 0.99);75console.log(76  `Known-noise latency z-interval: ${treatmentLatencyCiZ.mean.toFixed(1)} ms | 99% CI [${treatmentLatencyCiZ.lower.toFixed(1)}, ${treatmentLatencyCiZ.upper.toFixed(1)}]`77);7879// ============================================================================80// Part 2: Conversion-rate intervals81// ============================================================================82console.log("\n✅ Part 2: Proportion Intervals");83console.log("-".repeat(72));8485const controlConversionCi = proportionConfidenceInterval(86  controlConversions.successes,87  controlConversions.total,88  0.9589);90const treatmentConversionCi = proportionConfidenceInterval(91  treatmentConversions.successes,92  treatmentConversions.total,93  0.9594);9596console.log(97  `Control conversion rate:   ${(controlConversionCi.mean * 100).toFixed(2)}% | 95% CI [${(controlConversionCi.lower * 100).toFixed(2)}%, ${(controlConversionCi.upper * 100).toFixed(2)}%]`98);99console.log(100  `Treatment conversion rate: ${(treatmentConversionCi.mean * 100).toFixed(2)}% | 95% CI [${(treatmentConversionCi.lower * 100).toFixed(2)}%, ${(treatmentConversionCi.upper * 100).toFixed(2)}%]`101);102103// ============================================================================104// Part 3: Bootstrap the mean uplift105// ============================================================================106console.log("\n♻️  Part 3: Bootstrap Uncertainty");107console.log("-".repeat(72));108109const upliftByMatchedCell = treatmentRevenuePerSession.map(110  (value, index) => value - (controlRevenuePerSession[index] ?? 0)111);112const upliftBootstrap = bootstrap(113  upliftByMatchedCell,114  (sample) => sample.reduce((sum, value) => sum + value, 0) / sample.length,115  {116    nResamples: 4000,117    seed: 48,118    confidenceLevel: 0.95,119  }120);121122console.log(123  `Bootstrap mean uplift estimate: ${upliftBootstrap.estimate.toFixed(2)} | percentile CI [${upliftBootstrap.ci[0].toFixed(2)}, ${upliftBootstrap.ci[1].toFixed(2)}]`124);125console.log(`Bootstrap resamples generated:  ${upliftBootstrap.samples.length}`);126127// ============================================================================128// Part 4: Gaussian KDE and density diagnostics129// ============================================================================130console.log("\n📈 Part 4: Gaussian KDE");131console.log("-".repeat(72));132133const controlKde = gaussian_kde(controlRevenuePerSession, { bw_method: "silverman" });134const treatmentKde = gaussian_kde(treatmentRevenuePerSession, { bw_method: "silverman" });135const probePoints = [90, 95, 100];136const controlDensity = Array.from(controlKde.evaluate(probePoints), (value) => value.toFixed(5));137const treatmentDensity = Array.from(treatmentKde.evaluate(probePoints), (value) =>138  value.toFixed(5)139);140141console.log(142  `KDE bandwidths -> control=${controlKde.bandwidth.toFixed(3)}, treatment=${treatmentKde.bandwidth.toFixed(3)}`143);144console.log(`Control KDE at [90, 95, 100]:   ${controlDensity.join(", ")}`);145console.log(`Treatment KDE at [90, 95, 100]: ${treatmentDensity.join(", ")}`);146147const densityFigure = figure({ width: 840, height: 520 });148kdeplot(tensor(controlRevenuePerSession), {149  color: "#1d4ed8",150  label: "control revenue/session",151  bw_method: "silverman",152});153kdeplot(tensor(treatmentRevenuePerSession), {154  color: "#059669",155  label: "treatment revenue/session",156  bw_method: "silverman",157});158legend();159await saveFig(`${OUTPUT_DIR}/revenue-density.svg`, { figure: densityFigure });160console.log(`Saved revenue density plot: ${OUTPUT_DIR}/revenue-density.svg`);161162// ============================================================================163// Part 5: Correct for multiple comparisons164// ============================================================================165console.log("\n🧪 Part 5: Multiple Comparisons");166console.log("-".repeat(72));167168const metricTests = [169  {170    metric: "revenue_per_session",171    pvalue: ttest_ind(tensor(controlRevenuePerSession), tensor(treatmentRevenuePerSession)).pvalue,172  },173  {174    metric: "latency_ms",175    pvalue: ttest_ind(tensor(controlLatencyMs), tensor(treatmentLatencyMs)).pvalue,176  },177  {178    metric: "csat",179    pvalue: ttest_ind(tensor(controlCsat), tensor(treatmentCsat)).pvalue,180  },181  {182    metric: "resolution_rate",183    pvalue: ttest_ind(tensor(controlResolutionRate), tensor(treatmentResolutionRate)).pvalue,184  },185];186187const rawPvalues = metricTests.map((test) => test.pvalue);188const bhCorrection = benjaminiHochberg(rawPvalues, 0.05);189const bonfCorrection = bonferroni(rawPvalues, 0.05);190191for (const [index, test] of metricTests.entries()) {192  console.log(193    `${test.metric.padEnd(20)} raw=${test.pvalue.toFixed(6)} | BH=${bhCorrection.corrected[index]?.toFixed(6)} | Bonferroni=${bonfCorrection.corrected[index]?.toFixed(6)}`194  );195}196197const correctionFigure = figure({ width: 840, height: 520 });198groupedBar(tensor([1, 2, 3, 4]), [tensor(rawPvalues), tensor(Array.from(bhCorrection.corrected))], {199  colors: ["#2563eb", "#f97316"],200  labels: ["raw p-values", "BH corrected"],201});202axhline(0.05, { color: "#991b1b", linewidth: 2, label: "alpha = 0.05" });203legend();204await saveFig(`${OUTPUT_DIR}/multiple-comparisons.svg`, { figure: correctionFigure });205console.log(`Saved multiple-comparison plot: ${OUTPUT_DIR}/multiple-comparisons.svg`);206207// ============================================================================208// Part 6: Power planning209// ============================================================================210console.log("\n🎯 Part 6: Power Analysis");211console.log("-".repeat(72));212213const observedEffectSize = Math.abs(cohenD(controlRevenuePerSession, treatmentRevenuePerSession));214const currentPower = tTestPower({215  effectSize: observedEffectSize,216  nObs: controlRevenuePerSession.length,217  alpha: 0.05,218});219const requiredSample = tTestPower({220  effectSize: observedEffectSize,221  alpha: 0.05,222  power: 0.9,223});224225console.log(`Observed effect size (|Cohen's d|): ${observedEffectSize.toFixed(3)}`);226console.log(227  `Current power at n=${controlRevenuePerSession.length}: ${currentPower.power.toFixed(3)}`228);229console.log(`Per-arm sample size for 90% power:   ${requiredSample.nObs}`);230231// ============================================================================232// Summary233// ============================================================================234console.log("\n💡 Key Takeaways");235console.log("-".repeat(72));236console.log(237  "• Use t-based confidence intervals for sample means when the population variance is unknown."238);239console.log(240  "• Use z-intervals when instrumentation or historical monitoring gives you a known noise level."241);242console.log(243  "• Bootstrap resampling gives a robust uncertainty estimate when you want fewer distributional assumptions."244);245console.log(246  "• Gaussian KDE is useful for inspecting overlap, skew, and multi-modal behavior before rollout decisions."247);248console.log(249  "• Correcting p-values matters once you inspect multiple metrics in the same experiment."250);251console.log(252  "• Power analysis turns an observed effect into a concrete follow-up sample-size plan."253);254255console.log("\n✅ Statistical Inference Playbook Complete!");256console.log("=".repeat(72));257