use statrs::distribution::{ContinuousCDF, Normal};
use statrs::function::gamma::{gamma, gamma_lr};
#[inline]
pub fn ecld_pdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
if (lambda - 1.0).abs() < 1e-12 {
return gaussian_pdf(x, mu, sigma);
}
let beta = 2.0 / lambda;
let z = ((x - mu) / sigma).abs();
let norm = beta / (2.0 * sigma * gamma(1.0 / beta));
norm * (-z.powf(beta)).exp()
}
#[inline]
pub fn ecld_log_pdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
ecld_pdf(x, mu, sigma, lambda).ln()
}
#[inline]
pub fn ecld_cdf(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
if (lambda - 1.0).abs() < 1e-12 {
return Normal::new(mu, sigma).map(|n| n.cdf(x)).unwrap_or(0.5);
}
let beta = 2.0 / lambda;
let u = (x - mu) / sigma;
let p = gamma_lr(1.0 / beta, u.abs().powf(beta));
if u >= 0.0 {
0.5 + 0.5 * p
} else {
0.5 - 0.5 * p
}
}
#[inline]
pub fn ecld_variance(sigma: f64, lambda: f64) -> f64 {
if (lambda - 1.0).abs() < 1e-12 {
return sigma * sigma;
}
let beta = 2.0 / lambda;
let s2 = sigma * sigma;
s2 * gamma(3.0 / beta) / gamma(1.0 / beta)
}
#[inline]
fn gaussian_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
let variance = sigma * sigma;
let denom = (2.0 * std::f64::consts::PI * variance).sqrt();
let exponent = -((x - mu).powi(2)) / (2.0 * variance);
exponent.exp() / denom
}
pub fn work_to_natural(mu: f64, log_sigma: f64, log_lambda: f64) -> (f64, f64, f64) {
(mu, log_sigma.exp(), log_lambda.exp())
}
pub fn natural_to_work(mu: f64, sigma: f64, lambda: f64) -> (f64, f64, f64) {
(mu, sigma.ln(), lambda.ln())
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn lambda_one_matches_gaussian() {
let x = 0.005;
let mu = 0.008;
let sigma = 0.018;
let g = gaussian_pdf(x, mu, sigma);
let l = ecld_pdf(x, mu, sigma, 1.0);
assert_relative_eq!(l, g, epsilon = 1e-15);
}
#[test]
fn leptokurtic_lambda_increases_center_density() {
let x = 0.008; let mu = 0.008;
let sigma = 0.018;
let g = ecld_pdf(x, mu, sigma, 1.0);
let l = ecld_pdf(x, mu, sigma, 1.3);
assert!(l > g);
}
}