Skip to main content

stats_claw/distributions/symmetric/
cauchy.rs

1//! Cauchy (Lorentz) distribution numerics, for the [`CauchyDistribution`] parameter struct.
2//!
3//! Equivalent to `scipy.stats.cauchy(loc=location, scale=scale)`: a symmetric,
4//! heavy-tailed distribution with no finite moments. The CDF/quantile are the
5//! closed-form arctangent/tangent pair, and sampling is the inverse-CDF transform
6//! of a single uniform draw.
7//!
8//! # Examples
9//!
10//! ```
11//! use stats_claw::distributions::Pdf;
12//! use stats_claw::distributions::CauchyDistribution;
13//! use std::f64::consts::PI;
14//!
15//! let d = CauchyDistribution { location: 0.0, scale: 1.0, ..Default::default() };
16//! // The standard Cauchy peaks at 1/π at the location.
17//! assert!((d.pdf(0.0) - 1.0 / PI).abs() < 1e-12, "peak was {}", d.pdf(0.0));
18//! ```
19
20use super::super::{Cdf, Moments, Pdf, Quantile, Sample};
21use crate::distributions::CauchyDistribution;
22use crate::rng::SplitMix64;
23use std::f64::consts::PI;
24
25impl Pdf for CauchyDistribution {
26    fn pdf(&self, x: f64) -> f64 {
27        let z = (x - self.location) / self.scale;
28        1.0 / (PI * self.scale * z.mul_add(z, 1.0))
29    }
30}
31
32impl Cdf for CauchyDistribution {
33    fn cdf(&self, x: f64) -> f64 {
34        let z = (x - self.location) / self.scale;
35        0.5 + z.atan() / PI
36    }
37}
38
39impl Quantile for CauchyDistribution {
40    fn quantile(&self, p: f64) -> f64 {
41        self.scale.mul_add((PI * (p - 0.5)).tan(), self.location)
42    }
43}
44
45impl Moments for CauchyDistribution {
46    /// Undefined — the Cauchy integral for the mean does not converge.
47    fn mean(&self) -> Option<f64> {
48        None
49    }
50    /// Undefined — the Cauchy second moment does not converge.
51    fn variance(&self) -> Option<f64> {
52        None
53    }
54}
55
56impl Sample for CauchyDistribution {
57    fn sample(&self, rng: &mut SplitMix64) -> f64 {
58        self.quantile(rng.next_f64())
59    }
60}
61
62#[cfg(test)]
63mod tests {
64    use super::*;
65
66    /// The standard Cauchy peaks at `1/π` at the location.
67    #[test]
68    fn peak_density_at_location() {
69        let d = CauchyDistribution {
70            location: 0.0,
71            scale: 1.0,
72            ..Default::default()
73        };
74        assert!(
75            (d.pdf(0.0) - 1.0 / PI).abs() < 1e-12,
76            "peak was {}",
77            d.pdf(0.0)
78        );
79    }
80
81    /// Both moments are undefined regardless of parameters.
82    #[test]
83    fn moments_are_none() {
84        let d = CauchyDistribution {
85            location: 3.0,
86            scale: 2.0,
87            ..Default::default()
88        };
89        assert!(d.mean().is_none() && d.variance().is_none());
90    }
91}