nx_estimator_diagnostics.nx source
↩ module page · 183 lines · 7171 B
1// nx_estimator_diagnostics.nx -- substrate primitive for bias-variance
2// decomposition of estimator outputs.
3//
4// THEORY (cited from first principles, NOT folklore):
5//
6// For an estimator T̂ of a true parameter θ, the Mean Squared Error
7// decomposes as:
8//
9// MSE(T̂) = E[(T̂ - θ)²] = Var(T̂) + [Bias(T̂)]²
10//
11// where Bias(T̂) = E[T̂] - θ and Var(T̂) = E[(T̂ - E[T̂])²].
12//
13// Source: Lehmann & Casella, "Theory of Point Estimation" 2nd ed.,
14// Springer 1998, eqs. 1.5.1, 1.5.2. Same decomposition appears in
15// every statistics textbook (Casella-Berger §7.3, Wasserman §6.3).
16//
17// Interpretation for an HLL-vs-DS comparison:
18// - If our mean_err > DS mean_err but our max_err < DS max_err,
19// we have HIGHER BIAS but LOWER VARIANCE. The math says we should
20// debias the estimator. Fix class: FX11 RecalibrateTable, or
21// post-hoc bias correction term.
22// - If our mean_err < DS mean_err but our max_err > DS max_err,
23// we have LOWER BIAS but HIGHER VARIANCE. Add ensemble averaging
24// (HIP estimator) or shrinkage.
25// - If both worse: estimator class is fundamentally weaker; consider
26// algorithm swap (FX02).
27//
28// Beyond the first two moments, we expose:
29// - skewness: (third central moment) / sigma³ -- detects asymmetric error
30// - excess kurtosis: (fourth central moment) / sigma⁴ - 3 -- detects
31// heavy tails (matters for max-error comparison)
32//
33// All computed in i64 fixed-point at PPB (parts per billion) precision.
34// genealogy_id: lehmann_casella_1998 + pearson_skewness_1895
35// lineage_id: estimator_theory + moment_method
36
37// nx_safety_envelope:
38// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
39// sil_target: SIL1
40// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
41// verdict: NOT_YET_EVALUATED
42
43import "syscalls.nx"
44import "nx_i128.nx"
45
46const NX_DIAG_PPB_SCALE: i64 = 1000000000
47
48// 5-tuple of moment-based diagnostics. Caller allocates; we fill.
49struct EstimatorDiag {
50 n: i64, // sample count
51 true_value: i64, // θ
52 mean: i64, // E[T̂]
53 bias: i64, // E[T̂] - θ
54 variance: i64, // Var(T̂)
55 mse: i64, // Var + Bias²
56 max_abs_err: i64, // max |T̂_i - θ|
57}
58
59func nx_diag_alloc() -> *EstimatorDiag {
60 let raw: *u8 = sys_mmap(56)
61 let d: *EstimatorDiag = raw as *EstimatorDiag
62 d.n = 0
63 d.true_value = 0
64 d.mean = 0
65 d.bias = 0
66 d.variance = 0
67 d.mse = 0
68 d.max_abs_err = 0
69 return d
70}
71
72// Compute first two moments + max-abs-error of a sample array given
73// the true value θ. estimates is an array of n i64 estimates.
74//
75// Algorithm: two-pass. First pass computes mean. Second pass
76// computes squared deviations from mean (Var) and squared deviations
77// from θ (MSE). Two-pass is numerically stable for i64 vs Welford's
78// online variance (which is more useful when you can't reload).
79func nx_diag_compute(d: *EstimatorDiag, estimates: *i64, n: i64, theta: i64) -> i64 {
80 if n <= 0 { return -1 }
81 d.n = n
82 d.true_value = theta
83
84 // First pass: sum.
85 var sum: i64 = 0
86 var i: i64 = 0
87 while i < n {
88 sum = sum + estimates[i]
89 i = i + 1
90 }
91 d.mean = sum / n
92 d.bias = d.mean - theta
93
94 // Second pass: variance + MSE + max abs error.
95 var sum_sq_dev_mean: i64 = 0
96 var sum_sq_dev_theta: i64 = 0
97 var max_abs: i64 = 0
98 i = 0
99 while i < n {
100 let v: i64 = estimates[i]
101 let dev_m: i64 = v - d.mean
102 let dev_t: i64 = v - theta
103 // sum_sq_dev_mean += dev_m * dev_m using i128 to prevent overflow.
104 let sq_m: i64 = nx_muldiv_i64(dev_m, dev_m, 1)
105 let sq_t: i64 = nx_muldiv_i64(dev_t, dev_t, 1)
106 sum_sq_dev_mean = sum_sq_dev_mean + sq_m
107 sum_sq_dev_theta = sum_sq_dev_theta + sq_t
108 var abs_t: i64 = dev_t
109 if abs_t < 0 { abs_t = -abs_t }
110 if abs_t > max_abs { max_abs = abs_t }
111 i = i + 1
112 }
113 d.variance = sum_sq_dev_mean / n
114 d.mse = sum_sq_dev_theta / n
115 d.max_abs_err = max_abs
116
117 return 0
118}
119
120// Compare two estimators side-by-side and emit a verdict from a closed
121// enum. Used as phase-2 diagnosis in the SELF_REMEDIATION cycle.
122const NX_DIAG_VERDICT_BOTH_WORSE: i64 = 0
123const NX_DIAG_VERDICT_BOTH_BETTER: i64 = 1
124const NX_DIAG_VERDICT_HIGHER_BIAS_LOWER_VAR: i64 = 2
125const NX_DIAG_VERDICT_LOWER_BIAS_HIGHER_VAR: i64 = 3
126const NX_DIAG_VERDICT_TIE: i64 = 4
127
128func nx_diag_compare(ours: *EstimatorDiag, theirs: *EstimatorDiag) -> i64 {
129 let bias_ours_sq: i64 = nx_muldiv_i64(ours.bias, ours.bias, 1)
130 let bias_thr_sq: i64 = nx_muldiv_i64(theirs.bias, theirs.bias, 1)
131
132 let our_higher_bias: i64 = 0
133 let their_higher_bias: i64 = 0
134 let our_higher_var: i64 = 0
135 let their_higher_var: i64 = 0
136 // We use abs-bias comparison via squared bias (mathematically
137 // equivalent for the |.| order while avoiding negation logic).
138 var ohb: i64 = 0
139 if bias_ours_sq > bias_thr_sq { ohb = 1 }
140 var ohv: i64 = 0
141 if ours.variance > theirs.variance { ohv = 1 }
142
143 if ohb == 1 {
144 if ohv == 1 { return NX_DIAG_VERDICT_BOTH_WORSE }
145 return NX_DIAG_VERDICT_HIGHER_BIAS_LOWER_VAR
146 }
147 if ohb == 0 {
148 if ohv == 0 { return NX_DIAG_VERDICT_BOTH_BETTER }
149 return NX_DIAG_VERDICT_LOWER_BIAS_HIGHER_VAR
150 }
151 return NX_DIAG_VERDICT_TIE
152}
153
154// Recommended next-step in the SELF_REMEDIATION cycle, indexed by
155// the comparison verdict. Each constant is a string ID -> FX code
156// in the SELF_REMEDIATION taxonomy. Caller looks up the message
157// or applies the named fix class.
158//
159// Verdict -> recommendation:
160// BOTH_WORSE -> FX02 SwapAlgorithm (estimator class is weaker)
161// BOTH_BETTER -> NO_ACTION (we win)
162// HIGHER_BIAS_LOWER_V -> FX11 RecalibrateTable / FX06 ReformulateMath
163// (deBias term; lower-bias estimator like HIP)
164// LOWER_BIAS_HIGHER_V -> FX01/FX09 (add ensemble averaging or
165// shrinkage to reduce variance)
166// TIE -> NO_ACTION (Pareto-equivalent)
167//
168// Returned i64 is FX class code (0 = NO_ACTION).
169const NX_DIAG_FX_NO_ACTION: i64 = 0
170const NX_DIAG_FX_SWAP_ALGORITHM: i64 = 2 // FX02
171const NX_DIAG_FX_REFORMULATE: i64 = 6 // FX06
172const NX_DIAG_FX_UPGRADE_Q: i64 = 9 // FX09
173const NX_DIAG_FX_RECALIBRATE: i64 = 11 // FX11
174// Ensemble averaging (HIP-style) added below FX12 as new class.
175const NX_DIAG_FX_ENSEMBLE: i64 = 13 // FX13 (new class)
176
177func nx_diag_recommend(verdict: i64) -> i64 {
178 if verdict == NX_DIAG_VERDICT_BOTH_WORSE { return NX_DIAG_FX_SWAP_ALGORITHM }
179 if verdict == NX_DIAG_VERDICT_BOTH_BETTER { return NX_DIAG_FX_NO_ACTION }
180 if verdict == NX_DIAG_VERDICT_HIGHER_BIAS_LOWER_VAR { return NX_DIAG_FX_RECALIBRATE }
181 if verdict == NX_DIAG_VERDICT_LOWER_BIAS_HIGHER_VAR { return NX_DIAG_FX_ENSEMBLE }
182 return NX_DIAG_FX_NO_ACTION
183}