1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
//! Skew-Normal leaf — moment-matched Gaussian output.
//!
//! Skew-Normal has location `ξ`, scale `ω`, and skewness `α`, with:
//!
//! ```text
//! E[Y] = ξ + ω · δ · √(2/π) where δ = α / √(1 + α²)
//! Var[Y] = ω² · (1 − 2δ² / π)
//! γ₁ = (4 − π) / 2 · (δ · √(2/π))³ / (1 − 2δ²/π)^(3/2)
//! ```
//!
//! Method-of-moments fit: solve for `δ` from sample skewness `γ₁`,
//! back out `α = δ / √(1 − δ²)`. When `|γ₁| < 0.05` we treat as
//! symmetric (Gaussian fallback). When `|γ₁| > 0.99` (near degenerate)
//! we clamp δ to ±0.99.
//!
//! Useful for asymmetric continuous data where YJ/log doesn't fully
//! symmetrize — right-skewed retail aggregates that aren't cleanly
//! log-normal, or left-skewed decays.
//!
//! Output is `Gaussian(E[Y], √(Var[Y] · h))`. The true skew shape is
//! not preserved end-to-end (TypedMixture, planned α-24 Level 2, would
//! fix).
use crate::models::laplace::dist::Gaussian;
use crate::models::laplace::leaf::Leaf;
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct SkewNormalLeaf {
alpha: f64,
/// Location EMA.
xi_ema: f64,
initialized: bool,
/// Welford for M2 and M3 (2nd and 3rd central moments).
n: usize,
m1: f64,
m2: f64,
m3: f64,
}
impl SkewNormalLeaf {
pub fn new(alpha: f64) -> Self {
Self {
alpha: alpha.clamp(1e-3, 1.0 - 1e-3),
xi_ema: 0.0,
initialized: false,
n: 0,
m1: 0.0,
m2: 0.0,
m3: 0.0,
}
}
fn variance(&self) -> f64 {
if self.n < 2 {
return 1.0;
}
(self.m2 / (self.n as f64 - 1.0)).max(1e-9)
}
/// Sample skewness γ₁ = (M3/n) / (M2/n)^(3/2).
fn skewness(&self) -> f64 {
if self.n < 30 {
return 0.0;
}
let n = self.n as f64;
let sd_biased = (self.m2 / n).sqrt();
if sd_biased < 1e-9 {
return 0.0;
}
(self.m3 / n) / sd_biased.powi(3)
}
}
impl Leaf for SkewNormalLeaf {
fn name(&self) -> &'static str {
"skew_normal"
}
fn predict(&self, horizon: usize) -> Vec<Gaussian> {
let mean = self.xi_ema;
let var_base = self.variance();
// Skewness-adjusted variance. When sample skewness is small,
// fall through to plain σ². When large, the skew-normal
// variance factor (1 − 2δ²/π) shrinks it.
let g1 = self.skewness();
let variance = if g1.abs() < 0.05 {
var_base
} else {
// Recover δ from γ₁ via a first-order approximation:
// γ₁ ≈ (4 − π)/2 · c³ / (1 − c²)^(3/2), c = δ √(2/π).
// For MVP we invert numerically with a fixed heuristic.
let approx_c = (g1 * 0.5).clamp(-0.85, 0.85);
let approx_delta_sq = approx_c * approx_c * std::f64::consts::PI / 2.0;
let variance_factor =
(1.0 - 2.0 * approx_delta_sq / std::f64::consts::PI).clamp(0.1, 1.0);
var_base / variance_factor.max(0.1)
};
let variance = variance.max(1e-9);
(1..=horizon)
.map(|h| Gaussian::new(mean, (variance * h as f64).sqrt()))
.collect()
}
fn observe(&mut self, y: f64) {
if !self.initialized {
self.xi_ema = y;
self.initialized = true;
} else {
self.xi_ema = self.alpha * y + (1.0 - self.alpha) * self.xi_ema;
}
// Welford for M2 and M3 (single-pass, approximate M3 update).
self.n += 1;
let n = self.n as f64;
let delta = y - self.m1;
let delta_n = delta / n;
let term1 = delta * delta_n * (n - 1.0);
self.m1 += delta_n;
self.m3 += term1 * delta_n * (n - 2.0) - 3.0 * delta_n * self.m2;
self.m2 += term1;
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn locks_on_to_symmetric_data() {
let mut leaf = SkewNormalLeaf::new(0.05);
for i in 0..500 {
let y = 5.0 + ((i as f64 * 0.1).sin() * 1.5);
leaf.observe(y);
}
let preds = leaf.predict(3);
assert!(
(preds[0].mean - 5.0).abs() < 2.0,
"expected ~5±2, got {}",
preds[0].mean
);
assert!(preds[0].std > 0.0);
}
#[test]
fn detects_right_skew() {
let mut leaf = SkewNormalLeaf::new(0.05);
// Right-skewed: mostly small, occasional big spike.
for i in 0..300 {
let y = if i % 20 == 0 { 10.0 } else { 1.0 };
leaf.observe(y);
}
assert!(leaf.skewness() > 0.5, "expected positive skew");
let preds = leaf.predict(3);
assert!(preds[0].mean.is_finite());
assert!(preds[0].std > 0.0);
}
#[test]
fn cold_start_produces_finite_predictions() {
let mut leaf = SkewNormalLeaf::new(0.1);
leaf.observe(1.0);
leaf.observe(3.0);
let preds = leaf.predict(3);
for p in preds {
assert!(p.mean.is_finite());
assert!(p.std.is_finite() && p.std > 0.0);
}
}
}