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 31import "syscalls.nx" 32import "sketch_types.nx" 33 34const NX_MK_MIN_W: i64 = 3 35const NX_MK_MAX_W: i64 = 100000 36 37const NX_MK_TREND_NONE: i64 = 0 38const NX_MK_TREND_UP: i64 = 1 39const NX_MK_TREND_DOWN: i64 = 2 40 41struct MannKendall { 42 buffer: *i64, // circular ring of last W samples 43 window: i64, 44 count: i64, 45 head: i64, // next write position 46 s: i64, // running Mann-Kendall S over current window 47} 48 49// === construction ================================================= 50 51func nx_mk_alloc(window: i64) -> *MannKendall { 52 if window < NX_MK_MIN_W { return 0 as *MannKendall } 53 if window > NX_MK_MAX_W { return 0 as *MannKendall } 54 let raw: *u8 = sys_mmap(48) 55 let m: *MannKendall = raw as *MannKendall 56 m.buffer = sys_mmap(window * 8) as *i64 57 var i: i64 = 0 58 while i < window { 59 m.buffer[i] = 0 60 i = i + 1 61 } 62 m.window = window 63 m.count = 0 64 m.head = 0 65 m.s = 0 66 return m 67} 68 69// === sign ========================================================= 70 71func nx_mk_sign(x: i64) -> i64 { 72 if x > 0 { return 1 } 73 if x < 0 { return -1 } 74 return 0 75} 76 77// === recompute S over current window (used after eviction) ======== 78// 79// O(W²) -- called only when an eviction invalidates the running S 80// incrementally. v1 always recomputes; v2 can incrementally update 81// by subtracting all pairs involving the evicted sample. 82 83func nx_mk_recompute_s(m: *MannKendall) -> i64 { 84 var new_s: i64 = 0 85 // Buffer order: oldest at (head - count) mod window; newest at head-1. 86 var i: i64 = 0 87 while i < m.count { 88 var j: i64 = i + 1 89 while j < m.count { 90 // Physical positions in ring. 91 let pi: i64 = (m.head - m.count + i + m.window) % m.window 92 let pj: i64 = (m.head - m.count + j + m.window) % m.window 93 new_s = new_s + nx_mk_sign(m.buffer[pj] - m.buffer[pi]) 94 j = j + 1 95 } 96 i = i + 1 97 } 98 m.s = new_s 99 return 0 100} 101 102// === push ========================================================== 103// 104// O(W): compare new sample against all current samples (O(W)) -- add 105// each comparison to S. If buffer was full, evicted sample's 106// contributions must be subtracted (O(W) recompute). v1 just 107// recomputes on overflow for simplicity. 108 109func nx_mk_push(m: *MannKendall, value: i64) -> i64 { 110 if m.count >= m.window { 111 // Evict oldest and recompute (simpler than incremental subtract). 112 m.buffer[m.head] = value 113 m.head = (m.head + 1) % m.window 114 nx_mk_recompute_s(m) 115 return 0 116 } 117 // Under-cap: append, update S incrementally. 118 // For each EXISTING sample (oldest..newest), S += sign(value - sample). 119 var i: i64 = 0 120 while i < m.count { 121 let pos: i64 = (m.head - m.count + i + m.window) % m.window 122 m.s = m.s + nx_mk_sign(value - m.buffer[pos]) 123 i = i + 1 124 } 125 m.buffer[m.head] = value 126 m.head = (m.head + 1) % m.window 127 m.count = m.count + 1 128 return 0 129} 130 131// === variance bound ============================================== 132// 133// Var(S) = n(n-1)(2n+5)/18 under null (no trend, no ties). 134 135func nx_mk_var_s(m: *MannKendall) -> i64 { 136 let n: i64 = m.count 137 if n < 2 { return 0 } 138 return (n * (n - 1) * (2 * n + 5)) / 18 139} 140 141// === verdict ====================================================== 142// 143// Use a critical-value-style cutoff: |S| > 1.96 * sqrt(Var(S)) for 144// 95% confidence (two-sided). We approximate this as |S|² > 4 * Var(S). 145// Returns sealed enum: NONE / UP / DOWN. 146 147func nx_mk_verdict(m: *MannKendall) -> i64 { 148 if m.count < 3 { return NX_MK_TREND_NONE } 149 let vs: i64 = nx_mk_var_s(m) 150 if vs == 0 { return NX_MK_TREND_NONE } 151 let s_sq: i64 = m.s * m.s 152 // |S|² > 4 * Var(S) approximates |Z| > 2 (~95% confidence two-sided). 153 if s_sq <= 4 * vs { return NX_MK_TREND_NONE } 154 if m.s > 0 { return NX_MK_TREND_UP } 155 return NX_MK_TREND_DOWN 156} 157 158// === introspection ================================================ 159 160func nx_mk_s(m: *MannKendall) -> i64 { return m.s } 161func nx_mk_count(m: *MannKendall) -> i64 { return m.count } 162 163func nx_mk_query(m: *MannKendall) -> *ApproxI64 { 164 let v: i64 = nx_mk_verdict(m) 165 return nx_approx_new(v, NX_ENV_ABS, 0, 950000000, 166 NX_MATURITY_PRODUCTION, 167 NX_ADV_HONEST) 168} 169 170func nx_mk_memory_bytes(m: *MannKendall) -> i64 { 171 return 48 + m.window * 8 172}