nx_sketch_moments.nx source
↩ module page · 208 lines · 6669 B
1// sketch_moments.nx -- streaming higher moments (skewness + kurtosis).
2//
3// Extends sketch_stream_stats with 3rd and 4th central moments via a
4// naive cumulative-power approach:
5// sum_x = Σ x_i
6// sum_x2 = Σ x_i²
7// sum_x3 = Σ x_i³
8// sum_x4 = Σ x_i⁴
9//
10// From these, central moments via the algebraic identities:
11// mean = sum_x / n
12// var = sum_x2/n - mean²
13// m3 = sum_x3/n - 3·mean·var - mean³ (third central)
14// m4 = sum_x4/n - 4·mean·(sum_x3/n) + 6·mean²·(sum_x2/n) - 3·mean⁴
15// skewness = m3 / σ³ (Fisher-Pearson)
16// kurtosis = m4 / σ⁴ - 3 (excess kurtosis; normal distribution = 0)
17//
18// OVERFLOW BUDGET (CRITICAL):
19// sum_x4 grows like n · |x|⁴. For |x| up to 2^16 and n up to 2^30,
20// sum_x4 < 2^30 · 2^64 → overflows i64 (which caps at 2^63).
21// Tight bound: n · |x|⁴ < 2^62. E.g. |x| <= 2^11 (~2K) and n <= 2^18 (~260K)
22// are safe. Caller responsibility. nx_mom_safe_p flags unsafe values.
23//
24// FOR PROPER OBSERVABILITY ENGINEERING: use this with normalized /
25// quantized inputs (e.g. milliseconds with values < 2K, sample counts
26// < 100K). For broader ranges, use sketch_reservoir + caller-side
27// computation on the sample.
28//
29// LOSSLESS-LANGUAGE DISCIPLINE: skewness/kurtosis envelope is
30// NX_ENV_REL_STDDEV with param_a ~ 1/sqrt(n) (standard error of
31// higher-moment estimators). Caller treats these as bounded-error
32// statistics over the integer-exact accumulators.
33
34// nx_safety_envelope:
35// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
36// sil_target: SIL1
37// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
38// verdict: NOT_YET_EVALUATED
39
40import "nx_syscalls.nx"
41import "nx_sketch_stream_stats.nx"
42import "nx_sketch_types.nx"
43import "nx_vecmath.nx"
44
45const NX_MOM_X_LIMIT: i64 = 2048 // |x| < 2^11 for safe x^4
46
47struct Moments {
48 count: i64,
49 sum_x: i64,
50 sum_x2: i64,
51 sum_x3: i64,
52 sum_x4: i64,
53}
54
55// === construction =================================================
56
57func nx_mom_alloc() -> *Moments {
58 let raw: *u8 = sys_mmap(40)
59 let m: *Moments = raw as *Moments
60 m.count = 0
61 m.sum_x = 0
62 m.sum_x2 = 0
63 m.sum_x3 = 0
64 m.sum_x4 = 0
65 return m
66}
67
68// === overflow safety check ========================================
69
70func nx_mom_safe_p(value: i64) -> i64 {
71 var v: i64 = value
72 if v < 0 { v = -v }
73 if v >= NX_MOM_X_LIMIT { return 0 }
74 return 1
75}
76
77// === add ==========================================================
78
79func nx_mom_add(m: *Moments, value: i64) -> i64 {
80 if nx_mom_safe_p(value) == 0 { return -1 }
81 let v2: i64 = value * value
82 let v3: i64 = v2 * value
83 let v4: i64 = v2 * v2
84 m.count = m.count + 1
85 m.sum_x = m.sum_x + value
86 m.sum_x2 = m.sum_x2 + v2
87 m.sum_x3 = m.sum_x3 + v3
88 m.sum_x4 = m.sum_x4 + v4
89 return 0
90}
91
92// === queries ======================================================
93
94func nx_mom_mean(m: *Moments) -> i64 {
95 if m.count == 0 { return 0 }
96 return m.sum_x / m.count
97}
98
99func nx_mom_variance(m: *Moments) -> i64 {
100 if m.count == 0 { return 0 }
101 let mn: i64 = nx_mom_mean(m)
102 let e_sq: i64 = m.sum_x2 / m.count
103 let mn_sq: i64 = mn * mn
104 if e_sq < mn_sq { return 0 }
105 return e_sq - mn_sq
106}
107
108// Third central moment, in (value units)^3.
109func nx_mom_m3(m: *Moments) -> i64 {
110 if m.count == 0 { return 0 }
111 let mn: i64 = nx_mom_mean(m)
112 let e1: i64 = m.sum_x3 / m.count
113 let e2: i64 = (3 * mn * m.sum_x2) / m.count
114 let mn3: i64 = mn * mn * mn
115 return e1 - e2 + 2 * mn3
116}
117
118// Fourth central moment, in (value units)^4.
119func nx_mom_m4(m: *Moments) -> i64 {
120 if m.count == 0 { return 0 }
121 let mn: i64 = nx_mom_mean(m)
122 let e1: i64 = m.sum_x4 / m.count
123 let e2: i64 = (4 * mn * m.sum_x3) / m.count
124 let e3: i64 = (6 * mn * mn * m.sum_x2) / m.count
125 let mn4: i64 = mn * mn * mn * mn
126 return e1 - e2 + e3 - 3 * mn4
127}
128
129// === skewness + kurtosis (in PPM) ================================
130//
131// skewness_ppm = m3 / sigma^3 * 1_000_000
132// sigma^3 = isqrt(variance)^3 (approximate; integer math)
133// kurtosis_ppm = m4 / sigma^4 * 1_000_000 - 3_000_000 (excess form)
134
135func nx_mom_isqrt(x: i64) -> i64 { return vm_isqrt(x) }
136
137func nx_mom_skewness_ppm(m: *Moments) -> i64 {
138 let var_val: i64 = nx_mom_variance(m)
139 if var_val == 0 { return 0 }
140 let sigma: i64 = nx_mom_isqrt(var_val)
141 if sigma == 0 { return 0 }
142 let sigma3: i64 = sigma * sigma * sigma
143 if sigma3 == 0 { return 0 }
144 let m3: i64 = nx_mom_m3(m)
145 return (m3 * 1000000) / sigma3
146}
147
148func nx_mom_kurtosis_ppm(m: *Moments) -> i64 {
149 let var_val: i64 = nx_mom_variance(m)
150 if var_val == 0 { return 0 }
151 let sigma: i64 = nx_mom_isqrt(var_val)
152 let sigma4: i64 = sigma * sigma * sigma * sigma
153 if sigma4 == 0 { return 0 }
154 let m4: i64 = nx_mom_m4(m)
155 return (m4 * 1000000) / sigma4 - 3000000
156}
157
158// === typed envelope ===============================================
159//
160// Standard error of the skewness estimator under normal distribution:
161// se(skew) ≈ sqrt(6/n)
162// For n=1000: se=0.0775, in ppm: 77_460.
163//
164// We declare param_a = 1/sqrt(n) ppb (rough bound, ignoring distribution-
165// specific factors).
166
167func nx_mom_stderr_ppb(count: i64) -> i64 {
168 if count < 1 { return 1000000000 }
169 let isq: i64 = nx_mom_isqrt(count)
170 if isq == 0 { return 1000000000 }
171 return 1000000000 / isq
172}
173
174func nx_mom_query_skewness(m: *Moments) -> *ApproxI64 {
175 let s: i64 = nx_mom_skewness_ppm(m)
176 return nx_approx_new(s, NX_ENV_REL_STDDEV,
177 nx_mom_stderr_ppb(m.count),
178 682700000,
179 NX_MATURITY_REFERENCE_IMPL,
180 NX_ADV_HONEST)
181}
182
183func nx_mom_query_kurtosis(m: *Moments) -> *ApproxI64 {
184 let k: i64 = nx_mom_kurtosis_ppm(m)
185 return nx_approx_new(k, NX_ENV_REL_STDDEV,
186 nx_mom_stderr_ppb(m.count),
187 682700000,
188 NX_MATURITY_REFERENCE_IMPL,
189 NX_ADV_HONEST)
190}
191
192// === merge (Chan 1979 generalization) =============================
193//
194// Combine two independent moment accumulators. All four sums add exactly.
195
196func nx_mom_merge(a: *Moments, b: *Moments) -> *Moments {
197 let out: *Moments = nx_mom_alloc()
198 out.count = a.count + b.count
199 out.sum_x = a.sum_x + b.sum_x
200 out.sum_x2 = a.sum_x2 + b.sum_x2
201 out.sum_x3 = a.sum_x3 + b.sum_x3
202 out.sum_x4 = a.sum_x4 + b.sum_x4
203 return out
204}
205
206func nx_mom_memory_bytes(m: *Moments) -> i64 {
207 return 40
208}