nx_multivariate.nx source
↩ module page · 320 lines · 12153 B
1// nx_multivariate.nx -- multi-axis Monte Carlo verifier with
2// Bonferroni family-wise correction.
3//
4// User: "inputs from classifiers lead to known outputs all with
5// statistical measurement i mean monte carlo or other multivariate
6// systems prove this".
7//
8// Closes the multivariate half of the predictability cardinal.
9// Single-axis MC verification (nx_monte_carlo.nx) handles ONE
10// measurable property at a time. Real generation needs to verify
11// MANY properties simultaneously: age AND gender AND ethnicity AND
12// pose AND lighting AND style all within their target distributions.
13//
14// ===== Statistical method =========================================
15//
16// Two ways to extend single-axis t-test to multivariate:
17//
18// (a) Hotelling T-squared (Hotelling 1931): joint test with full
19// covariance matrix. Optimal power; requires N x N matrix
20// inverse for N axes.
21//
22// (b) Per-axis t-tests + Bonferroni correction (Bonferroni 1936
23// / Dunn 1961): independent test per axis; correct family-
24// wise alpha by dividing by N_axes. Conservative but no
25// matrix inverse; works for arbitrary axis count.
26//
27// We ship (b) for v1. Hotelling T² queued (nx_matrix has det /
28// inv for 2x2 + 3x3; arbitrary N needs LU-decomp primitive).
29//
30// Bonferroni: to maintain family-wise alpha = 0.05 across N axes,
31// each per-axis test uses alpha_individual = 0.05 / N. Critical
32// t-value increases accordingly. Conservative because it assumes
33// worst-case correlation; actual independence makes it overpowered,
34// but the worst case is what statistical guarantees need.
35//
36// Decision rule:
37// For each axis i:
38// verify with nx_monte_carlo using corrected critical t.
39// If ANY axis returns OFF_TARGET: overall OFF_TARGET.
40// If ALL axes return TARGET_MET: overall TARGET_MET.
41// If any axis INSUFFICIENT_SAMPLES: overall INSUFFICIENT.
42//
43// Per the bits-up + bounded-loop cardinals.
44//
45// genealogy_id: bonferroni_1936_inequality + dunn_1961_multiple_comparisons +
46// hotelling_1931_t_squared + mahalanobis_1936_distance +
47// shaffer_1995_multiple_hypothesis_testing
48// lineage_id: substrate_multivariate_v1_bonferroni
49
50// nx_safety_envelope:
51// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
52// sil_target: SIL1
53// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
54// verdict: NOT_YET_EVALUATED
55
56import "nx_syscalls.nx"
57import "nx_tier.nx"
58import "nx_loop.nx"
59import "nx_monte_carlo.nx"
60
61const NX_MV_Q10: nx_int = 1024
62
63// Critical t-values for Bonferroni-corrected per-axis tests at
64// FAMILY-WISE alpha = 0.05. Indexed by axis count:
65// 1 axis: alpha_each = 0.05 -> t* = 1.96 -> Q10 = 2007
66// 2 axes: alpha_each = 0.025 -> t* = 2.24 -> Q10 = 2294
67// 3 axes: alpha_each = 0.0167-> t* = 2.39 -> Q10 = 2447
68// 4 axes: alpha_each = 0.0125-> t* = 2.50 -> Q10 = 2560
69// 5 axes: alpha_each = 0.01 -> t* = 2.58 -> Q10 = 2641
70// 6 axes: alpha_each = 0.0083-> t* = 2.64 -> Q10 = 2703
71// 8 axes: alpha_each = 0.00625->t* = 2.73 -> Q10 = 2795
72// 10 axes: alpha_each = 0.005 -> t* = 2.81 -> Q10 = 2877
73// 16 axes: alpha_each = 0.003125-> t* = 2.97 -> Q10 = 3041
74//
75// Beyond 16 axes the table caps; callers picking that many axes
76// should consider Hotelling T² (cleaner power) or PCA dimension
77// reduction.
78
79func _mv_bonferroni_critical_t_q10(n_axes: nx_int) -> nx_int {
80 if n_axes <= 1 { return 2007 }
81 if n_axes == 2 { return 2294 }
82 if n_axes == 3 { return 2447 }
83 if n_axes == 4 { return 2560 }
84 if n_axes == 5 { return 2641 }
85 if n_axes == 6 { return 2703 }
86 if n_axes == 7 { return 2754 }
87 if n_axes == 8 { return 2795 }
88 if n_axes == 9 { return 2839 }
89 if n_axes == 10 { return 2877 }
90 if n_axes == 12 { return 2940 }
91 if n_axes == 14 { return 2993 }
92 if n_axes == 16 { return 3041 }
93 return 3041
94}
95
96// ===== Sealed-enum: MultivariateVerdict ===========================
97
98const NX_MV_TARGET_MET: nx_int = 0
99const NX_MV_OFF_TARGET: nx_int = 1
100const NX_MV_INSUFFICIENT: nx_int = 2
101const NX_MV_ERR_BAD_PARAMS: nx_int = 3
102const NX_MV_N_VERDICTS: nx_int = 4
103
104func nx_mv_verdict_is_valid(v: nx_int) -> nx_int {
105 if v < 0 { return 0 }
106 if v >= NX_MV_N_VERDICTS { return 0 }
107 return 1
108}
109
110// ===== Result envelope ============================================
111//
112// Per-axis results stored as flat i64 arrays so callers can read
113// any individual axis's measurement without struct walking.
114//
115// axis_verdicts[i] = NX_MC_TARGET_MET / OFF / INSUFFICIENT for axis i
116// axis_means[i] = empirical mean for axis i
117// axis_stds[i] = std dev for axis i
118
119struct NxMultivariateResult {
120 overall_verdict: nx_int, // NX_MV_*
121 n_axes: nx_int,
122 n_samples: nx_int,
123 axis_verdicts: *i64, // [n_axes]
124 axis_means: *i64, // [n_axes]
125 axis_stds: *i64, // [n_axes]
126 failing_axis: nx_int // -1 if all OK, else first failing axis index
127}
128
129const NX_MV_RESULT_BYTES: nx_int = 56
130
131// ===== Multi-axis sampler signature ==============================
132//
133// Caller's sampler returns ONE measurement vector at a time.
134// Signature: sample_fn(prng, out_vec, n_axes) -> caller status.
135// The substrate fills out_vec[0..n_axes) with the per-axis values.
136
137// ===== Main multivariate verify ===================================
138//
139// sample_fn: fills out_vec [n_axes] with one measurement vector
140// n_samples: how many draws per axis
141// n_axes: number of measurement axes
142// targets_q10: [n_axes] target means
143// tols_q10: [n_axes] per-axis tolerances
144// seed: PRNG seed for reproducibility
145//
146// Per-axis the substrate applies nx_monte_carlo with the Bonferroni-
147// corrected critical t. Overall verdict aggregates per-axis.
148
149func nx_multivariate_verify(
150 sample_fn: func(*i64, *i64, nx_int) -> i64,
151 n_samples: nx_int, n_axes: nx_int,
152 targets_q10: *i64, tols_q10: *i64,
153 seed: i64) -> *NxMultivariateResult {
154
155 let r: *NxMultivariateResult = sys_mmap(NX_MV_RESULT_BYTES) as *NxMultivariateResult
156 r.overall_verdict = NX_MV_ERR_BAD_PARAMS
157 r.n_axes = n_axes
158 r.n_samples = 0
159 r.failing_axis = 0 - 1
160
161 if n_samples <= 0 { return r }
162 if n_axes <= 0 { return r }
163
164 r.axis_verdicts = sys_mmap(n_axes * 8) as *i64
165 r.axis_means = sys_mmap(n_axes * 8) as *i64
166 r.axis_stds = sys_mmap(n_axes * 8) as *i64
167
168 // Per-axis Welford accumulators. We process samples in vector
169 // form (all axes together per draw) for cache efficiency.
170 let counts: *i64 = sys_mmap(n_axes * 8) as *i64
171 let means: *i64 = sys_mmap(n_axes * 8) as *i64
172 let M2s: *i64 = sys_mmap(n_axes * 8) as *i64
173 var k: nx_int = 0
174 while k < n_axes {
175 counts[k] = 0; means[k] = 0; M2s[k] = 0
176 k = k + 1
177 }
178
179 let prng: *i64 = sys_mmap(8) as *i64
180 nx_prng_init(prng, seed)
181 let vec: *i64 = sys_mmap(n_axes * 8) as *i64
182
183 var iter: nx_int = 0
184 var verdict: nx_int = NX_LOOP_RUNNING
185 let BUDGET: nx_int = n_samples
186 while verdict == NX_LOOP_RUNNING && iter < BUDGET {
187 let v_s: i64 = sample_fn(prng, vec, n_axes)
188 if v_s != 0 { verdict = NX_LOOP_ABORTED }
189 if verdict == NX_LOOP_RUNNING {
190 var ai: nx_int = 0
191 while ai < n_axes {
192 let n_new: nx_int = counts[ai] + 1
193 let delta: i64 = vec[ai] - means[ai]
194 let mean_new: i64 = means[ai] + delta / n_new
195 let delta2: i64 = vec[ai] - mean_new
196 M2s[ai] = M2s[ai] + delta * delta2
197 means[ai] = mean_new
198 counts[ai] = n_new
199 ai = ai + 1
200 }
201 }
202 iter = iter + 1
203 }
204 r.n_samples = counts[0] // all axes get the same count
205
206 if counts[0] < NX_MC_MIN_N {
207 r.overall_verdict = NX_MV_INSUFFICIENT
208 return r
209 }
210
211 // Bonferroni-corrected critical t.
212 let crit_t: nx_int = _mv_bonferroni_critical_t_q10(n_axes)
213
214 // Per-axis decision.
215 var any_off: nx_int = 0
216 var fail_idx: nx_int = 0 - 1
217 var ax: nx_int = 0
218 while ax < n_axes {
219 let n: nx_int = counts[ax]
220 let variance: i64 = M2s[ax] / (n - 1)
221 let std: i64 = nx_isqrt_q10(variance + 1)
222 r.axis_means[ax] = means[ax]
223 r.axis_stds[ax] = std
224
225 let diff: i64 = means[ax] - targets_q10[ax]
226 var abs_diff: i64 = diff
227 if abs_diff < 0 { abs_diff = 0 - abs_diff }
228
229 var axis_ok: nx_int = 0
230 if abs_diff <= tols_q10[ax] {
231 // Within tolerance; t-stat must also confirm.
232 var t_stat: i64 = 0
233 if std > 0 {
234 let sqrt_n_q10: i64 = nx_isqrt_q10(n * NX_MV_Q10)
235 t_stat = (abs_diff * sqrt_n_q10) / std
236 }
237 if t_stat < crit_t { axis_ok = 1 }
238 }
239 if axis_ok == 1 { r.axis_verdicts[ax] = NX_MC_TARGET_MET }
240 if axis_ok == 0 {
241 r.axis_verdicts[ax] = NX_MC_OFF_TARGET
242 if any_off == 0 { fail_idx = ax }
243 any_off = 1
244 }
245 ax = ax + 1
246 }
247
248 r.failing_axis = fail_idx
249 if any_off == 0 { r.overall_verdict = NX_MV_TARGET_MET }
250 if any_off == 1 { r.overall_verdict = NX_MV_OFF_TARGET }
251 return r
252}
253
254// ===== Self-test ==================================================
255//
256// Smoke: 3-axis sampler returning fixed-mean vectors with controlled
257// dispersion. Targets matched: TARGET_MET. One target far off:
258// OFF_TARGET with the failing-axis index pointing at it.
259//
260// Closed-form invariants:
261// (a) All axes within target tol -> NX_MV_TARGET_MET
262// (b) One axis far off target -> NX_MV_OFF_TARGET; failing_axis
263// = that index
264// (c) Few samples -> NX_MV_INSUFFICIENT
265// (d) Verdict gate
266
267func _mv_test_sampler(prng: *i64, out_vec: *i64, n_axes: nx_int) -> i64 {
268 // Axis 0: uniform[1000, 1099] -> mean ~1049.5
269 // Axis 1: uniform[2000, 2099] -> mean ~2049.5
270 // Axis 2: uniform[ 500, 599] -> mean ~ 549.5
271 out_vec[0] = nx_prng_range(prng, 100) + 1000
272 if n_axes > 1 { out_vec[1] = nx_prng_range(prng, 100) + 2000 }
273 if n_axes > 2 { out_vec[2] = nx_prng_range(prng, 100) + 500 }
274 return 0
275}
276
277func main() -> i64 {
278 // --- (a) 3-axis all within tolerance ---
279 let targets: *i64 = sys_mmap(3 * 8) as *i64
280 let tols: *i64 = sys_mmap(3 * 8) as *i64
281 targets[0] = 1050; targets[1] = 2050; targets[2] = 550
282 tols[0] = 50; tols[1] = 50; tols[2] = 50
283
284 let r1: *NxMultivariateResult = nx_multivariate_verify(
285 _mv_test_sampler, 500, 3, targets, tols, 0xdeadbeef)
286 if r1.overall_verdict != NX_MV_TARGET_MET { return 10 }
287 if r1.failing_axis != -1 { return 11 }
288 if r1.n_axes != 3 { return 12 }
289
290 // --- (b) One axis far off ---
291 targets[1] = 5000 // far from actual ~2050 mean
292 let r2: *NxMultivariateResult = nx_multivariate_verify(
293 _mv_test_sampler, 500, 3, targets, tols, 0xfeedface)
294 if r2.overall_verdict != NX_MV_OFF_TARGET { return 20 }
295 if r2.failing_axis != 1 { return 21 }
296 // Axis 0 and 2 should have TARGET_MET; axis 1 OFF_TARGET.
297 if r2.axis_verdicts[0] != NX_MC_TARGET_MET { return 22 }
298 if r2.axis_verdicts[1] != NX_MC_OFF_TARGET { return 23 }
299 if r2.axis_verdicts[2] != NX_MC_TARGET_MET { return 24 }
300
301 // --- (c) Insufficient samples ---
302 targets[1] = 2050 // restore good target
303 let r3: *NxMultivariateResult = nx_multivariate_verify(
304 _mv_test_sampler, 10, 3, targets, tols, 0xc0ffee)
305 if r3.overall_verdict != NX_MV_INSUFFICIENT { return 30 }
306
307 // --- (d) Bad params (n_axes = 0) ---
308 let r4: *NxMultivariateResult = nx_multivariate_verify(
309 _mv_test_sampler, 500, 0, targets, tols, 0xbadbad)
310 if r4.overall_verdict != NX_MV_ERR_BAD_PARAMS { return 40 }
311
312 // --- (e) Verdict gate ---
313 var vi: nx_int = 0
314 while vi < NX_MV_N_VERDICTS {
315 if nx_mv_verdict_is_valid(vi) != 1 { return 50 + vi }
316 vi = vi + 1
317 }
318
319 return 0
320}