quantwave_core/regimes/
ecld.rs

1//! Symmetric lambda (exponential power / generalized normal) emission density.
2//!
3//! ldhmm uses `ecld.pdf` via `gnorm::dgnorm` with `beta = 2/λ`. At λ=1 this is Gaussian.
4//!
5//! Source: references/ldhmm/ssrn-2979516.pdf §2; ldhmm ecld.* functions.
6
7use statrs::distribution::{ContinuousCDF, Normal};
8use statrs::function::gamma::{gamma, gamma_lr};
9
10/// PDF of the symmetric lambda distribution (ecld / generalized normal).
11///
12/// `pdf(x) = (β / (2σ Γ(1/β))) exp(-|(x-μ)/σ|^β)` where `β = 2/λ`.
13/// λ=1 delegates to the standard normal density for bit-identical Gaussian parity.
14#[inline]
15pub fn ecld_pdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
16    if (lambda - 1.0).abs() < 1e-12 {
17        return gaussian_pdf(x, mu, sigma);
18    }
19    let beta = 2.0 / lambda;
20    let z = ((x - mu) / sigma).abs();
21    let norm = beta / (2.0 * sigma * gamma(1.0 / beta));
22    norm * (-z.powf(beta)).exp()
23}
24
25#[inline]
26pub fn ecld_log_pdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
27    ecld_pdf(x, mu, sigma, lambda).ln()
28}
29
30/// CDF of the symmetric lambda distribution (ldhmm `ecld.cdf` / gnorm).
31#[inline]
32pub fn ecld_cdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
33    if (lambda - 1.0).abs() < 1e-12 {
34        return Normal::new(mu, sigma).map(|n| n.cdf(x)).unwrap_or(0.5);
35    }
36    let beta = 2.0 / lambda;
37    let u = (x - mu) / sigma;
38    let p = gamma_lr(1.0 / beta, u.abs().powf(beta));
39    if u >= 0.0 {
40        0.5 + 0.5 * p
41    } else {
42        0.5 - 0.5 * p
43    }
44}
45
46/// Emission variance of the symmetric lambda distribution: `σ² Γ(3/β) / Γ(1/β)`, β=2/λ.
47#[inline]
48pub fn ecld_variance(sigma: f64, lambda: f64) -> f64 {
49    if (lambda - 1.0).abs() < 1e-12 {
50        return sigma * sigma;
51    }
52    let beta = 2.0 / lambda;
53    let s2 = sigma * sigma;
54    s2 * gamma(3.0 / beta) / gamma(1.0 / beta)
55}
56
57#[inline]
58fn gaussian_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
59    let variance = sigma * sigma;
60    let denom = (2.0 * std::f64::consts::PI * variance).sqrt();
61    let exponent = -((x - mu).powi(2)) / (2.0 * variance);
62    exponent.exp() / denom
63}
64
65/// Natural-parameter style mapping helpers (ldhmm `n2w` / `w2n` stability).
66///
67/// Work parameters: `(μ, log σ, log λ)` ↔ natural-ish `(μ, σ, λ)`.
68pub fn work_to_natural(mu: f64, log_sigma: f64, log_lambda: f64) -> (f64, f64, f64) {
69    (mu, log_sigma.exp(), log_lambda.exp())
70}
71
72pub fn natural_to_work(mu: f64, sigma: f64, lambda: f64) -> (f64, f64, f64) {
73    (mu, sigma.ln(), lambda.ln())
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79    use approx::assert_relative_eq;
80
81    #[test]
82    fn lambda_one_matches_gaussian() {
83        let x = 0.005;
84        let mu = 0.008;
85        let sigma = 0.018;
86        let g = gaussian_pdf(x, mu, sigma);
87        let l = ecld_pdf(x, mu, sigma, 1.0);
88        assert_relative_eq!(l, g, epsilon = 1e-15);
89    }
90
91    #[test]
92    fn leptokurtic_lambda_increases_center_density() {
93        let x = 0.008; // at mean
94        let mu = 0.008;
95        let sigma = 0.018;
96        let g = ecld_pdf(x, mu, sigma, 1.0);
97        let l = ecld_pdf(x, mu, sigma, 1.3);
98        assert!(l > g);
99    }
100}
quantwave_core/regimes/ecld.rs

quantwave_core/regimes/
ecld.rs
