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