nx_monte_carlo.nx source
↩ module page · 287 lines · 10201 B
1// nx_monte_carlo.nx -- Monte Carlo verifier with statistical bounds.
2//
3// First substrate primitive grounded in CLASSICAL STATISTICAL
4// THEORY per cardinal feedback-loras-and-negatives-are-patches-not-
5// systems:
6//
7// User: "inputs from classifiers lead to known outputs all with
8// statistical measurement i mean monte carlo or other multivariate
9// systems prove this"
10//
11// The substrate doesn't invent predictability. It APPLIES the
12// techniques (Monte Carlo, multivariate statistics, classification
13// theory) the AI field abandoned for deep-learning opacity.
14//
15// ===== What this primitive does ==================================
16//
17// Given:
18// - sample_fn : generates one measurement (i64 scalar)
19// - n_samples : how many to draw
20// - target_q10 : the target mean (in caller's chosen units)
21// - tol_q10 : tolerance bound (output mean must be within
22// target +/- tol)
23// - seed : PRNG seed for reproducibility
24//
25// Returns:
26// - empirical mean
27// - empirical variance (Welford 1962 online algorithm)
28// - verdict: TARGET_MET / OFF_TARGET / INSUFFICIENT_SAMPLES
29// - the t-statistic for the test (caller can compute confidence
30// intervals + decide their own threshold if 95% default is
31// wrong for their use)
32//
33// ===== Statistical method =========================================
34//
35// One-sample t-test (Student 1908, Welch 1947 variant for unequal
36// variance not needed here since we test against fixed target):
37//
38// t = (mean - target) / (sd / sqrt(n))
39//
40// |t| < critical_value for given confidence + (n-1) dof -> mean is
41// within target with that confidence. We default to 95% CI which
42// has critical_t ~ 1.96 for large n (infinite dof) and ~2.05 for
43// n=30. In Q10 the critical value at 95% is ~2010 for n >= 30.
44//
45// Welford 1962 online mean/variance avoids numerical-precision
46// problems that naive accumulation of sum + sum_sq has. Each
47// sample contributes:
48//
49// delta = x - mean
50// mean += delta / n
51// delta2 = x - mean (post-update)
52// M2 += delta * delta2
53// var = M2 / (n - 1) (Bessel-corrected sample variance)
54//
55// Per the bits-up + bounded-loop cardinals.
56//
57// genealogy_id: metropolis_ulam_1949_monte_carlo +
58// student_1908_t_distribution +
59// welford_1962_online_variance +
60// shewhart_1924_statistical_process_control
61// lineage_id: substrate_monte_carlo_v1
62
63// nx_safety_envelope:
64// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
65// sil_target: SIL1
66// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
67// verdict: NOT_YET_EVALUATED
68
69import "nx_syscalls.nx"
70import "nx_tier.nx"
71import "nx_loop.nx"
72import "nx_isqrt.nx"
73import "nx_prng.nx"
74
75// ===== Sealed-enum: MonteCarloVerdict =============================
76
77const NX_MC_TARGET_MET: nx_int = 0
78const NX_MC_OFF_TARGET: nx_int = 1
79const NX_MC_INSUFFICIENT_SAMPLES: nx_int = 2
80const NX_MC_ERR_BAD_PARAMS: nx_int = 3
81const NX_MC_N_VERDICTS: nx_int = 4
82
83func nx_mc_verdict_is_valid(v: nx_int) -> nx_int {
84 if v < 0 { return 0 }
85 if v >= NX_MC_N_VERDICTS { return 0 }
86 return 1
87}
88
89// ===== Result envelope ============================================
90
91struct NxMonteCarloResult {
92 verdict: nx_int, // NX_MC_TARGET_MET / OFF / INSUFFICIENT / ERR
93 n_samples: nx_int, // actual samples drawn
94 mean_q10: nx_int, // empirical mean (Q10)
95 variance_q10: nx_int, // empirical variance (Q20 / Q10 = Q10)
96 std_q10: nx_int, // sqrt(variance), Q10
97 t_stat_q10: nx_int, // t-statistic Q10
98 target_q10: nx_int // echo target for caller convenience
99}
100
101const NX_MC_RESULT_BYTES: nx_int = 56 // 7 fields * 8
102
103// Critical t-value for 95% CI at "large" sample count (n >= 30):
104// 1.96 in Q10 = 2007. Substrate keeps this as a constant. Stricter
105// confidence (99%, 99.9%) is a caller decision; expose constants if
106// needed.
107
108const NX_MC_T_95: nx_int = 2007 // 1.96 Q10
109const NX_MC_T_99: nx_int = 2640 // 2.58 Q10 (large-n approx)
110const NX_MC_Q10: nx_int = 1024
111
112// Minimum samples for the test to be statistically valid. Below
113// this the t-distribution is too wide to make a meaningful claim.
114const NX_MC_MIN_N: nx_int = 30
115
116// ===== Welford online stats helper ===============================
117//
118// Caller maintains (n, mean, M2) across samples; we update each.
119
120func _mc_welford_update(n_old: nx_int, mean_old: i64, M2_old: i64,
121 sample: i64,
122 out_n: *i64, out_mean: *i64, out_M2: *i64) -> nx_int {
123 let n_new: nx_int = n_old + 1
124 let delta: i64 = sample - mean_old
125 let mean_new: i64 = mean_old + delta / n_new
126 let delta2: i64 = sample - mean_new
127 let M2_new: i64 = M2_old + delta * delta2
128 out_n[0] = n_new
129 out_mean[0] = mean_new
130 out_M2[0] = M2_new
131 return 0
132}
133
134// ===== Main verify entrypoint =====================================
135//
136// sample_fn: caller's sampler; takes a *i64 prng_state, returns
137// one measurement (i64).
138// n_samples: how many to draw
139// target_q10: target value the mean should match
140// tol_q10: +/- tolerance for "match"
141// seed: PRNG seed
142// confidence: NX_MC_T_95 or NX_MC_T_99
143//
144// Returns *NxMonteCarloResult.
145
146func nx_monte_carlo_verify(
147 sample_fn: func(*i64) -> i64,
148 n_samples: nx_int,
149 target_q10: nx_int, tol_q10: nx_int,
150 seed: i64, confidence_t_q10: nx_int) -> *NxMonteCarloResult {
151
152 let r: *NxMonteCarloResult = sys_mmap(NX_MC_RESULT_BYTES) as *NxMonteCarloResult
153 r.verdict = NX_MC_ERR_BAD_PARAMS
154 r.n_samples = 0
155 r.target_q10 = target_q10
156
157 if n_samples <= 0 { return r }
158 if tol_q10 < 0 { return r }
159
160 let prng: *i64 = sys_mmap(8) as *i64
161 nx_prng_init(prng, seed)
162
163 // Welford accumulators.
164 let n_p: *i64 = sys_mmap(8) as *i64
165 let mean_p: *i64 = sys_mmap(8) as *i64
166 let M2_p: *i64 = sys_mmap(8) as *i64
167 n_p[0] = 0
168 mean_p[0] = 0
169 M2_p[0] = 0
170
171 var iter: nx_int = 0
172 var verdict: nx_int = NX_LOOP_RUNNING
173 let BUDGET: nx_int = n_samples
174 while verdict == NX_LOOP_RUNNING && iter < BUDGET {
175 let s: i64 = sample_fn(prng)
176 _mc_welford_update(n_p[0], mean_p[0], M2_p[0], s, n_p, mean_p, M2_p)
177 iter = iter + 1
178 }
179
180 let n: nx_int = n_p[0]
181 r.n_samples = n
182 r.mean_q10 = mean_p[0]
183
184 if n < NX_MC_MIN_N {
185 r.verdict = NX_MC_INSUFFICIENT_SAMPLES
186 return r
187 }
188
189 // Bessel-corrected sample variance: var = M2 / (n - 1).
190 let variance: i64 = M2_p[0] / (n - 1)
191 r.variance_q10 = variance
192
193 // Standard deviation: sqrt(variance). Variance is in Q20 if
194 // samples are in Q10; sqrt gives Q10 back via nx_isqrt_q10
195 // shape (the Q-scale-preserving sqrt).
196 let std: i64 = nx_isqrt_q10(variance + 1)
197 r.std_q10 = std
198
199 // t-statistic: t = (mean - target) * sqrt(n) / std.
200 // In Q10: numerator = (mean - target) * sqrt(n)_Q10 (we leave
201 // sqrt(n) integer since n is i64 and we want Q10 t).
202 let diff: i64 = mean_p[0] - target_q10
203 var abs_diff: i64 = diff
204 if abs_diff < 0 { abs_diff = 0 - abs_diff }
205
206 // sqrt(n) in Q10: nx_isqrt(n) gives integer sqrt; we scale to Q10
207 // by multiplying by Q10 first: sqrt(n * Q10^2) = sqrt(n) * Q10.
208 let sqrt_n_q10: i64 = nx_isqrt_q10(n * NX_MC_Q10)
209 // t = (abs_diff * sqrt_n_q10) / std. Both numerator factors Q10
210 // -> Q20; divide by std (Q10) -> Q10.
211 var t_stat: i64 = 0
212 if std > 0 {
213 t_stat = (abs_diff * sqrt_n_q10) / std
214 }
215 r.t_stat_q10 = t_stat
216
217 // Decision: if mean is within target +/- tol AND t < critical,
218 // accept. If outside tol OR t >= critical, reject.
219 let within_tol: nx_int = 0
220 if abs_diff <= tol_q10 {
221 // Within tol -- check if t-statistic confirms.
222 if t_stat < confidence_t_q10 {
223 r.verdict = NX_MC_TARGET_MET
224 }
225 if t_stat >= confidence_t_q10 {
226 r.verdict = NX_MC_OFF_TARGET
227 }
228 }
229 if abs_diff > tol_q10 {
230 r.verdict = NX_MC_OFF_TARGET
231 }
232 return r
233}
234
235// ===== Self-test ==================================================
236//
237// Sampler that returns the PRNG output mod 200 + 1000. Expected
238// mean is roughly 1099.5 Q10 (uniform [1000, 1199]). We verify the
239// MC verifier detects this distribution centred at ~1100 Q10.
240//
241// Closed-form invariants:
242// (a) Sampler with target = empirical mean (within tol)
243// -> TARGET_MET
244// (b) Sampler with target FAR from empirical mean (tol = 1)
245// -> OFF_TARGET
246// (c) Sample count < 30 -> INSUFFICIENT_SAMPLES
247// (d) Verdict gate
248
249func _mc_test_sampler(prng: *i64) -> i64 {
250 let r: i64 = nx_prng_range(prng, 200)
251 return r + 1000
252}
253
254func main() -> i64 {
255 // --- (a) Sampler around mean ~1099; target = 1100 with tol 50. ---
256 let r1: *NxMonteCarloResult = nx_monte_carlo_verify(
257 _mc_test_sampler, 500, 1100, 50, 0xdeadbeef, NX_MC_T_95)
258 if r1.verdict != NX_MC_TARGET_MET { return 10 }
259 if r1.n_samples != 500 { return 11 }
260 // Mean should be around 1099 +/- 5.
261 if r1.mean_q10 < 1080 { return 12 }
262 if r1.mean_q10 > 1120 { return 13 }
263
264 // --- (b) Target far off (mean is ~1099, target = 500). ---
265 let r2: *NxMonteCarloResult = nx_monte_carlo_verify(
266 _mc_test_sampler, 500, 500, 10, 0xfeedface, NX_MC_T_95)
267 if r2.verdict != NX_MC_OFF_TARGET { return 20 }
268
269 // --- (c) Tiny n -> INSUFFICIENT. ---
270 let r3: *NxMonteCarloResult = nx_monte_carlo_verify(
271 _mc_test_sampler, 10, 1100, 50, 0xc0ffee, NX_MC_T_95)
272 if r3.verdict != NX_MC_INSUFFICIENT_SAMPLES { return 30 }
273
274 // --- (d) Bad params (n_samples = 0). ---
275 let r4: *NxMonteCarloResult = nx_monte_carlo_verify(
276 _mc_test_sampler, 0, 1100, 50, 0xbadbad, NX_MC_T_95)
277 if r4.verdict != NX_MC_ERR_BAD_PARAMS { return 40 }
278
279 // --- (e) Verdict gate ---
280 var vi: nx_int = 0
281 while vi < NX_MC_N_VERDICTS {
282 if nx_mc_verdict_is_valid(vi) != 1 { return 50 + vi }
283 vi = vi + 1
284 }
285
286 return 0
287}