nx_sketch_moments.nx source
↩ module page · 222 lines · 6963 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"
43
44const NX_MOM_X_LIMIT: i64 = 2048 // |x| < 2^11 for safe x^4
45
46struct Moments {
47 count: i64,
48 sum_x: i64,
49 sum_x2: i64,
50 sum_x3: i64,
51 sum_x4: i64,
52}
53
54// === construction =================================================
55
56func nx_mom_alloc() -> *Moments {
57 let raw: *u8 = sys_mmap(40)
58 let m: *Moments = raw as *Moments
59 m.count = 0
60 m.sum_x = 0
61 m.sum_x2 = 0
62 m.sum_x3 = 0
63 m.sum_x4 = 0
64 return m
65}
66
67// === overflow safety check ========================================
68
69func nx_mom_safe_p(value: i64) -> i64 {
70 var v: i64 = value
71 if v < 0 { v = -v }
72 if v >= NX_MOM_X_LIMIT { return 0 }
73 return 1
74}
75
76// === add ==========================================================
77
78func nx_mom_add(m: *Moments, value: i64) -> i64 {
79 if nx_mom_safe_p(value) == 0 { return -1 }
80 let v2: i64 = value * value
81 let v3: i64 = v2 * value
82 let v4: i64 = v2 * v2
83 m.count = m.count + 1
84 m.sum_x = m.sum_x + value
85 m.sum_x2 = m.sum_x2 + v2
86 m.sum_x3 = m.sum_x3 + v3
87 m.sum_x4 = m.sum_x4 + v4
88 return 0
89}
90
91// === queries ======================================================
92
93func nx_mom_mean(m: *Moments) -> i64 {
94 if m.count == 0 { return 0 }
95 return m.sum_x / m.count
96}
97
98func nx_mom_variance(m: *Moments) -> i64 {
99 if m.count == 0 { return 0 }
100 let mn: i64 = nx_mom_mean(m)
101 let e_sq: i64 = m.sum_x2 / m.count
102 let mn_sq: i64 = mn * mn
103 if e_sq < mn_sq { return 0 }
104 return e_sq - mn_sq
105}
106
107// Third central moment, in (value units)^3.
108func nx_mom_m3(m: *Moments) -> i64 {
109 if m.count == 0 { return 0 }
110 let mn: i64 = nx_mom_mean(m)
111 let e1: i64 = m.sum_x3 / m.count
112 let e2: i64 = (3 * mn * m.sum_x2) / m.count
113 let mn3: i64 = mn * mn * mn
114 return e1 - e2 + 2 * mn3
115}
116
117// Fourth central moment, in (value units)^4.
118func nx_mom_m4(m: *Moments) -> i64 {
119 if m.count == 0 { return 0 }
120 let mn: i64 = nx_mom_mean(m)
121 let e1: i64 = m.sum_x4 / m.count
122 let e2: i64 = (4 * mn * m.sum_x3) / m.count
123 let e3: i64 = (6 * mn * mn * m.sum_x2) / m.count
124 let mn4: i64 = mn * mn * mn * mn
125 return e1 - e2 + e3 - 3 * mn4
126}
127
128// === skewness + kurtosis (in PPM) ================================
129//
130// skewness_ppm = m3 / sigma^3 * 1_000_000
131// sigma^3 = isqrt(variance)^3 (approximate; integer math)
132// kurtosis_ppm = m4 / sigma^4 * 1_000_000 - 3_000_000 (excess form)
133
134func nx_mom_isqrt(x: i64) -> i64 {
135 if x < 0 { return 0 }
136 if x == 0 { return 0 }
137 if x < 4 { return 1 }
138 var g: i64 = (x >> 1) + 1
139 var iter: i64 = 0
140 while iter < 64 {
141 let next_g: i64 = (g + x / g) / 2
142 if next_g >= g { iter = 64 }
143 if next_g < g {
144 g = next_g
145 iter = iter + 1
146 }
147 }
148 return g
149}
150
151func nx_mom_skewness_ppm(m: *Moments) -> i64 {
152 let var_val: i64 = nx_mom_variance(m)
153 if var_val == 0 { return 0 }
154 let sigma: i64 = nx_mom_isqrt(var_val)
155 if sigma == 0 { return 0 }
156 let sigma3: i64 = sigma * sigma * sigma
157 if sigma3 == 0 { return 0 }
158 let m3: i64 = nx_mom_m3(m)
159 return (m3 * 1000000) / sigma3
160}
161
162func nx_mom_kurtosis_ppm(m: *Moments) -> i64 {
163 let var_val: i64 = nx_mom_variance(m)
164 if var_val == 0 { return 0 }
165 let sigma: i64 = nx_mom_isqrt(var_val)
166 let sigma4: i64 = sigma * sigma * sigma * sigma
167 if sigma4 == 0 { return 0 }
168 let m4: i64 = nx_mom_m4(m)
169 return (m4 * 1000000) / sigma4 - 3000000
170}
171
172// === typed envelope ===============================================
173//
174// Standard error of the skewness estimator under normal distribution:
175// se(skew) ≈ sqrt(6/n)
176// For n=1000: se=0.0775, in ppm: 77_460.
177//
178// We declare param_a = 1/sqrt(n) ppb (rough bound, ignoring distribution-
179// specific factors).
180
181func nx_mom_stderr_ppb(count: i64) -> i64 {
182 if count < 1 { return 1000000000 }
183 let isq: i64 = nx_mom_isqrt(count)
184 if isq == 0 { return 1000000000 }
185 return 1000000000 / isq
186}
187
188func nx_mom_query_skewness(m: *Moments) -> *ApproxI64 {
189 let s: i64 = nx_mom_skewness_ppm(m)
190 return nx_approx_new(s, NX_ENV_REL_STDDEV,
191 nx_mom_stderr_ppb(m.count),
192 682700000,
193 NX_MATURITY_REFERENCE_IMPL,
194 NX_ADV_HONEST)
195}
196
197func nx_mom_query_kurtosis(m: *Moments) -> *ApproxI64 {
198 let k: i64 = nx_mom_kurtosis_ppm(m)
199 return nx_approx_new(k, NX_ENV_REL_STDDEV,
200 nx_mom_stderr_ppb(m.count),
201 682700000,
202 NX_MATURITY_REFERENCE_IMPL,
203 NX_ADV_HONEST)
204}
205
206// === merge (Chan 1979 generalization) =============================
207//
208// Combine two independent moment accumulators. All four sums add exactly.
209
210func nx_mom_merge(a: *Moments, b: *Moments) -> *Moments {
211 let out: *Moments = nx_mom_alloc()
212 out.count = a.count + b.count
213 out.sum_x = a.sum_x + b.sum_x
214 out.sum_x2 = a.sum_x2 + b.sum_x2
215 out.sum_x3 = a.sum_x3 + b.sum_x3
216 out.sum_x4 = a.sum_x4 + b.sum_x4
217 return out
218}
219
220func nx_mom_memory_bytes(m: *Moments) -> i64 {
221 return 40
222}