Skip to main content

stats_claw/distributions/symmetric/
cauchy.rs

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