statrs/distribution/
cauchy.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::statistics::*;
3use crate::{Result, StatsError};
4use rand::Rng;
5use std::f64;
6
7/// Implements the [Cauchy](https://en.wikipedia.org/wiki/Cauchy_distribution)
8/// distribution, also known as the Lorentz distribution.
9///
10/// # Examples
11///
12/// ```
13/// use statrs::distribution::{Cauchy, Continuous};
14/// use statrs::statistics::Mode;
15///
16/// let n = Cauchy::new(0.0, 1.0).unwrap();
17/// assert_eq!(n.mode().unwrap(), 0.0);
18/// assert_eq!(n.pdf(1.0), 0.1591549430918953357689);
19/// ```
20#[derive(Debug, Copy, Clone, PartialEq)]
21pub struct Cauchy {
22    location: f64,
23    scale: f64,
24}
25
26impl Cauchy {
27    /// Constructs a new cauchy distribution with the given
28    /// location and scale.
29    ///
30    /// # Errors
31    ///
32    /// Returns an error if location or scale are `NaN` or `scale <= 0.0`
33    ///
34    /// # Examples
35    ///
36    /// ```
37    /// use statrs::distribution::Cauchy;
38    ///
39    /// let mut result = Cauchy::new(0.0, 1.0);
40    /// assert!(result.is_ok());
41    ///
42    /// result = Cauchy::new(0.0, -1.0);
43    /// assert!(result.is_err());
44    /// ```
45    pub fn new(location: f64, scale: f64) -> Result<Cauchy> {
46        if location.is_nan() || scale.is_nan() || scale <= 0.0 {
47            Err(StatsError::BadParams)
48        } else {
49            Ok(Cauchy { location, scale })
50        }
51    }
52
53    /// Returns the location of the cauchy distribution
54    ///
55    /// # Examples
56    ///
57    /// ```
58    /// use statrs::distribution::Cauchy;
59    ///
60    /// let n = Cauchy::new(0.0, 1.0).unwrap();
61    /// assert_eq!(n.location(), 0.0);
62    /// ```
63    pub fn location(&self) -> f64 {
64        self.location
65    }
66
67    /// Returns the scale of the cauchy distribution
68    ///
69    /// # Examples
70    ///
71    /// ```
72    /// use statrs::distribution::Cauchy;
73    ///
74    /// let n = Cauchy::new(0.0, 1.0).unwrap();
75    /// assert_eq!(n.scale(), 1.0);
76    /// ```
77    pub fn scale(&self) -> f64 {
78        self.scale
79    }
80}
81
82impl ::rand::distributions::Distribution<f64> for Cauchy {
83    fn sample<R: Rng + ?Sized>(&self, r: &mut R) -> f64 {
84        self.location + self.scale * (f64::consts::PI * (r.gen::<f64>() - 0.5)).tan()
85    }
86}
87
88impl ContinuousCDF<f64, f64> for Cauchy {
89    /// Calculates the cumulative distribution function for the
90    /// cauchy distribution at `x`
91    ///
92    /// # Formula
93    ///
94    /// ```ignore
95    /// (1 / π) * arctan((x - x_0) / γ) + 0.5
96    /// ```
97    ///
98    /// where `x_0` is the location and `γ` is the scale
99    fn cdf(&self, x: f64) -> f64 {
100        (1.0 / f64::consts::PI) * ((x - self.location) / self.scale).atan() + 0.5
101    }
102
103    /// Calculates the survival function for the
104    /// cauchy distribution at `x`
105    ///
106    /// # Formula
107    ///
108    /// ```ignore
109    /// (1 / π) * arctan(-(x - x_0) / γ) + 0.5
110    /// ```
111    ///
112    /// where `x_0` is the location and `γ` is the scale.
113    /// note that this is identical to the cdf except for
114    /// the negative argument to the arctan function
115    fn sf(&self, x: f64) -> f64 {
116        (1.0 / f64::consts::PI) * ((self.location - x) / self.scale).atan() + 0.5
117    }
118}
119
120impl Min<f64> for Cauchy {
121    /// Returns the minimum value in the domain of the cauchy
122    /// distribution representable by a double precision float
123    ///
124    /// # Formula
125    ///
126    /// ```ignore
127    /// NEG_INF
128    /// ```
129    fn min(&self) -> f64 {
130        f64::NEG_INFINITY
131    }
132}
133
134impl Max<f64> for Cauchy {
135    /// Returns the maximum value in the domain of the cauchy
136    /// distribution representable by a double precision float
137    ///
138    /// # Formula
139    ///
140    /// ```ignore
141    /// INF
142    /// ```
143    fn max(&self) -> f64 {
144        f64::INFINITY
145    }
146}
147
148impl Distribution<f64> for Cauchy {
149    /// Returns the entropy of the cauchy distribution
150    ///
151    /// # Formula
152    ///
153    /// ```ignore
154    /// ln(γ) + ln(4π)
155    /// ```
156    ///
157    /// where `γ` is the scale
158    fn entropy(&self) -> Option<f64> {
159        Some((4.0 * f64::consts::PI * self.scale).ln())
160    }
161}
162
163impl Median<f64> for Cauchy {
164    /// Returns the median of the cauchy distribution
165    ///
166    /// # Formula
167    ///
168    /// ```ignore
169    /// x_0
170    /// ```
171    ///
172    /// where `x_0` is the location
173    fn median(&self) -> f64 {
174        self.location
175    }
176}
177
178impl Mode<Option<f64>> for Cauchy {
179    /// Returns the mode of the cauchy distribution
180    ///
181    /// # Formula
182    ///
183    /// ```ignore
184    /// x_0
185    /// ```
186    ///
187    /// where `x_0` is the location
188    fn mode(&self) -> Option<f64> {
189        Some(self.location)
190    }
191}
192
193impl Continuous<f64, f64> for Cauchy {
194    /// Calculates the probability density function for the cauchy
195    /// distribution at `x`
196    ///
197    /// # Formula
198    ///
199    /// ```ignore
200    /// 1 / (πγ * (1 + ((x - x_0) / γ)^2))
201    /// ```
202    ///
203    /// where `x_0` is the location and `γ` is the scale
204    fn pdf(&self, x: f64) -> f64 {
205        1.0 / (f64::consts::PI
206            * self.scale
207            * (1.0 + ((x - self.location) / self.scale) * ((x - self.location) / self.scale)))
208    }
209
210    /// Calculates the log probability density function for the cauchy
211    /// distribution at `x`
212    ///
213    /// # Formula
214    ///
215    /// ```ignore
216    /// ln(1 / (πγ * (1 + ((x - x_0) / γ)^2)))
217    /// ```
218    ///
219    /// where `x_0` is the location and `γ` is the scale
220    fn ln_pdf(&self, x: f64) -> f64 {
221        -(f64::consts::PI
222            * self.scale
223            * (1.0 + ((x - self.location) / self.scale) * ((x - self.location) / self.scale)))
224            .ln()
225    }
226}
227
228#[rustfmt::skip]
229#[cfg(all(test, feature = "nightly"))]
230mod tests {
231    use crate::statistics::*;
232    use crate::distribution::{ContinuousCDF, Continuous, Cauchy};
233    use crate::distribution::internal::*;
234    use crate::consts::ACC;
235
236    fn try_create(location: f64, scale: f64) -> Cauchy {
237        let n = Cauchy::new(location, scale);
238        assert!(n.is_ok());
239        n.unwrap()
240    }
241
242    fn create_case(location: f64, scale: f64) {
243        let n = try_create(location, scale);
244        assert_eq!(location, n.location());
245        assert_eq!(scale, n.scale());
246    }
247
248    fn bad_create_case(location: f64, scale: f64) {
249        let n = Cauchy::new(location, scale);
250        assert!(n.is_err());
251    }
252
253    fn test_case<F>(location: f64, scale: f64, expected: f64, eval: F)
254        where F: Fn(Cauchy) -> f64
255    {
256        let n = try_create(location, scale);
257        let x = eval(n);
258        assert_eq!(expected, x);
259    }
260
261    fn test_almost<F>(location: f64, scale: f64, expected: f64, acc: f64, eval: F)
262        where F: Fn(Cauchy) -> f64
263    {
264        let n = try_create(location, scale);
265        let x = eval(n);
266        assert_almost_eq!(expected, x, acc);
267    }
268
269    #[test]
270    fn test_create() {
271        create_case(0.0, 0.1);
272        create_case(0.0, 1.0);
273        create_case(0.0, 10.0);
274        create_case(10.0, 11.0);
275        create_case(-5.0, 100.0);
276        create_case(0.0, f64::INFINITY);
277    }
278
279    #[test]
280    fn test_bad_create() {
281        bad_create_case(f64::NAN, 1.0);
282        bad_create_case(1.0, f64::NAN);
283        bad_create_case(f64::NAN, f64::NAN);
284        bad_create_case(1.0, 0.0);
285    }
286
287    #[test]
288    fn test_entropy() {
289        let entropy = |x: Cauchy| x.entropy().unwrap();
290        test_case(0.0, 2.0, 3.224171427529236102395, entropy);
291        test_case(0.1, 4.0, 3.917318608089181411812, entropy);
292        test_case(1.0, 10.0, 4.833609339963336476996, entropy);
293        test_case(10.0, 11.0, 4.92891951976766133704, entropy);
294    }
295
296    #[test]
297    fn test_mode() {
298        let mode = |x: Cauchy| x.mode().unwrap();
299        test_case(0.0, 2.0, 0.0, mode);
300        test_case(0.1, 4.0, 0.1, mode);
301        test_case(1.0, 10.0, 1.0, mode);
302        test_case(10.0, 11.0, 10.0, mode);
303        test_case(0.0, f64::INFINITY, 0.0, mode);
304    }
305
306    #[test]
307    fn test_median() {
308        let median = |x: Cauchy| x.median();
309        test_case(0.0, 2.0, 0.0, median);
310        test_case(0.1, 4.0, 0.1, median);
311        test_case(1.0, 10.0, 1.0, median);
312        test_case(10.0, 11.0, 10.0, median);
313        test_case(0.0, f64::INFINITY, 0.0, median);
314    }
315
316    #[test]
317    fn test_min_max() {
318        let min = |x: Cauchy| x.min();
319        let max = |x: Cauchy| x.max();
320        test_case(0.0, 1.0, f64::NEG_INFINITY, min);
321        test_case(0.0, 1.0, f64::INFINITY, max);
322    }
323
324    #[test]
325    fn test_pdf() {
326        let pdf = |arg: f64| move |x: Cauchy| x.pdf(arg);
327        test_case(0.0, 0.1, 0.001272730452554141029739, pdf(-5.0));
328        test_case(0.0, 0.1, 0.03151583031522679916216, pdf(-1.0));
329        test_almost(0.0, 0.1, 3.183098861837906715378, 1e-14, pdf(0.0));
330        test_case(0.0, 0.1, 0.03151583031522679916216, pdf(1.0));
331        test_case(0.0, 0.1, 0.001272730452554141029739, pdf(5.0));
332        test_almost(0.0, 1.0, 0.01224268793014579505914, 1e-17, pdf(-5.0));
333        test_case(0.0, 1.0, 0.1591549430918953357689, pdf(-1.0));
334        test_case(0.0, 1.0, 0.3183098861837906715378, pdf(0.0));
335        test_case(0.0, 1.0, 0.1591549430918953357689, pdf(1.0));
336        test_almost(0.0, 1.0, 0.01224268793014579505914, 1e-17, pdf(5.0));
337        test_case(0.0, 10.0, 0.02546479089470325372302, pdf(-5.0));
338        test_case(0.0, 10.0, 0.03151583031522679916216, pdf(-1.0));
339        test_case(0.0, 10.0, 0.03183098861837906715378, pdf(0.0));
340        test_case(0.0, 10.0, 0.03151583031522679916216, pdf(1.0));
341        test_case(0.0, 10.0, 0.02546479089470325372302, pdf(5.0));
342        test_case(-5.0, 100.0, 0.003183098861837906715378, pdf(-5.0));
343        test_almost(-5.0, 100.0, 0.003178014039374906864395, 1e-17, pdf(-1.0));
344        test_case(-5.0, 100.0, 0.003175160959439308444267, pdf(0.0));
345        test_case(-5.0, 100.0, 0.003171680810918599756255, pdf(1.0));
346        test_almost(-5.0, 100.0, 0.003151583031522679916216, 1e-17, pdf(5.0));
347        test_case(0.0, f64::INFINITY, 0.0, pdf(-5.0));
348        test_case(0.0, f64::INFINITY, 0.0, pdf(-1.0));
349        test_case(0.0, f64::INFINITY, 0.0, pdf(0.0));
350        test_case(0.0, f64::INFINITY, 0.0, pdf(1.0));
351        test_case(0.0, f64::INFINITY, 0.0, pdf(5.0));
352        test_case(f64::INFINITY, 1.0, 0.0, pdf(-5.0));
353        test_case(f64::INFINITY, 1.0, 0.0, pdf(-1.0));
354        test_case(f64::INFINITY, 1.0, 0.0, pdf(0.0));
355        test_case(f64::INFINITY, 1.0, 0.0, pdf(1.0));
356        test_case(f64::INFINITY, 1.0, 0.0, pdf(5.0));
357    }
358
359    #[test]
360    fn test_ln_pdf() {
361        let ln_pdf = |arg: f64| move |x: Cauchy| x.ln_pdf(arg);
362        test_case(0.0, 0.1, -6.666590723732973542744, ln_pdf(-5.0));
363        test_almost(0.0, 0.1, -3.457265309696613941009, 1e-14, ln_pdf(-1.0));
364        test_case(0.0, 0.1, 1.157855207144645509875, ln_pdf(0.0));
365        test_almost(0.0, 0.1, -3.457265309696613941009, 1e-14, ln_pdf(1.0));
366        test_case(0.0, 0.1, -6.666590723732973542744, ln_pdf(5.0));
367        test_case(0.0, 1.0, -4.402826423870882219615, ln_pdf(-5.0));
368        test_almost(0.0, 1.0, -1.837877066409345483561, 1e-15, ln_pdf(-1.0));
369        test_case(0.0, 1.0, -1.144729885849400174143, ln_pdf(0.0));
370        test_almost(0.0, 1.0, -1.837877066409345483561, 1e-15, ln_pdf(1.0));
371        test_case(0.0, 1.0, -4.402826423870882219615, ln_pdf(5.0));
372        test_case(0.0, 10.0, -3.670458530157655613928, ln_pdf(-5.0));
373        test_almost(0.0, 10.0, -3.457265309696613941009, 1e-14, ln_pdf(-1.0));
374        test_case(0.0, 10.0, -3.447314978843445858161, ln_pdf(0.0));
375        test_almost(0.0, 10.0, -3.457265309696613941009, 1e-14, ln_pdf(1.0));
376        test_case(0.0, 10.0, -3.670458530157655613928, ln_pdf(5.0));
377        test_case(-5.0, 100.0, -5.749900071837491542179, ln_pdf(-5.0));
378        test_case(-5.0, 100.0, -5.751498793201188569872, ln_pdf(-1.0));
379        test_case(-5.0, 100.0, -5.75239695203607874116, ln_pdf(0.0));
380        test_case(-5.0, 100.0, -5.75349360734762171285, ln_pdf(1.0));
381        test_case(-5.0, 100.0, -5.759850402690659625027, ln_pdf(5.0));
382        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-5.0));
383        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-1.0));
384        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.0));
385        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(1.0));
386        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(5.0));
387        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(-5.0));
388        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(-1.0));
389        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
390        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(1.0));
391        test_case(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(5.0));
392    }
393
394    #[test]
395    fn test_cdf() {
396        let cdf = |arg: f64| move |x: Cauchy| x.cdf(arg);
397        test_almost(0.0, 0.1, 0.006365349100972796679298, 1e-16, cdf(-5.0));
398        test_almost(0.0, 0.1, 0.03172551743055356951498, 1e-16, cdf(-1.0));
399        test_case(0.0, 0.1, 0.5, cdf(0.0));
400        test_case(0.0, 0.1, 0.968274482569446430485, cdf(1.0));
401        test_case(0.0, 0.1, 0.9936346508990272033207, cdf(5.0));
402        test_almost(0.0, 1.0, 0.06283295818900118381375, 1e-16, cdf(-5.0));
403        test_case(0.0, 1.0, 0.25, cdf(-1.0));
404        test_case(0.0, 1.0, 0.5, cdf(0.0));
405        test_case(0.0, 1.0, 0.75, cdf(1.0));
406        test_case(0.0, 1.0, 0.9371670418109988161863, cdf(5.0));
407        test_case(0.0, 10.0, 0.3524163823495667258246, cdf(-5.0));
408        test_case(0.0, 10.0, 0.468274482569446430485, cdf(-1.0));
409        test_case(0.0, 10.0, 0.5, cdf(0.0));
410        test_case(0.0, 10.0, 0.531725517430553569515, cdf(1.0));
411        test_case(0.0, 10.0, 0.6475836176504332741754, cdf(5.0));
412        test_case(-5.0, 100.0, 0.5, cdf(-5.0));
413        test_case(-5.0, 100.0, 0.5127256113479918307809, cdf(-1.0));
414        test_case(-5.0, 100.0, 0.5159022512561763751816, cdf(0.0));
415        test_case(-5.0, 100.0, 0.5190757242358362337495, cdf(1.0));
416        test_case(-5.0, 100.0, 0.531725517430553569515, cdf(5.0));
417        test_case(0.0, f64::INFINITY, 0.5, cdf(-5.0));
418        test_case(0.0, f64::INFINITY, 0.5, cdf(-1.0));
419        test_case(0.0, f64::INFINITY, 0.5, cdf(0.0));
420        test_case(0.0, f64::INFINITY, 0.5, cdf(1.0));
421        test_case(0.0, f64::INFINITY, 0.5, cdf(5.0));
422        test_case(f64::INFINITY, 1.0, 0.0, cdf(-5.0));
423        test_case(f64::INFINITY, 1.0, 0.0, cdf(-1.0));
424        test_case(f64::INFINITY, 1.0, 0.0, cdf(0.0));
425        test_case(f64::INFINITY, 1.0, 0.0, cdf(1.0));
426        test_case(f64::INFINITY, 1.0, 0.0, cdf(5.0));
427    }
428
429    #[test]
430    fn test_sf() {
431        let sf = |arg: f64| move |x: Cauchy| x.sf(arg);
432        test_almost(0.0, 0.1, 0.9936346508990272, 1e-16, sf(-5.0));
433        test_almost(0.0, 0.1, 0.9682744825694465, 1e-16, sf(-1.0));
434        test_case(0.0, 0.1, 0.5, sf(0.0));
435        test_case(0.0, 0.1, 0.03172551743055352, sf(1.0));
436        test_case(0.0, 0.1, 0.006365349100972806, sf(5.0));
437        test_almost(0.0, 1.0, 0.9371670418109989, 1e-16, sf(-5.0));
438        test_case(0.0, 1.0, 0.75, sf(-1.0));
439        test_case(0.0, 1.0, 0.5, sf(0.0));
440        test_case(0.0, 1.0, 0.25, sf(1.0));
441        test_case(0.0, 1.0, 0.06283295818900114, sf(5.0));
442        test_case(0.0, 10.0, 0.6475836176504333, sf(-5.0));
443        test_case(0.0, 10.0, 0.5317255174305535, sf(-1.0));
444        test_case(0.0, 10.0, 0.5, sf(0.0));
445        test_case(0.0, 10.0, 0.4682744825694464, sf(1.0));
446        test_case(0.0, 10.0, 0.35241638234956674, sf(5.0));
447        test_case(-5.0, 100.0, 0.5, sf(-5.0));
448        test_case(-5.0, 100.0, 0.4872743886520082, sf(-1.0));
449        test_case(-5.0, 100.0, 0.4840977487438236, sf(0.0));
450        test_case(-5.0, 100.0, 0.48092427576416374, sf(1.0));
451        test_case(-5.0, 100.0, 0.4682744825694464, sf(5.0));
452        test_case(0.0, f64::INFINITY, 0.5, sf(-5.0));
453        test_case(0.0, f64::INFINITY, 0.5, sf(-1.0));
454        test_case(0.0, f64::INFINITY, 0.5, sf(0.0));
455        test_case(0.0, f64::INFINITY, 0.5, sf(1.0));
456        test_case(0.0, f64::INFINITY, 0.5, sf(5.0));
457        test_case(f64::INFINITY, 1.0, 1.0, sf(-5.0));
458        test_case(f64::INFINITY, 1.0, 1.0, sf(-1.0));
459        test_case(f64::INFINITY, 1.0, 1.0, sf(0.0));
460        test_case(f64::INFINITY, 1.0, 1.0, sf(1.0));
461        test_case(f64::INFINITY, 1.0, 1.0, sf(5.0));
462    }
463
464    #[test]
465    fn test_continuous() {
466        test::check_continuous_distribution(&try_create(-1.2, 3.4), -1500.0, 1500.0);
467        test::check_continuous_distribution(&try_create(-4.5, 6.7), -5000.0, 5000.0);
468    }
469}
statrs/distribution/cauchy.rs

statrs/distribution/
cauchy.rs
