nx_abstat_lib.nx source
↩ module page · 688 lines · 38479 B
1// nx_abstat_lib.nx -- THE A/B CONTRAST CORE: two arms, replicated, per-axis, family-wise corrected.
2//
3// WHY THIS EXISTS (operator 2026-09-03: "validating if we are improving in ab and multivariate methodology").
4// The three-party acceptance ledger (nx_accept) answers "has everyone signed?". It CANNOT answer "are we
5// improving?", because a referee row is a SINGLE-ARM, SINGLE-REPLICATE point grade: one capture, one number,
6// no control, no variance. percept=214 floor=151 says the frame cleared the bar. It says NOTHING about
7// whether the frame is better than the one before it, and a bar-clearing frame can be a REGRESSION from a
8// much better predecessor.
9//
10// MEASURED, NOT ASSERTED -- the near-miss that motivated this lib, 2026-09-03: four point measurements of the
11// SAME live world (/world/beach, unchanged build) returned q=5,1,5,5 and worst_ever=13,250,13,12. The spread
12// was the measuring harness own concurrent load, not the subject. A single-arm grade would have published
13// whichever sample the seat happened to take. Only replication plus a control arm separates the two.
14//
15// TWO QUESTIONS, NEVER CONFLATED. This lib answers exactly one of them and refuses to answer the other:
16// "is it above the bar?" -> referee percept vs floor. nx_accept owns it. NOT computed here.
17// "is it better than before?" -> this lib.
18// A subject can be improving and below floor, or above floor and regressing. Neither verdict is derivable
19// from the other, so a consumer must print both or say which one it is showing.
20//
21// WHY BONFERRONI. Acceptance is judged on MANY axes at once (face, skin, hair, garments, locomotion; frame
22// p50/p95/p99; draws; memory). Five axes each tested at alpha=0.05 carry a ~23% chance that at least one
23// reads "improved" by chance alone. ab_crit_t_q10 composes nx_multivariate _mv_bonferroni_critical_t_q10
24// so a multi-axis "we improved" is a claim rather than a lottery ticket. The table is NOT copied here --
25// one Bonferroni ruler in the estate, and this lib is a caller of it.
26//
27// NO SQUARE ROOT, AND THAT IS EXACT, NOT A SHORTCUT. The two-sample test is t = delta / sqrt(se2). Comparing
28// t against a critical value is algebraically identical to comparing delta^2 against crit^2 * se2 whenever
29// delta carries the sign we are testing, so this lib squares instead of rooting: no isqrt primitive, no
30// truncation of an irrational, and the comparison is exact integer arithmetic.
31//
32// Q10 IS APPLIED BY THIS LIB, AND HERE IS WHY. It composes the Welford recurrence, whose running mean uses
33// INTEGER division (mean + delta/n). On raw small integers that is lossy: pushing 4 then 5 yields mean 4,
34// not 4.5. Pushing the SAME samples scaled by Q10 (4096, 5120) yields 4608 = 4.5 exactly. So ab_push takes
35// RAW samples and scales them itself -- the caller cannot forget. The residual imprecision is one part in
36// 1024 of a sample unit, named here rather than left for a reader to discover.
37//
38// WASM-COMPATIBLE BY CONSTRUCTION: no allocation, no syscalls, no floating point in the arithmetic path. The
39// caller hands in a region of AB_WORDS i64 (an engine arena offset, or a native mmap), as nx_perf_lib does.
40// license_tier: ORIGINAL No hw writes (Rule 26).
41import "nx_multivariate.nx"
42import "nx_intlog.nx" // ilog2_1024: the estate ONE integer logarithm, for the sequential test
43
44// ---- region layout: one accumulator per ARM -----------------------------------------------------------
45const AB_A_N: i64 = 0 // samples pushed into this arm
46const AB_A_MEAN: i64 = 1 // Welford running mean, Q10
47const AB_A_M2: i64 = 2 // Welford sum of squared deviations, Q10 squared
48const AB_A_MIN: i64 = 3 // Q10
49const AB_A_MAX: i64 = 4 // Q10
50const AB_A_WORDS: i64 = 5
51
52// an EXPERIMENT region is two arms (B baseline first, then A candidate) per axis, laid out axis-major
53const AB_ARM_B: i64 = 0
54const AB_ARM_A: i64 = 1
55const AB_ARMS: i64 = 2
56const AB_MAX_AXES: i64 = 16 // _mv_bonferroni_critical_t_q10 caps its table at 16; past that it repeats
57const AB_WORDS: i64 = 160 // AB_MAX_AXES * AB_ARMS * AB_A_WORDS
58
59const AB_Q10: i64 = 1024 // NX_MC_Q10, restated as a local name only for arithmetic readability
60const AB_UNOBSERVED: i64 = 0 - 1
61const AB_HUGE: i64 = 3037000499 // floor(sqrt(i64max)): the squaring fence, so a product never wraps
62
63// ---- per-axis verdicts. FIVE, because four would force a lie ------------------------------------------
64const AB_AX_UP: i64 = 0 // significantly better AND at or above the pre-declared effect floor
65const AB_AX_DOWN: i64 = 1 // significantly worse
66const AB_AX_FLAT: i64 = 2 // the contrast does not clear the corrected critical value
67const AB_AX_BELOW_MDE: i64 = 3 // significant, but the effect is smaller than the caller declared to matter
68const AB_AX_INSUFF: i64 = 4 // cannot be tested at all: too few samples, or no within-arm variance
69const AB_AX_N: i64 = 5
70
71// ---- family verdicts ----------------------------------------------------------------------------------
72const AB_V_IMPROVED: i64 = 0
73const AB_V_REGRESSED: i64 = 1
74const AB_V_NO_CHANGE: i64 = 2
75const AB_V_INSUFFICIENT: i64 = 3
76const AB_V_N: i64 = 4
77
78// ---- named reasons an axis could not be tested. A bare INSUFF hides which one it was ------------------
79const AB_R_OK: i64 = 0
80const AB_R_FEW_B: i64 = 1 // baseline arm under the sample floor
81const AB_R_FEW_A: i64 = 2 // candidate arm under the sample floor
82const AB_R_ZEROVAR: i64 = 3 // BOTH arms have zero within-arm variance: the fixture cannot fail
83const AB_R_OVERFLOW: i64 = 4 // the squared comparison would exceed i64: refuse, never wrap
84const AB_R_NOTBINARY: i64 = 5 // a proportion test was asked about an axis holding non-0/1 samples
85const AB_R_N: i64 = 6
86
87func ab_words() -> i64 { return AB_WORDS }
88func ab_max_axes() -> i64 { return AB_MAX_AXES }
89func ab_min_n() -> i64 { return NX_MC_MIN_N } // the sample floor is the incumbent, not a fresh number
90
91// the Bonferroni-corrected critical t for k axes, from nx_multivariate table. One ruler; this is a caller.
92func ab_crit_t_q10(n_axes: i64) -> i64 { return _mv_bonferroni_critical_t_q10(n_axes) }
93
94func ab_init(r: *i64) -> i64 {
95 var i: i64 = 0
96 while i < AB_WORDS { r[i] = 0; i = i + 1 }
97 return 0
98}
99
100// slot base for (axis, arm)
101func ab_base(axis: i64, arm: i64) -> i64 {
102 return ((axis * AB_ARMS) + arm) * AB_A_WORDS
103}
104
105func ab_valid(axis: i64, arm: i64) -> i64 {
106 if axis < 0 { return 0 }
107 if axis >= AB_MAX_AXES { return 0 }
108 if arm < 0 { return 0 }
109 if arm >= AB_ARMS { return 0 }
110 return 1
111}
112
113// push ONE raw sample into (axis, arm). Scaling to Q10 happens HERE so a caller cannot forget it.
114// The Welford recurrence is inlined over the SAME arithmetic _mc_welford_update implements, because that
115// function writes through three out-pointers while this lib keeps its state in one flat caller region.
116// ab_welford_matches_incumbent proves the two agree, and the gate RUNS that proof rather than trusting
117// this comment -- an equivalence asserted in prose is not an equivalence.
118func ab_push(r: *i64, axis: i64, arm: i64, raw: i64) -> i64 {
119 if ab_valid(axis, arm) == 0 { return 0 - 1 }
120 let b: i64 = ab_base(axis, arm)
121 let s: i64 = raw * AB_Q10
122 let n_old: i64 = r[b + AB_A_N]
123 let mean_old: i64 = r[b + AB_A_MEAN]
124 let n_new: i64 = n_old + 1
125 let delta: i64 = s - mean_old
126 let mean_new: i64 = mean_old + (delta / n_new)
127 let delta2: i64 = s - mean_new
128 r[b + AB_A_N] = n_new
129 r[b + AB_A_MEAN] = mean_new
130 r[b + AB_A_M2] = r[b + AB_A_M2] + (delta * delta2)
131 if n_old == 0 { r[b + AB_A_MIN] = s; r[b + AB_A_MAX] = s; return 0 }
132 if s < r[b + AB_A_MIN] { r[b + AB_A_MIN] = s }
133 if s > r[b + AB_A_MAX] { r[b + AB_A_MAX] = s }
134 return 0
135}
136
137// THE EQUIVALENCE PROOF, callable: run the incumbent _mc_welford_update over the same two samples and
138// report whether it lands on the same (n, mean, M2). Returns 1 on agreement, 0 on divergence. A gate calls
139// this so "the arithmetic is identical" is a measurement rather than a claim in a comment.
140func ab_welford_matches_incumbent(s1: i64, s2: i64) -> i64 {
141 let on: *i64 = sys_mmap(8)
142 let om: *i64 = sys_mmap(8)
143 let o2: *i64 = sys_mmap(8)
144 _mc_welford_update(0, 0, 0, s1, on, om, o2)
145 let n1: i64 = on[0]
146 let m1: i64 = om[0]
147 let q1: i64 = o2[0]
148 _mc_welford_update(n1, m1, q1, s2, on, om, o2)
149 let reg: *i64 = sys_mmap(AB_WORDS * 8)
150 ab_init(reg)
151 ab_push(reg, 0, AB_ARM_B, s1 / AB_Q10)
152 ab_push(reg, 0, AB_ARM_B, s2 / AB_Q10)
153 let b: i64 = ab_base(0, AB_ARM_B)
154 if reg[b + AB_A_N] != on[0] { return 0 }
155 if reg[b + AB_A_MEAN] != om[0] { return 0 }
156 if reg[b + AB_A_M2] != o2[0] { return 0 }
157 return 1
158}
159
160func ab_n(r: *i64, axis: i64, arm: i64) -> i64 {
161 if ab_valid(axis, arm) == 0 { return 0 }
162 return r[ab_base(axis, arm) + AB_A_N]
163}
164
165func ab_mean_q10(r: *i64, axis: i64, arm: i64) -> i64 {
166 if ab_valid(axis, arm) == 0 { return AB_UNOBSERVED }
167 let b: i64 = ab_base(axis, arm)
168 if r[b + AB_A_N] == 0 { return AB_UNOBSERVED } // an empty arm answers UNOBSERVED, never 0
169 return r[b + AB_A_MEAN]
170}
171
172// sample variance, Q10. Bessel-corrected (n-1). UNOBSERVED below two samples -- one sample has no spread,
173// and reporting 0 there would read as "perfectly consistent" instead of "not measured".
174func ab_var_q10(r: *i64, axis: i64, arm: i64) -> i64 {
175 if ab_valid(axis, arm) == 0 { return AB_UNOBSERVED }
176 let b: i64 = ab_base(axis, arm)
177 let n: i64 = r[b + AB_A_N]
178 if n < 2 { return AB_UNOBSERVED }
179 return r[b + AB_A_M2] / ((n - 1) * AB_Q10)
180}
181
182// the contrast: candidate mean minus baseline mean, Q10. Positive = candidate is HIGHER (not yet "better").
183func ab_delta_q10(r: *i64, axis: i64) -> i64 {
184 let ma: i64 = ab_mean_q10(r, axis, AB_ARM_A)
185 let mb: i64 = ab_mean_q10(r, axis, AB_ARM_B)
186 if ma == AB_UNOBSERVED { return AB_UNOBSERVED }
187 if mb == AB_UNOBSERVED { return AB_UNOBSERVED }
188 return ma - mb
189}
190
191// why an axis cannot be tested, or AB_R_OK. Separated from the verdict so the reason survives into the report.
192func ab_axis_reason(r: *i64, axis: i64) -> i64 {
193 let nb: i64 = ab_n(r, axis, AB_ARM_B)
194 let na: i64 = ab_n(r, axis, AB_ARM_A)
195 if nb < NX_MC_MIN_N { return AB_R_FEW_B }
196 if na < NX_MC_MIN_N { return AB_R_FEW_A }
197 let vb: i64 = ab_var_q10(r, axis, AB_ARM_B)
198 let va: i64 = ab_var_q10(r, axis, AB_ARM_A)
199 // ANTI-VACUITY: zero spread in BOTH arms means the nuisance factor never varied. The test cannot fail,
200 // so it is not a test. This must never be reported as NO-CHANGE, which would read as a measurement.
201 if vb == 0 { if va == 0 { return AB_R_ZEROVAR } }
202 return AB_R_OK
203}
204
205// squared standard error of the difference, Q10: var_a/n_a + var_b/n_b (Welch, unequal variances allowed).
206func ab_se2_q10(r: *i64, axis: i64) -> i64 {
207 let vb: i64 = ab_var_q10(r, axis, AB_ARM_B)
208 let va: i64 = ab_var_q10(r, axis, AB_ARM_A)
209 if vb == AB_UNOBSERVED { return AB_UNOBSERVED }
210 if va == AB_UNOBSERVED { return AB_UNOBSERVED }
211 let nb: i64 = ab_n(r, axis, AB_ARM_B)
212 let na: i64 = ab_n(r, axis, AB_ARM_A)
213 if nb <= 0 { return AB_UNOBSERVED }
214 if na <= 0 { return AB_UNOBSERVED }
215 return (va / na) + (vb / nb)
216}
217
218// ---- the per-axis decision ----------------------------------------------------------------------------
219// higher_is_better: 1 when a LARGER value is an improvement (a percept grade), 0 when SMALLER is better
220// (a frame time in ms). Getting this wrong inverts the verdict, so it is an explicit argument with no
221// default -- an organ that guesses the direction of its own metric is the defect this parameter prevents.
222// mde_q10: the pre-declared minimum effect that MATTERS, declared BEFORE the run, never after seeing it.
223func ab_axis_verdict(r: *i64, axis: i64, n_axes: i64, higher_is_better: i64, mde_q10: i64) -> i64 {
224 if ab_axis_reason(r, axis) != AB_R_OK { return AB_AX_INSUFF }
225 let d_raw: i64 = ab_delta_q10(r, axis)
226 if d_raw == AB_UNOBSERVED { return AB_AX_INSUFF }
227 let se2: i64 = ab_se2_q10(r, axis)
228 if se2 == AB_UNOBSERVED { return AB_AX_INSUFF }
229 return ab_decide(d_raw, se2, n_axes, higher_is_better, mde_q10)
230}
231
232// THE ONE DECISION RULE -- orient, fence, square, compare, classify -- shared by the plain contrast above and
233// the CUPED contrast below, so the two estimators cannot drift onto two bars. Extracted 2026-09-04 from the
234// body of ab_axis_verdict without changing a comparison; nx_abstat_gate and nx_frameab_gate re-prove it.
235func ab_decide(d_raw: i64, se2: i64, n_axes: i64, higher_is_better: i64, mde_q10: i64) -> i64 {
236 // orient the contrast so POSITIVE always means "better", whichever way the metric runs
237 var d: i64 = d_raw
238 if higher_is_better == 0 { d = 0 - d_raw }
239 let crit: i64 = ab_crit_t_q10(n_axes)
240 var mag: i64 = d
241 if mag < 0 { mag = 0 - mag }
242 // OVERFLOW FENCE: refuse rather than wrap. A wrapped comparison answers confidently and wrongly.
243 if mag > AB_HUGE { return AB_AX_INSUFF }
244 if se2 > AB_HUGE { return AB_AX_INSUFF }
245 // t > crit <=> d^2 > (crit^2 / Q10) * se2, for d carrying the tested sign. Exact; no root taken.
246 let lhs: i64 = mag * mag
247 let rhs: i64 = ((crit * crit) / AB_Q10) * se2
248 if lhs <= rhs { return AB_AX_FLAT } // does not clear the corrected bar
249 if d < 0 { return AB_AX_DOWN } // clears it, in the WORSE direction
250 if mag < mde_q10 { return AB_AX_BELOW_MDE } // real, but smaller than the caller declared to matter
251 return AB_AX_UP
252}
253
254// ---- THE RATE AXIS NEEDS A DIFFERENT TEST, AND USING THE WRONG ONE ABSTAINS ON REAL FINDINGS -----------
255// MEASURED on nx_frameab's own fixture: every baseline frame over budget and no candidate frame over is a
256// true 100-percent-to-0 shift, and ab_axis_verdict correctly answered INSUFFICIENT -- because with BOTH arms
257// constant the two-sample t has an undefined denominator. Abstaining there is the safe direction, but it is
258// the WRONG INSTRUMENT, not a limit of the data: a 0/1 axis is a PROPORTION, and a proportion's spread is a
259// function of its mean rather than something estimated separately from it, which is exactly why the
260// degenerate case that defeats the t-test is ordinary for the two-proportion test.
261// pooled p = (k_a + k_b) / (n_a + n_b), se2 = p(1-p)(1/n_a + 1/n_b), z = (p_a - p_b) / sqrt(se2)
262// Compared by squaring, as everywhere else in this lib, so no root is taken.
263// * AN AXIS MUST DECLARE WHICH TEST IT IS: a rate graded by the continuous test silently abstains on its
264// strongest findings, and a continuous axis graded by the proportion test is simply invalid arithmetic.
265// That is why this is a separate named function and not a flag inside ab_axis_verdict.
266
267// a proportion axis must actually hold 0/1 samples; min and max are already tracked, so this is free.
268// A caller that pushed continuous data into a rate axis gets a REFUSAL, never a confident wrong number.
269func ab_is_binary(r: *i64, axis: i64, arm: i64) -> i64 {
270 if ab_valid(axis, arm) == 0 { return 0 }
271 let b: i64 = ab_base(axis, arm)
272 if r[b + AB_A_N] == 0 { return 0 }
273 let lo: i64 = r[b + AB_A_MIN]
274 let hi: i64 = r[b + AB_A_MAX]
275 if lo != 0 { if lo != AB_Q10 { return 0 } }
276 if hi != 0 { if hi != AB_Q10 { return 0 } }
277 return 1
278}
279
280// pooled proportion of the two arms, Q10. UNOBSERVED when either arm is empty.
281func ab_pooled_p_q10(r: *i64, axis: i64) -> i64 {
282 let nb: i64 = ab_n(r, axis, AB_ARM_B)
283 let na: i64 = ab_n(r, axis, AB_ARM_A)
284 if nb <= 0 { return AB_UNOBSERVED }
285 if na <= 0 { return AB_UNOBSERVED }
286 let pb: i64 = ab_mean_q10(r, axis, AB_ARM_B)
287 let pa: i64 = ab_mean_q10(r, axis, AB_ARM_A)
288 if pb == AB_UNOBSERVED { return AB_UNOBSERVED }
289 if pa == AB_UNOBSERVED { return AB_UNOBSERVED }
290 return ((pa * na) + (pb * nb)) / (na + nb)
291}
292
293// squared standard error of the difference in proportions, Q10.
294func ab_prop_se2_q10(r: *i64, axis: i64) -> i64 {
295 let p: i64 = ab_pooled_p_q10(r, axis)
296 if p == AB_UNOBSERVED { return AB_UNOBSERVED }
297 let nb: i64 = ab_n(r, axis, AB_ARM_B)
298 let na: i64 = ab_n(r, axis, AB_ARM_A)
299 if nb <= 0 { return AB_UNOBSERVED }
300 if na <= 0 { return AB_UNOBSERVED }
301 // p(1-p) x (n_a + n_b) / (n_a x n_b), all in Q10; the multiply is ordered to keep the numerator large
302 // so a small proportion does not truncate its own variance to zero before the division.
303 return (p * (AB_Q10 - p) * (na + nb)) / (AB_Q10 * na * nb)
304}
305
306// THE REASON MUST COME FROM THE TEST THAT WAS ACTUALLY RUN. ab_axis_reason answers for the CONTINUOUS test,
307// and printing it beside a proportion verdict produces a line that contradicts itself -- measured live on
308// nx_frameab's own fixture, where an axis read "verdict=UP reason=ZERO-VARIANCE-BOTH-ARMS-NOT-A-TEST" and a
309// reader could not tell which half to believe.
310// * A VERDICT AND ITS REASON DERIVED FROM DIFFERENT INSTRUMENTS IS A COMPOUND ASSERTION THAT NAMES THE WRONG
311// SUBJECT, AND IT IS WORSE THAN NO REASON AT ALL, BECAUSE THE READER TRUSTS THE FIELD.
312func ab_prop_reason(r: *i64, axis: i64) -> i64 {
313 let nb: i64 = ab_n(r, axis, AB_ARM_B)
314 let na: i64 = ab_n(r, axis, AB_ARM_A)
315 if nb < NX_MC_MIN_N { return AB_R_FEW_B }
316 if na < NX_MC_MIN_N { return AB_R_FEW_A }
317 if ab_is_binary(r, axis, AB_ARM_B) == 0 { return AB_R_NOTBINARY }
318 if ab_is_binary(r, axis, AB_ARM_A) == 0 { return AB_R_NOTBINARY }
319 let se2: i64 = ab_prop_se2_q10(r, axis)
320 if se2 == AB_UNOBSERVED { return AB_R_FEW_B }
321 // the one genuine abstention for a proportion: both arms the SAME constant, so there is no contrast
322 if se2 <= 0 { if ab_delta_q10(r, axis) == 0 { return AB_R_ZEROVAR } }
323 return AB_R_OK
324}
325
326// the per-axis decision for a PROPORTION axis. Same five states, same pre-declared MDE, same orientation
327// argument -- only the variance model differs, and a non-binary axis is REFUSED rather than approximated.
328func ab_prop_verdict(r: *i64, axis: i64, n_axes: i64, higher_is_better: i64, mde_q10: i64) -> i64 {
329 let nb: i64 = ab_n(r, axis, AB_ARM_B)
330 let na: i64 = ab_n(r, axis, AB_ARM_A)
331 if nb < NX_MC_MIN_N { return AB_AX_INSUFF }
332 if na < NX_MC_MIN_N { return AB_AX_INSUFF }
333 if ab_is_binary(r, axis, AB_ARM_B) == 0 { return AB_AX_INSUFF }
334 if ab_is_binary(r, axis, AB_ARM_A) == 0 { return AB_AX_INSUFF }
335 let d_raw: i64 = ab_delta_q10(r, axis)
336 if d_raw == AB_UNOBSERVED { return AB_AX_INSUFF }
337 var d: i64 = d_raw
338 if higher_is_better == 0 { d = 0 - d_raw }
339 let se2: i64 = ab_prop_se2_q10(r, axis)
340 if se2 == AB_UNOBSERVED { return AB_AX_INSUFF }
341 // BOTH ARMS IDENTICAL AND CONSTANT is the one case a proportion test still cannot speak to: the pooled
342 // p is 0 or 1, se2 is 0, and there is genuinely no contrast. That is a real abstention, not a wrong tool.
343 if se2 <= 0 { if d == 0 { return AB_AX_INSUFF } }
344 let crit: i64 = ab_crit_t_q10(n_axes)
345 var mag: i64 = d
346 if mag < 0 { mag = 0 - mag }
347 if mag > AB_HUGE { return AB_AX_INSUFF }
348 if se2 > AB_HUGE { return AB_AX_INSUFF }
349 let lhs: i64 = mag * mag
350 let rhs: i64 = ((crit * crit) / AB_Q10) * se2
351 if lhs <= rhs { return AB_AX_FLAT }
352 if d < 0 { return AB_AX_DOWN }
353 if mag < mde_q10 { return AB_AX_BELOW_MDE }
354 return AB_AX_UP
355}
356
357// ---- the family verdict -------------------------------------------------------------------------------
358// REGRESSED DOMINATES IMPROVED, deliberately. A change that lifts four axes and drops one is a regression
359// until someone adjudicates the dropped one; the fail direction of this organ is to withhold good news.
360// An axis that could not be tested does NOT silently vanish: if no axis was testable the family answers
361// INSUFFICIENT, because a verdict computed over zero testable axes is the pass-on-the-empty-set defect.
362func ab_family_verdict(axv: *i64, n_axes: i64) -> i64 {
363 if n_axes <= 0 { return AB_V_INSUFFICIENT }
364 var tested: i64 = 0
365 var up: i64 = 0
366 var down: i64 = 0
367 var i: i64 = 0
368 while i < n_axes {
369 let v: i64 = axv[i]
370 if v != AB_AX_INSUFF { tested = tested + 1 }
371 if v == AB_AX_UP { up = up + 1 }
372 if v == AB_AX_DOWN { down = down + 1 }
373 i = i + 1
374 }
375 if tested == 0 { return AB_V_INSUFFICIENT }
376 if down > 0 { return AB_V_REGRESSED }
377 if up > 0 { return AB_V_IMPROVED }
378 return AB_V_NO_CHANGE
379}
380
381// how many axes were genuinely testable -- the denominator every consumer must print beside the verdict,
382// or an abstaining family reads as a confident NO-CHANGE.
383func ab_tested_axes(axv: *i64, n_axes: i64) -> i64 {
384 var t: i64 = 0
385 var i: i64 = 0
386 while i < n_axes { if axv[i] != AB_AX_INSUFF { t = t + 1 } i = i + 1 }
387 return t
388}
389
390// ---- A-PRIORI SAMPLE SIZE, the half of an honest improvement claim that must be printed BEFORE the arms are read
391// SOURCE, read from the mirrored bytes (Kohavi, Henne, Sommerfield, Practical Guide to Controlled Experiments on the
392// Web, KDD 2007, ref ge-expguide): for 95 percent confidence and 90 percent power the number of users per variant is
393// approximately n = (4 r sigma / Delta)^2, r the number of variants, sigma the standard deviation of the OEC and Delta
394// the minimum difference to detect, in the SAME unit as sigma. The paper's own worked examples are the known answers
395// this function is gated against: 4 x 2 x 30 / (3.75 x 0.05) squared is over 1.6 million users for a 5 percent revenue
396// change at a 30 dollar standard deviation on a 3.75 dollar mean; the checkout example 4 x 2 x 0.5 / (0.5 x 0.05)
397// squared is 25,600. The quotient is rounded to nearest before squaring (floor under-states n, which is the
398// permissive direction); a quotient past AB_HUGE returns AB_UNOBSERVED rather than wrapping. Unit-agnostic: the
399// caller passes sigma and Delta in one integer unit fine enough to make both integral.
400const AB_SS_COEF: i64 = 4 // the paper's coefficient for 95 percent confidence and 90 percent power
401func ab_sample_size_aa(r_variants: i64, sigma_units: i64, mde_units: i64) -> i64 {
402 if r_variants <= 0 { return AB_UNOBSERVED }
403 if sigma_units <= 0 { return AB_UNOBSERVED }
404 if mde_units <= 0 { return AB_UNOBSERVED }
405 let num: i64 = AB_SS_COEF * r_variants * sigma_units
406 let q: i64 = (num * 2 + mde_units) / (2 * mde_units) // round to nearest
407 if q > AB_HUGE { return AB_UNOBSERVED }
408 return q * q
409}
410// The A/A half of the same contract is NOT a library function: a rejection RATE is a measurement over many
411// synthetic experiments driven by a seeded generator, which belongs in the gate that calibrates this lib
412// (nx_abstat_aa_gate) -- the paper's bar is that the null is rejected about 5 percent of the time at 95 percent
413// confidence, and that bar is asserted there with a binomial envelope, never here as a constant.
414
415// =====================================================================================================
416// (2) ALWAYS-VALID INFERENCE -- the mixture sequential probability ratio test (Robbins 1970) as deployed by
417// Johari, Pekelis and Walsh (ref ge-alwaysvalid). Read from the mirrored paper: the mSPRT statistic is the
418// mixture likelihood ratio Lambda_n = integral over theta of (f_theta(s_n)/f_0(s_n))^n dH(theta), eq (7); the
419// test stops the first time Lambda_n reaches 1/alpha, eq (8); and Theorem 1 turns any such sequential test
420// into an always-valid p-value process, p_n = inf{alpha : the alpha-test has stopped by n}, i.e. the running
421// minimum of 1/Lambda_n, which controls Type I error at ANY data-dependent stopping time -- the property the
422// paper measures the fixed-horizon p-value as lacking under continuous monitoring ("even with 10000 samples
423// Type I error can easily increase fivefold").
424// THE NORMAL CLOSED FORM, derived here because the paper fixes the family (single-parameter exponential; the
425// normal case G = N(0, sigma^2), H = N(0, tau^2) in Section 5.5) and leaves the Gaussian integral to the
426// reader. For a contrast estimate d with squared standard error V (the paper's sigma^2/n), under H0 E[d] = 0,
427// the likelihood ratio at effect theta is exp((theta d - theta^2/2)/V); integrating against N(0, tau^2)
428// completes the square:
429// ln Lambda = (1/2) ln( V / (V + tau^2) ) + tau^2 d^2 / ( 2 V (V + tau^2) ).
430// The first term is the mixture's price for not knowing the effect -- negative, so Lambda starts below one and
431// an empty record never rejects; the second grows with the standardised effect and tends to d^2/(2V) as tau^2
432// dominates V. Reject at level alpha when ln Lambda >= ln(1/alpha).
433// THIS IS THE PER-LOOK STATISTIC. The always-valid p-value is the running minimum of 1/Lambda, so a caller
434// that monitors continuously keeps the running MAXIMUM of this value across its looks and compares that to
435// ab_av_bar_q10 -- a p-value that could rise again after falling is not a valid one (Theorem 1).
436// tau^2 IS DERIVED, NEVER TYPED: the paper chooses the mixing distribution to match the effects the user
437// cares to detect, and this estate already declares that effect as the pre-declared MDE, so
438// ab_av_tau2_from_mde_q10 squares it. Integer Q10 throughout; the natural log is the estate's integer log2
439// times ln 2, whose rounding (1 part in 710) is the imprecision this lives with, named here.
440// =====================================================================================================
441const AB_LN2_Q10: i64 = 710 // ln 2 = 0.693147 x 1024 = 709.78, rounded: the ONE ln 2 in this lib
442const AB_LOG2_ONE_Q10: i64 = 10240 // ilog2_1024(1024): log2 of a Q10 one, in Q10
443const AB_I64_MAX: i64 = 9223372036854775807 // the exact product fence: a x b is safe iff a <= MAX / b
444const AB_PERMIL: i64 = 1000
445
446// num / den in Q10 for den > 0, without forming num x Q10 (which overflows first): integer part, then the
447// remainder scaled. Truncates toward zero like every division in this lib. UNOBSERVED when den <= 0.
448func ab_ratio_q10(num: i64, den: i64) -> i64 {
449 if den <= 0 { return AB_UNOBSERVED }
450 var n: i64 = num
451 var d: i64 = den
452 // scale both operands down together until the remainder can carry a Q10 factor: a quotient of two large
453 // moments loses nothing a Q10 result can express, and a refusal here would blind the estimator on exactly
454 // the big samples where it matters (measured: 400,000 pairs put the co-moment past the plain fence)
455 while d > AB_I64_MAX / AB_Q10 { n = n / 2; d = d / 2 }
456 return (n / d) * AB_Q10 + ((n % d) * AB_Q10) / d
457}
458
459// a x b for a Q10 a and any Q10-scaled b, result at b's scale: split a into whole and fractional Q10 parts so
460// the product never forms a x b directly (the whole part is fenced by AB_I64_MAX / |b|).
461func ab_mul_q10(a: i64, b: i64) -> i64 {
462 var bm: i64 = b
463 if bm < 0 { bm = 0 - b }
464 if bm == 0 { return 0 }
465 var aw: i64 = a / AB_Q10
466 if aw < 0 { aw = 0 - aw }
467 if aw > AB_I64_MAX / bm { return AB_UNOBSERVED }
468 return (a / AB_Q10) * b + ((a % AB_Q10) * b) / AB_Q10
469}
470
471// natural log of a POSITIVE Q10 value, Q10 out: ln(x) = (log2(x_q10) - log2(1024)) x ln 2.
472func ab_ln_q10(x_q10: i64) -> i64 {
473 if x_q10 <= 0 { return AB_UNOBSERVED }
474 return ((ilog2_1024(x_q10) - AB_LOG2_ONE_Q10) * AB_LN2_Q10) / AB_Q10
475}
476
477// the mixing variance from the pre-declared minimum detectable effect: tau^2 = mde^2. Q10 in, Q10 out.
478func ab_av_tau2_from_mde_q10(mde_q10: i64) -> i64 {
479 if mde_q10 <= 0 { return AB_UNOBSERVED }
480 if mde_q10 > AB_HUGE { return AB_UNOBSERVED }
481 return (mde_q10 * mde_q10) / AB_Q10
482}
483
484// the rejection bar ln(1/alpha) in Q10, alpha in permil: 50 -> ln 20 = 2.9957 -> 3068.
485func ab_av_bar_q10(alpha_permil: i64) -> i64 {
486 if alpha_permil <= 0 { return AB_UNOBSERVED }
487 if alpha_permil >= AB_PERMIL { return AB_UNOBSERVED }
488 return ab_ln_q10((AB_PERMIL * AB_Q10) / alpha_permil)
489}
490
491// ab_always_valid: ln Lambda for ONE look at the contrast. delta_q10 from ab_delta_q10 (or the CUPED delta),
492// se2_q10 from ab_se2_q10 (or the CUPED se2), tau2_q10 from ab_av_tau2_from_mde_q10. Q10 in, Q10 out.
493// UNOBSERVED when a variance is non-positive or a product would exceed i64 -- refuse, never wrap, because a
494// wrapped statistic rejects confidently.
495func ab_always_valid(delta_q10: i64, se2_q10: i64, tau2_q10: i64) -> i64 {
496 if se2_q10 <= 0 { return AB_UNOBSERVED }
497 if tau2_q10 <= 0 { return AB_UNOBSERVED }
498 let vt: i64 = se2_q10 + tau2_q10
499 let ln_v: i64 = ab_ln_q10(se2_q10)
500 let ln_vt: i64 = ab_ln_q10(vt)
501 if ln_v == AB_UNOBSERVED { return AB_UNOBSERVED }
502 if ln_vt == AB_UNOBSERVED { return AB_UNOBSERVED }
503 let lnterm: i64 = (ln_v - ln_vt) / 2
504 var d: i64 = delta_q10
505 if d < 0 { d = 0 - delta_q10 }
506 if d == 0 { return lnterm }
507 // numerator tau^2 d^2, Q10, two fenced products
508 if tau2_q10 > AB_I64_MAX / d { return AB_UNOBSERVED }
509 let a: i64 = (tau2_q10 * d) / AB_Q10
510 if a > AB_I64_MAX / d { return AB_UNOBSERVED }
511 let num: i64 = (a * d) / AB_Q10
512 // denominator 2 V (V + tau^2), Q10, one fenced product
513 if se2_q10 > AB_I64_MAX / vt { return AB_UNOBSERVED }
514 let den: i64 = 2 * ((se2_q10 * vt) / AB_Q10)
515 if den <= 0 { return AB_UNOBSERVED }
516 let term: i64 = ab_ratio_q10(num, den)
517 if term == AB_UNOBSERVED { return AB_UNOBSERVED }
518 return lnterm + term
519}
520
521// =====================================================================================================
522// (3) CUPED -- Controlled-experiment Using Pre-Experiment Data (Deng, Xu, Kohavi, Walker, WSDM 2013, ref
523// ge-cuped). Read from the mirrored paper, Section 3.2: the control-variate estimator is
524// Ycv = Y - theta X + theta E[X] (eq 3)
525// unbiased for E[Y] for ANY constant theta; its variance is minimised at
526// theta = cov(Y, X) / var(X) (eq 4)
527// where it equals var(Y)(1 - rho^2), rho the correlation of Y and X (eq 5)
528// -- "the variance is reduced by a factor of rho^2" -- and the same metric from the pre-experiment period is
529// the most effective X the paper found ("cut variance by about 50 percent" on Bing). THE SUBTLETY THE PAPER
530// NAMES AND THIS CODE ENFORCES BY CONSTRUCTION: "the same theta has to be used for both control and
531// treatment. The simplest way to estimate it is from the pooled population of control and treatment." So
532// theta is read from a POOLED accumulator that every push feeds, never from an arm, and the adjusted
533// contrast is delta_cv = (mean_Y_A - mean_Y_B) - theta (mean_X_A - mean_X_B) -- the E[X] terms cancel.
534// A bivariate Welford recurrence keeps n, both means, both second central moments and the co-moment per arm
535// and for the pool; Q10 like the rest of this lib, scaled at the push so a caller cannot forget.
536// =====================================================================================================
537const AB_CV_N: i64 = 0
538const AB_CV_MX: i64 = 1 // running mean of the covariate X, Q10
539const AB_CV_MY: i64 = 2 // running mean of the metric Y, Q10
540const AB_CV_M2X: i64 = 3 // sum of squared deviations of X, Q10 squared
541const AB_CV_M2Y: i64 = 4 // sum of squared deviations of Y, Q10 squared
542const AB_CV_CXY: i64 = 5 // co-moment sum (X - mx)(Y - my), Q10 squared
543const AB_CV_WORDS: i64 = 6
544const AB_CV_POOL: i64 = 2 // the third accumulator: both arms together, the paper's pooled theta
545const AB_CV_ACCS: i64 = 3 // B, A, POOL
546const AB_CV_TOTAL: i64 = 18 // AB_CV_ACCS x AB_CV_WORDS
547
548func ab_cuped_words() -> i64 { return AB_CV_TOTAL }
549func ab_cuped_init(cv: *i64) -> i64 {
550 var i: i64 = 0
551 while i < AB_CV_TOTAL { cv[i] = 0; i = i + 1 }
552 return 0
553}
554func ab_cv_base(acc: i64) -> i64 { return acc * AB_CV_WORDS }
555
556// one bivariate Welford step on accumulator acc over Q10-scaled (x, y): the mean recurrence is the one
557// ab_push inlines; the co-moment adds (x - mean_x_old)(y - mean_y_new), the pairwise form of the M2 update.
558func ab_cv_step(cv: *i64, acc: i64, x: i64, y: i64) -> i64 {
559 let b: i64 = ab_cv_base(acc)
560 let n_new: i64 = cv[b + AB_CV_N] + 1
561 let dx: i64 = x - cv[b + AB_CV_MX]
562 let dy: i64 = y - cv[b + AB_CV_MY]
563 let mx_new: i64 = cv[b + AB_CV_MX] + (dx / n_new)
564 let my_new: i64 = cv[b + AB_CV_MY] + (dy / n_new)
565 let dx2: i64 = x - mx_new
566 let dy2: i64 = y - my_new
567 cv[b + AB_CV_N] = n_new
568 cv[b + AB_CV_MX] = mx_new
569 cv[b + AB_CV_MY] = my_new
570 cv[b + AB_CV_M2X] = cv[b + AB_CV_M2X] + (dx * dx2)
571 cv[b + AB_CV_M2Y] = cv[b + AB_CV_M2Y] + (dy * dy2)
572 cv[b + AB_CV_CXY] = cv[b + AB_CV_CXY] + (dx * dy2)
573 return n_new
574}
575
576// ab_cuped: push ONE paired observation -- the metric y and its pre-experiment covariate x, RAW units -- into
577// an arm. Scaled to Q10 here, as ab_push does. Feeds the arm AND the pool. Returns the arm's new n; -1 on a
578// bad arm.
579func ab_cuped(cv: *i64, arm: i64, y_raw: i64, x_raw: i64) -> i64 {
580 if arm < 0 { return 0 - 1 }
581 if arm >= AB_ARMS { return 0 - 1 }
582 let x: i64 = x_raw * AB_Q10
583 let y: i64 = y_raw * AB_Q10
584 ab_cv_step(cv, AB_CV_POOL, x, y)
585 return ab_cv_step(cv, arm, x, y)
586}
587func ab_cuped_n(cv: *i64, acc: i64) -> i64 {
588 if acc < 0 { return 0 }
589 if acc >= AB_CV_ACCS { return 0 }
590 return cv[ab_cv_base(acc) + AB_CV_N]
591}
592
593// theta = cov(Y, X) / var(X) over the POOL, Q10 (eq 4). UNOBSERVED under two pooled samples or when X never
594// varied -- a constant covariate carries no information, and a division by zero would read as infinite gain.
595func ab_cuped_theta_q10(cv: *i64) -> i64 {
596 let b: i64 = ab_cv_base(AB_CV_POOL)
597 if cv[b + AB_CV_N] < 2 { return AB_UNOBSERVED }
598 let m2x: i64 = cv[b + AB_CV_M2X]
599 if m2x <= 0 { return AB_UNOBSERVED }
600 return ab_ratio_q10(cv[b + AB_CV_CXY], m2x)
601}
602
603// rho^2 = cov^2 / (var X var Y) over the pool, Q10 -- the paper's variance-reduction factor (eq 5), reported
604// so the reduction a consumer prints is measured from the same moments the adjustment used.
605func ab_cuped_rho2_q10(cv: *i64) -> i64 {
606 let b: i64 = ab_cv_base(AB_CV_POOL)
607 if cv[b + AB_CV_N] < 2 { return AB_UNOBSERVED }
608 let m2x: i64 = cv[b + AB_CV_M2X]
609 let m2y: i64 = cv[b + AB_CV_M2Y]
610 if m2x <= 0 { return AB_UNOBSERVED }
611 if m2y <= 0 { return AB_UNOBSERVED }
612 var c: i64 = cv[b + AB_CV_CXY]
613 if c < 0 { c = 0 - c }
614 let r1: i64 = ab_ratio_q10(c, m2x)
615 let r2: i64 = ab_ratio_q10(c, m2y)
616 if r1 == AB_UNOBSERVED { return AB_UNOBSERVED }
617 if r2 == AB_UNOBSERVED { return AB_UNOBSERVED }
618 return ab_mul_q10(r1, r2)
619}
620
621// CUPED-adjusted mean of an arm, Q10: mean_Y - theta (mean_X - mean_X_pool). E[X] is the pooled mean, so both
622// arms are referenced to the same point and the reference cancels in the contrast (eq 3).
623func ab_cuped_mean_q10(cv: *i64, arm: i64, theta_q10: i64) -> i64 {
624 if arm < 0 { return AB_UNOBSERVED }
625 if arm >= AB_ARMS { return AB_UNOBSERVED }
626 let b: i64 = ab_cv_base(arm)
627 if cv[b + AB_CV_N] == 0 { return AB_UNOBSERVED }
628 let dx: i64 = cv[b + AB_CV_MX] - cv[ab_cv_base(AB_CV_POOL) + AB_CV_MX]
629 let adj: i64 = ab_mul_q10(theta_q10, dx)
630 if adj == AB_UNOBSERVED { return AB_UNOBSERVED }
631 return cv[b + AB_CV_MY] - adj
632}
633
634// the CUPED contrast: adjusted candidate mean minus adjusted baseline mean, Q10
635func ab_cuped_delta_q10(cv: *i64, theta_q10: i64) -> i64 {
636 let ma: i64 = ab_cuped_mean_q10(cv, AB_ARM_A, theta_q10)
637 let mb: i64 = ab_cuped_mean_q10(cv, AB_ARM_B, theta_q10)
638 if ma == AB_UNOBSERVED { return AB_UNOBSERVED }
639 if mb == AB_UNOBSERVED { return AB_UNOBSERVED }
640 return ma - mb
641}
642
643// sample variance of the adjusted metric in an arm, Q10, Bessel-corrected:
644// var(Y - theta X) = (M2Y - 2 theta CXY + theta^2 M2X) / (n - 1), moments Q10 squared, theta Q10.
645// A rounding residue below zero is refused rather than published as a negative variance.
646func ab_cuped_var_q10(cv: *i64, arm: i64, theta_q10: i64) -> i64 {
647 if arm < 0 { return AB_UNOBSERVED }
648 if arm >= AB_ARMS { return AB_UNOBSERVED }
649 let b: i64 = ab_cv_base(arm)
650 let n: i64 = cv[b + AB_CV_N]
651 if n < 2 { return AB_UNOBSERVED }
652 let tc: i64 = ab_mul_q10(theta_q10, cv[b + AB_CV_CXY])
653 let t2: i64 = ab_mul_q10(theta_q10, theta_q10)
654 if tc == AB_UNOBSERVED { return AB_UNOBSERVED }
655 if t2 == AB_UNOBSERVED { return AB_UNOBSERVED }
656 let t2m: i64 = ab_mul_q10(t2, cv[b + AB_CV_M2X])
657 if t2m == AB_UNOBSERVED { return AB_UNOBSERVED }
658 let m2cv: i64 = cv[b + AB_CV_M2Y] - 2 * tc + t2m
659 if m2cv < 0 { return AB_UNOBSERVED }
660 return m2cv / ((n - 1) * AB_Q10)
661}
662
663// squared standard error of the CUPED contrast, Q10 (the Welch form ab_se2_q10 uses)
664func ab_cuped_se2_q10(cv: *i64, theta_q10: i64) -> i64 {
665 let va: i64 = ab_cuped_var_q10(cv, AB_ARM_A, theta_q10)
666 let vb: i64 = ab_cuped_var_q10(cv, AB_ARM_B, theta_q10)
667 if va == AB_UNOBSERVED { return AB_UNOBSERVED }
668 if vb == AB_UNOBSERVED { return AB_UNOBSERVED }
669 let na: i64 = ab_cuped_n(cv, AB_ARM_A)
670 let nb: i64 = ab_cuped_n(cv, AB_ARM_B)
671 if na <= 0 { return AB_UNOBSERVED }
672 if nb <= 0 { return AB_UNOBSERVED }
673 return (va / na) + (vb / nb)
674}
675
676// the CUPED verdict: ab_decide -- the SAME rule as ab_axis_verdict -- over the adjusted contrast. Returns AB_AX_*.
677// INSUFF below the sample floor, on an unobservable contrast, or when both adjusted arms have zero spread
678// (the anti-vacuity refusal ab_axis_reason applies to the plain contrast).
679func ab_cuped_verdict(cv: *i64, theta_q10: i64, n_axes: i64, higher_is_better: i64, mde_q10: i64) -> i64 {
680 if ab_cuped_n(cv, AB_ARM_B) < NX_MC_MIN_N { return AB_AX_INSUFF }
681 if ab_cuped_n(cv, AB_ARM_A) < NX_MC_MIN_N { return AB_AX_INSUFF }
682 let d_raw: i64 = ab_cuped_delta_q10(cv, theta_q10)
683 if d_raw == AB_UNOBSERVED { return AB_AX_INSUFF }
684 let se2: i64 = ab_cuped_se2_q10(cv, theta_q10)
685 if se2 == AB_UNOBSERVED { return AB_AX_INSUFF }
686 if se2 == 0 { return AB_AX_INSUFF }
687 return ab_decide(d_raw, se2, n_axes, higher_is_better, mde_q10)
688}