chronix_analytics/forecast/
util.rs

1//! Shared utilities for forecast models.
2
3use crate::forecast::error::ForecastError;
4
5/// Validate that all input values are finite (no NaN or Inf) and timestamps
6/// match values in length.
7#[inline]
8pub fn validate_input(timestamps: &[i64], values: &[f64]) -> Result<(), ForecastError> {
9    // Emptiness is checked here, in the one function every model calls,
10    // rather than in each model. It was checked in none of them: `SesModel`
11    // and `ArimaModel` survived an empty series by luck — their fitting loops
12    // are no-ops over it — and then answered `predict()` from an
13    // uninitialised level, which is a number with nothing behind it.
14    // `HoltWintersModel` indexed `values[0]` and **panicked**, which in an
15    // embedded database is the host process going down.
16    //
17    // A model fitted on no observations is not a model, so this is an error
18    // for all of them rather than a different accident in each.
19    if values.is_empty() {
20        return Err(ForecastError::InvalidInput(
21            "cannot fit a model to an empty series".to_string(),
22        ));
23    }
24    if timestamps.len() != values.len() {
25        return Err(ForecastError::InvalidInput(format!(
26            "timestamps length ({}) != values length ({})",
27            timestamps.len(),
28            values.len(),
29        )));
30    }
31    if let Some(pos) = values.iter().position(|v| !v.is_finite()) {
32        return Err(ForecastError::InvalidInput(format!(
33            "values[{pos}] is not finite ({})",
34            values[pos],
35        )));
36    }
37    Ok(())
38}
39
40/// Compute the median interval between sorted timestamps using O(n) selection.
41///
42/// Returns at least 1 to prevent division-by-zero downstream when all
43/// timestamps are identical (all deltas zero).
44pub fn median_interval(timestamps: &[i64]) -> i64 {
45    if timestamps.len() < 2 {
46        return 1;
47    }
48    let mut deltas: Vec<i64> = timestamps.windows(2).map(|w| w[1] - w[0]).collect();
49    let mid = deltas.len() / 2;
50    let median = *deltas.select_nth_unstable(mid).1;
51    median.max(1)
52}
53
54/// Standard-normal quantile function `Φ⁻¹(p)`, for `p ∈ (0, 1)`.
55///
56/// Acklam's rational approximation, whose relative error is below
57/// `1.15 × 10⁻⁹` over the whole open interval — several orders of magnitude
58/// finer than any prediction interval it is used to build. Returns `±∞` at the
59/// closed boundaries and `NaN` outside them, which is what the function is.
60#[must_use]
61pub fn normal_quantile(p: f64) -> f64 {
62    if !(0.0..=1.0).contains(&p) || p.is_nan() {
63        return f64::NAN;
64    }
65    if p == 0.0 {
66        return f64::NEG_INFINITY;
67    }
68    if p == 1.0 {
69        return f64::INFINITY;
70    }
71
72    const A: [f64; 6] = [
73        -3.969_683_028_665_376e1,
74        2.209_460_984_245_205e2,
75        -2.759_285_104_469_687e2,
76        1.383_577_518_672_69e2,
77        -3.066_479_806_614_716e1,
78        2.506_628_277_459_239e0,
79    ];
80    const B: [f64; 5] = [
81        -5.447_609_879_822_406e1,
82        1.615_858_368_580_409e2,
83        -1.556_989_798_598_866e2,
84        6.680_131_188_771_972e1,
85        -1.328_068_155_288_572e1,
86    ];
87    const C: [f64; 6] = [
88        -7.784_894_002_430_293e-3,
89        -3.223_964_580_411_365e-1,
90        -2.400_758_277_161_838e0,
91        -2.549_732_539_343_734e0,
92        4.374_664_141_464_968e0,
93        2.938_163_982_698_783e0,
94    ];
95    const D: [f64; 4] = [
96        7.784_695_709_041_462e-3,
97        3.224_671_290_700_398e-1,
98        2.445_134_137_142_996e0,
99        3.754_408_661_907_416e0,
100    ];
101    // Breakpoints between the tail and central rational fits.
102    const P_LOW: f64 = 0.024_25;
103    const P_HIGH: f64 = 1.0 - P_LOW;
104
105    if p < P_LOW {
106        let q = (-2.0 * p.ln()).sqrt();
107        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
108            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
109    } else if p <= P_HIGH {
110        let q = p - 0.5;
111        let r = q * q;
112        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
113            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
114    } else {
115        let q = (-2.0 * (1.0 - p).ln()).sqrt();
116        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
117            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
118    }
119}
120
121/// The `z` multiplier of a two-sided prediction interval at `level`.
122///
123/// `z_for_confidence(0.95) ≈ 1.96`. Returns `None` outside `(0, 1)`.
124#[must_use]
125pub fn z_for_confidence(level: f64) -> Option<f64> {
126    if !(level > 0.0 && level < 1.0) {
127        return None;
128    }
129    Some(normal_quantile(0.5 + level / 2.0))
130}
chronix_analytics/forecast/util.rs

chronix_analytics/forecast/
util.rs
