rv/dist/
cauchy.rs

1//! Cauchy distribution over x in (-∞, ∞)
2#[cfg(feature = "serde1")]
3use serde::{Deserialize, Serialize};
4
5use crate::consts::LN_PI;
6use crate::impl_display;
7use crate::traits::{
8    Cdf, ContinuousDistr, Entropy, HasDensity, InverseCdf, Median, Mode,
9    Parameterized, Sampleable, Scalable, Shiftable, Support,
10};
11use rand::Rng;
12use rand_distr::Cauchy as RCauchy;
13use std::f64::consts::{FRAC_1_PI, PI};
14use std::fmt;
15
16/// [Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution)
17/// over x in (-∞, ∞).
18///
19/// # Example
20/// ```
21/// use rv::prelude::*;
22///
23/// let cauchy = Cauchy::new(1.2, 3.4).expect("Invalid params");
24/// let ln_fx = cauchy.ln_pdf(&0.2_f64); // -2.4514716152673368
25///
26/// assert!((ln_fx + 2.4514716152673368).abs() < 1E-12);
27/// ```
28#[derive(Debug, Clone, PartialEq)]
29#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
30#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
31pub struct Cauchy {
32    /// location, x<sub>0</sub>, in (-∞, ∞)
33    loc: f64,
34    /// scale, γ, in (0, ∞)
35    scale: f64,
36}
37pub struct CauchyParameters {
38    pub loc: f64,
39    pub scale: f64,
40}
41
42impl Parameterized for Cauchy {
43    type Parameters = CauchyParameters;
44
45    fn emit_params(&self) -> Self::Parameters {
46        Self::Parameters {
47            loc: self.loc(),
48            scale: self.scale(),
49        }
50    }
51
52    fn from_params(params: Self::Parameters) -> Self {
53        Self::new_unchecked(params.loc, params.scale)
54    }
55}
56
57#[derive(Debug, Clone, PartialEq)]
58#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
59#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
60pub enum CauchyError {
61    /// Location parameter is infinite or NaN
62    LocNotFinite { loc: f64 },
63    /// Scale parameter is less than or equal to zero
64    ScaleTooLow { scale: f64 },
65    /// Scale parameter is infinite or Na]
66    ScaleNotFinite { scale: f64 },
67}
68
69impl Cauchy {
70    /// Creates a new Cauchy distribution
71    ///
72    /// # Arguments
73    /// - loc: location, x<sub>0</sub>, in (-∞, ∞)
74    /// - scale: scale, γ, in (0, ∞)
75    pub fn new(loc: f64, scale: f64) -> Result<Self, CauchyError> {
76        if !loc.is_finite() {
77            Err(CauchyError::LocNotFinite { loc })
78        } else if scale <= 0.0 {
79            Err(CauchyError::ScaleTooLow { scale })
80        } else if !scale.is_finite() {
81            Err(CauchyError::ScaleNotFinite { scale })
82        } else {
83            Ok(Cauchy { loc, scale })
84        }
85    }
86
87    /// Create a new Cauchy without checking whether the parameters are valid.
88    #[inline]
89    #[must_use]
90    pub fn new_unchecked(loc: f64, scale: f64) -> Self {
91        Cauchy { loc, scale }
92    }
93
94    /// Get the location parameter
95    ///
96    /// # Example
97    ///
98    /// ```rust
99    /// # use rv::dist::Cauchy;
100    /// let c = Cauchy::new(0.1, 1.0).unwrap();
101    /// assert_eq!(c.loc(), 0.1);
102    /// ```
103    #[inline]
104    #[must_use]
105    pub fn loc(&self) -> f64 {
106        self.loc
107    }
108
109    /// Set the location parameter
110    ///
111    /// # Example
112    ///
113    /// ```rust
114    /// # use rv::dist::Cauchy;
115    /// let mut c = Cauchy::new(0.1, 1.0).unwrap();
116    /// assert_eq!(c.loc(), 0.1);
117    ///
118    /// c.set_loc(2.0).unwrap();
119    /// assert_eq!(c.loc(), 2.0);
120    /// ```
121    /// Will error for invalid parameters
122    ///
123    /// ```rust
124    /// # use rv::dist::Cauchy;
125    /// # let mut c = Cauchy::new(0.1, 1.0).unwrap();
126    /// assert!(c.set_loc(2.0).is_ok());
127    /// assert!(c.set_loc(f64::INFINITY).is_err());
128    /// assert!(c.set_loc(f64::NEG_INFINITY).is_err());
129    /// assert!(c.set_loc(f64::NAN).is_err());
130    /// ```
131    #[inline]
132    pub fn set_loc(&mut self, loc: f64) -> Result<(), CauchyError> {
133        if loc.is_finite() {
134            self.set_loc_unchecked(loc);
135            Ok(())
136        } else {
137            Err(CauchyError::LocNotFinite { loc })
138        }
139    }
140
141    /// Set the location parameter without input validation
142    #[inline]
143    pub fn set_loc_unchecked(&mut self, loc: f64) {
144        self.loc = loc;
145    }
146
147    /// Get the scale parameter
148    ///
149    /// # Example
150    ///
151    /// ```rust
152    /// # use rv::dist::Cauchy;
153    /// let c = Cauchy::new(0.1, 1.0).unwrap();
154    /// assert_eq!(c.scale(), 1.0);
155    /// ```
156    #[inline]
157    #[must_use]
158    pub fn scale(&self) -> f64 {
159        self.scale
160    }
161
162    /// Set the scale parameter
163    ///
164    /// # Example
165    ///
166    /// ```rust
167    /// # use rv::dist::Cauchy;
168    /// let mut c = Cauchy::new(0.1, 1.0).unwrap();
169    /// assert_eq!(c.scale(), 1.0);
170    ///
171    /// c.set_scale(2.1).unwrap();
172    /// assert_eq!(c.scale(), 2.1);
173    /// ```
174    ///
175    /// Will error for invalid scale.
176    /// ```
177    /// # use rv::dist::Cauchy;
178    /// # let mut c = Cauchy::new(0.1, 1.0).unwrap();
179    /// assert!(c.set_scale(0.0).is_err());
180    /// assert!(c.set_scale(-1.0).is_err());
181    /// assert!(c.set_scale(f64::NAN).is_err());
182    /// assert!(c.set_scale(f64::INFINITY).is_err());
183    /// ```
184    #[inline]
185    pub fn set_scale(&mut self, scale: f64) -> Result<(), CauchyError> {
186        if !scale.is_finite() {
187            Err(CauchyError::ScaleNotFinite { scale })
188        } else if scale > 0.0 {
189            self.set_scale_unchecked(scale);
190            Ok(())
191        } else {
192            Err(CauchyError::ScaleTooLow { scale })
193        }
194    }
195
196    /// Set scale parameter without input validation
197    #[inline]
198    pub fn set_scale_unchecked(&mut self, scale: f64) {
199        self.scale = scale;
200    }
201}
202
203impl Default for Cauchy {
204    fn default() -> Self {
205        Cauchy::new_unchecked(0.0, 1.0)
206    }
207}
208
209impl From<&Cauchy> for String {
210    fn from(cauchy: &Cauchy) -> String {
211        format!("Cauchy(loc: {}, scale: {})", cauchy.loc, cauchy.scale)
212    }
213}
214
215impl_display!(Cauchy);
216use crate::misc::logaddexp;
217macro_rules! impl_traits {
218    ($kind:ty) => {
219        impl HasDensity<$kind> for Cauchy {
220            fn ln_f(&self, x: &$kind) -> f64 {
221                let ln_scale = self.scale.ln();
222                let term = 2.0_f64.mul_add(
223                    ((f64::from(*x) - self.loc).abs().ln() - ln_scale),
224                    ln_scale,
225                );
226                -logaddexp(ln_scale, term) - LN_PI
227            }
228        }
229
230        impl Sampleable<$kind> for Cauchy {
231            fn draw<R: Rng>(&self, rng: &mut R) -> $kind {
232                let cauchy = RCauchy::new(self.loc, self.scale).unwrap();
233                rng.sample(cauchy) as $kind
234            }
235
236            fn sample<R: Rng>(&self, n: usize, rng: &mut R) -> Vec<$kind> {
237                let cauchy = RCauchy::new(self.loc, self.scale).unwrap();
238                (0..n).map(|_| rng.sample(cauchy) as $kind).collect()
239            }
240        }
241
242        impl Support<$kind> for Cauchy {
243            fn supports(&self, x: &$kind) -> bool {
244                x.is_finite()
245            }
246        }
247
248        impl ContinuousDistr<$kind> for Cauchy {}
249
250        impl Cdf<$kind> for Cauchy {
251            fn cdf(&self, x: &$kind) -> f64 {
252                FRAC_1_PI.mul_add(
253                    ((f64::from(*x) - self.loc) / self.scale).atan(),
254                    0.5,
255                )
256            }
257        }
258
259        impl InverseCdf<$kind> for Cauchy {
260            fn invcdf(&self, p: f64) -> $kind {
261                self.scale.mul_add((PI * (p - 0.5)).tan(), self.loc) as $kind
262            }
263        }
264
265        impl Median<$kind> for Cauchy {
266            fn median(&self) -> Option<$kind> {
267                Some(self.loc as $kind)
268            }
269        }
270
271        impl Mode<$kind> for Cauchy {
272            fn mode(&self) -> Option<$kind> {
273                Some(self.loc as $kind)
274            }
275        }
276    };
277}
278
279impl Entropy for Cauchy {
280    fn entropy(&self) -> f64 {
281        (4.0 * PI * self.scale).ln()
282    }
283}
284
285impl_traits!(f64);
286impl_traits!(f32);
287
288impl std::error::Error for CauchyError {}
289
290#[cfg_attr(coverage_nightly, coverage(off))]
291impl fmt::Display for CauchyError {
292    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
293        match self {
294            Self::LocNotFinite { loc } => {
295                write!(f, "loc ({loc}) must be finite")
296            }
297            Self::ScaleTooLow { scale } => {
298                write!(f, "scale ({scale}) must be greater than zero")
299            }
300            Self::ScaleNotFinite { scale } => {
301                write!(f, "scale ({scale}) must be finite")
302            }
303        }
304    }
305}
306
307impl Shiftable for Cauchy {
308    type Output = Cauchy;
309    type Error = CauchyError;
310
311    fn shifted(self, shift: f64) -> Result<Self::Output, Self::Error>
312    where
313        Self: Sized,
314    {
315        Cauchy::new(self.loc() + shift, self.scale())
316    }
317
318    fn shifted_unchecked(self, shift: f64) -> Self::Output
319    where
320        Self: Sized,
321    {
322        Cauchy::new_unchecked(self.loc() + shift, self.scale())
323    }
324}
325impl Scalable for Cauchy {
326    type Output = Cauchy;
327    type Error = CauchyError;
328
329    fn scaled(self, scale: f64) -> Result<Self::Output, Self::Error>
330    where
331        Self: Sized,
332    {
333        Cauchy::new(self.loc() * scale, self.scale() * scale)
334    }
335
336    fn scaled_unchecked(self, scale: f64) -> Self::Output
337    where
338        Self: Sized,
339    {
340        Cauchy::new_unchecked(self.loc() * scale, self.scale() * scale)
341    }
342}
343
344#[cfg(test)]
345mod tests {
346    use super::*;
347    use crate::misc::ks_test;
348    use crate::test_basic_impls;
349    use std::f64;
350
351    const TOL: f64 = 1E-12;
352    const KS_PVAL: f64 = 0.2;
353    const N_TRIES: usize = 5;
354
355    test_basic_impls!(f64, Cauchy);
356
357    #[test]
358    fn ln_pdf_loc_zero() {
359        let c = Cauchy::new(0.0, 1.0).unwrap();
360        assert::close(c.ln_pdf(&0.2), -1.183_950_599_002_681_5, TOL);
361    }
362
363    #[test]
364    fn ln_pdf_loc_nonzero() {
365        let c = Cauchy::new(1.2, 3.4).unwrap();
366        assert::close(c.ln_pdf(&0.2), -2.451_471_615_267_337, TOL);
367    }
368
369    #[test]
370    fn cdf_at_loc() {
371        let c = Cauchy::new(1.2, 3.4).unwrap();
372        assert::close(c.cdf(&1.2), 0.5, TOL);
373    }
374
375    #[test]
376    fn cdf_off_loc() {
377        let c = Cauchy::new(1.2, 3.4).unwrap();
378        assert::close(c.cdf(&2.2), 0.591_053_001_855_748_8, TOL);
379        assert::close(c.cdf(&0.2), 1.0 - 0.591_053_001_855_748_8, TOL);
380    }
381
382    #[test]
383    fn inv_cdf_ident() {
384        let mut rng = rand::rng();
385        let c = Cauchy::default();
386        for _ in 0..100 {
387            let x: f64 = c.draw(&mut rng);
388            let p: f64 = c.cdf(&x);
389            let y: f64 = c.invcdf(p);
390            assert::close(x, y, 1E-8); // allow a little more error
391        }
392    }
393
394    #[test]
395    fn inv_cdf() {
396        let c = Cauchy::new(1.2, 3.4).unwrap();
397        let lower: f64 = c.invcdf(0.4);
398        let upper: f64 = c.invcdf(0.6);
399        assert::close(lower, 0.095_273_032_808_118_61, TOL);
400        assert::close(upper, 2.304_726_967_191_881_3, TOL);
401    }
402
403    #[test]
404    fn median_should_be_loc() {
405        let m: f64 = Cauchy::new(1.2, 3.4).unwrap().median().unwrap();
406        assert::close(m, 1.2, TOL);
407    }
408
409    #[test]
410    fn mode_should_be_loc() {
411        let m: f64 = Cauchy::new(-0.2, 3.4).unwrap().median().unwrap();
412        assert::close(m, -0.2, TOL);
413    }
414
415    #[test]
416    fn finite_numbers_should_be_in_support() {
417        let c = Cauchy::default();
418        assert!(c.supports(&0.0_f64));
419        assert!(c.supports(&f64::MIN_POSITIVE));
420        assert!(c.supports(&f64::MAX));
421        assert!(c.supports(&f64::MIN));
422    }
423
424    #[test]
425    fn non_finite_numbers_should_not_be_in_support() {
426        let c = Cauchy::default();
427        assert!(!c.supports(&f64::INFINITY));
428        assert!(!c.supports(&f64::NEG_INFINITY));
429        assert!(!c.supports(&f64::NAN));
430    }
431
432    #[test]
433    fn entropy() {
434        let c = Cauchy::new(1.2, 3.4).unwrap();
435        assert::close(c.entropy(), 3.754_799_678_591_406_4, TOL);
436    }
437
438    #[test]
439    fn loc_does_not_affect_entropy() {
440        let c1 = Cauchy::new(1.2, 3.4).unwrap();
441        let c2 = Cauchy::new(-99999.9, 3.4).unwrap();
442        assert::close(c1.entropy(), c2.entropy(), TOL);
443    }
444
445    #[test]
446    fn draw_test() {
447        let mut rng = rand::rng();
448        let c = Cauchy::new(1.2, 3.4).unwrap();
449        let cdf = |x: f64| c.cdf(&x);
450
451        // test is flaky, try a few times
452        let passes = (0..N_TRIES).fold(0, |acc, _| {
453            let xs: Vec<f64> = c.sample(1000, &mut rng);
454            let (_, p) = ks_test(&xs, cdf);
455            if p > KS_PVAL { acc + 1 } else { acc }
456        });
457
458        assert!(passes > 0);
459    }
460
461    use crate::test_shiftable_cdf;
462    use crate::test_shiftable_density;
463    use crate::test_shiftable_entropy;
464    use crate::test_shiftable_invcdf;
465    use crate::test_shiftable_method;
466
467    test_shiftable_method!(Cauchy::new(2.0, 4.0).unwrap(), median);
468    test_shiftable_density!(Cauchy::new(2.0, 4.0).unwrap());
469    test_shiftable_entropy!(Cauchy::new(2.0, 4.0).unwrap());
470    test_shiftable_cdf!(Cauchy::new(2.0, 4.0).unwrap());
471    test_shiftable_invcdf!(Cauchy::new(2.0, 4.0).unwrap());
472
473    use crate::test_scalable_cdf;
474    use crate::test_scalable_density;
475    use crate::test_scalable_entropy;
476    use crate::test_scalable_invcdf;
477    use crate::test_scalable_method;
478
479    test_scalable_method!(Cauchy::new(2.0, 4.0).unwrap(), median);
480    test_scalable_density!(Cauchy::new(2.0, 4.0).unwrap());
481    test_scalable_entropy!(Cauchy::new(2.0, 4.0).unwrap());
482    test_scalable_cdf!(Cauchy::new(2.0, 4.0).unwrap());
483    test_scalable_invcdf!(Cauchy::new(2.0, 4.0).unwrap());
484
485    #[test]
486    fn emit_and_from_params_are_identity() {
487        let dist_a = Cauchy::new(2.0, 4.0).unwrap();
488        let dist_b = Cauchy::from_params(dist_a.emit_params());
489        assert_eq!(dist_a, dist_b);
490    }
491}
rv/dist/cauchy.rs

rv/dist/
cauchy.rs
