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