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)
35            .expect("valid gaussian cdf params")
36            .cdf(x);
37    }
38    let beta = 2.0 / lambda;
39    let u = (x - mu) / sigma;
40    let p = gamma_lr(1.0 / beta, u.abs().powf(beta));
41    if u >= 0.0 {
42        0.5 + 0.5 * p
43    } else {
44        0.5 - 0.5 * p
45    }
46}
47
48/// Emission variance of the symmetric lambda distribution: `σ² Γ(3/β) / Γ(1/β)`, β=2/λ.
49#[inline]
50pub fn ecld_variance(sigma: f64, lambda: f64) -> f64 {
51    if (lambda - 1.0).abs() < 1e-12 {
52        return sigma * sigma;
53    }
54    let beta = 2.0 / lambda;
55    let s2 = sigma * sigma;
56    s2 * gamma(3.0 / beta) / gamma(1.0 / beta)
57}
58
59#[inline]
60fn gaussian_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
61    let variance = sigma * sigma;
62    let denom = (2.0 * std::f64::consts::PI * variance).sqrt();
63    let exponent = -((x - mu).powi(2)) / (2.0 * variance);
64    exponent.exp() / denom
65}
66
67/// Natural-parameter style mapping helpers (ldhmm `n2w` / `w2n` stability).
68///
69/// Work parameters: `(μ, log σ, log λ)` ↔ natural-ish `(μ, σ, λ)`.
70pub fn work_to_natural(mu: f64, log_sigma: f64, log_lambda: f64) -> (f64, f64, f64) {
71    (mu, log_sigma.exp(), log_lambda.exp())
72}
73
74pub fn natural_to_work(mu: f64, sigma: f64, lambda: f64) -> (f64, f64, f64) {
75    (mu, sigma.ln(), lambda.ln())
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81    use approx::assert_relative_eq;
82
83    #[test]
84    fn lambda_one_matches_gaussian() {
85        let x = 0.005;
86        let mu = 0.008;
87        let sigma = 0.018;
88        let g = gaussian_pdf(x, mu, sigma);
89        let l = ecld_pdf(x, mu, sigma, 1.0);
90        assert_relative_eq!(l, g, epsilon = 1e-15);
91    }
92
93    #[test]
94    fn leptokurtic_lambda_increases_center_density() {
95        let x = 0.008; // at mean
96        let mu = 0.008;
97        let sigma = 0.018;
98        let g = ecld_pdf(x, mu, sigma, 1.0);
99        let l = ecld_pdf(x, mu, sigma, 1.3);
100        assert!(l > g);
101    }
102}
quantwave_core/regimes/ecld.rs

quantwave_core/regimes/
ecld.rs
