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

rv/dist/
cauchy.rs
