nx_sketch_mann_kendall.nx source

1// sketch_mann_kendall.nx -- streaming Mann-Kendall trend test. 2// 3// Non-parametric trend detection over the last W samples. No distribution 4// assumption -- works on any ordinal data including heavy-tailed, 5// skewed, or with outliers (vs Pearson correlation which assumes 6// linearity + normality). 7// 8// MANN-KENDALL STATISTIC: 9// For each pair (i, j) with i < j in the window: 10// S += sign(x_j - x_i) where sign in {-1, 0, +1} 11// S > 0 = upward trend 12// S < 0 = downward trend 13// |S| > threshold = significant 14// 15// THEORETICAL VARIANCE (no ties): 16// Var(S) = n(n-1)(2n+5)/18 17// Z = S / sqrt(Var(S)) — standardized, approximately N(0,1) for large n. 18// 19// FOR THE STREAMING VARIANT: maintain a sliding window of W samples. 20// On each new sample: compare against all W-1 previous samples in O(W). 21// Total state: O(W). S incrementally updated. 22// 23// USE CASES: 24// - hydrology: long-term streamflow trends 25// - climatology: warming trend detection 26// - SRE: monotonic latency trend 27// - finance: directional momentum (vs mean-reverting) 28// 29// Production tier; exact integer arithmetic over the window. 30 31// nx_safety_envelope: 32// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 33// sil_target: SIL1 34// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 35// verdict: NOT_YET_EVALUATED 36 37import "nx_syscalls.nx" 38import "nx_sketch_types.nx" 39 40const NX_MK_MIN_W: i64 = 3 41const NX_MK_MAX_W: i64 = 100000 42 43const NX_MK_TREND_NONE: i64 = 0 44const NX_MK_TREND_UP: i64 = 1 45const NX_MK_TREND_DOWN: i64 = 2 46 47struct MannKendall { 48 buffer: *i64, // circular ring of last W samples 49 window: i64, 50 count: i64, 51 head: i64, // next write position 52 s: i64, // running Mann-Kendall S over current window 53} 54 55// === construction ================================================= 56 57func nx_mk_alloc(window: i64) -> *MannKendall { 58 if window < NX_MK_MIN_W { return 0 as *MannKendall } 59 if window > NX_MK_MAX_W { return 0 as *MannKendall } 60 let raw: *u8 = sys_mmap(48) 61 let m: *MannKendall = raw as *MannKendall 62 m.buffer = sys_mmap(window * 8) as *i64 63 var i: i64 = 0 64 while i < window { 65 m.buffer[i] = 0 66 i = i + 1 67 } 68 m.window = window 69 m.count = 0 70 m.head = 0 71 m.s = 0 72 return m 73} 74 75// === sign ========================================================= 76 77func nx_mk_sign(x: i64) -> i64 { 78 if x > 0 { return 1 } 79 if x < 0 { return -1 } 80 return 0 81} 82 83// === recompute S over current window (used after eviction) ======== 84// 85// O(W²) -- called only when an eviction invalidates the running S 86// incrementally. v1 always recomputes; v2 can incrementally update 87// by subtracting all pairs involving the evicted sample. 88 89func nx_mk_recompute_s(m: *MannKendall) -> i64 { 90 var new_s: i64 = 0 91 // Buffer order: oldest at (head - count) mod window; newest at head-1. 92 var i: i64 = 0 93 while i < m.count { 94 var j: i64 = i + 1 95 while j < m.count { 96 // Physical positions in ring. 97 let pi: i64 = (m.head - m.count + i + m.window) % m.window 98 let pj: i64 = (m.head - m.count + j + m.window) % m.window 99 new_s = new_s + nx_mk_sign(m.buffer[pj] - m.buffer[pi]) 100 j = j + 1 101 } 102 i = i + 1 103 } 104 m.s = new_s 105 return 0 106} 107 108// === push ========================================================== 109// 110// O(W): compare new sample against all current samples (O(W)) -- add 111// each comparison to S. If buffer was full, evicted sample's 112// contributions must be subtracted (O(W) recompute). v1 just 113// recomputes on overflow for simplicity. 114 115func nx_mk_push(m: *MannKendall, value: i64) -> i64 { 116 if m.count >= m.window { 117 // Evict oldest and recompute (simpler than incremental subtract). 118 m.buffer[m.head] = value 119 m.head = (m.head + 1) % m.window 120 nx_mk_recompute_s(m) 121 return 0 122 } 123 // Under-cap: append, update S incrementally. 124 // For each EXISTING sample (oldest..newest), S += sign(value - sample). 125 var i: i64 = 0 126 while i < m.count { 127 let pos: i64 = (m.head - m.count + i + m.window) % m.window 128 m.s = m.s + nx_mk_sign(value - m.buffer[pos]) 129 i = i + 1 130 } 131 m.buffer[m.head] = value 132 m.head = (m.head + 1) % m.window 133 m.count = m.count + 1 134 return 0 135} 136 137// === variance bound ============================================== 138// 139// Var(S) = n(n-1)(2n+5)/18 under null (no trend, no ties). 140 141func nx_mk_var_s(m: *MannKendall) -> i64 { 142 let n: i64 = m.count 143 if n < 2 { return 0 } 144 return (n * (n - 1) * (2 * n + 5)) / 18 145} 146 147// === verdict ====================================================== 148// 149// Use a critical-value-style cutoff: |S| > 1.96 * sqrt(Var(S)) for 150// 95% confidence (two-sided). We approximate this as |S|² > 4 * Var(S). 151// Returns sealed enum: NONE / UP / DOWN. 152 153func nx_mk_verdict(m: *MannKendall) -> i64 { 154 if m.count < 3 { return NX_MK_TREND_NONE } 155 let vs: i64 = nx_mk_var_s(m) 156 if vs == 0 { return NX_MK_TREND_NONE } 157 let s_sq: i64 = m.s * m.s 158 // |S|² > 4 * Var(S) approximates |Z| > 2 (~95% confidence two-sided). 159 if s_sq <= 4 * vs { return NX_MK_TREND_NONE } 160 if m.s > 0 { return NX_MK_TREND_UP } 161 return NX_MK_TREND_DOWN 162} 163 164// === introspection ================================================ 165 166func nx_mk_s(m: *MannKendall) -> i64 { return m.s } 167func nx_mk_count(m: *MannKendall) -> i64 { return m.count } 168 169func nx_mk_query(m: *MannKendall) -> *ApproxI64 { 170 let v: i64 = nx_mk_verdict(m) 171 return nx_approx_new(v, NX_ENV_ABS, 0, 950000000, 172 NX_MATURITY_PRODUCTION, 173 NX_ADV_HONEST) 174} 175 176func nx_mk_memory_bytes(m: *MannKendall) -> i64 { 177 return 48 + m.window * 8 178}