1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::statistics::*;
3use std::f64;
4
5#[derive(Copy, Clone, PartialEq, Debug)]
19pub struct Cauchy {
20 location: f64,
21 scale: f64,
22}
23
24#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
26#[non_exhaustive]
27pub enum CauchyError {
28 LocationInvalid,
30
31 ScaleInvalid,
33}
34
35impl std::fmt::Display for CauchyError {
36 #[cfg_attr(coverage_nightly, coverage(off))]
37 fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
38 match self {
39 CauchyError::LocationInvalid => write!(f, "Location is NaN"),
40 CauchyError::ScaleInvalid => write!(f, "Scale is NaN, zero or less than zero"),
41 }
42 }
43}
44
45impl std::error::Error for CauchyError {}
46
47impl Cauchy {
48 pub fn new(location: f64, scale: f64) -> Result<Cauchy, CauchyError> {
67 if location.is_nan() {
68 return Err(CauchyError::LocationInvalid);
69 }
70
71 if scale.is_nan() || scale <= 0.0 {
72 return Err(CauchyError::ScaleInvalid);
73 }
74
75 Ok(Cauchy { location, scale })
76 }
77
78 pub fn location(&self) -> f64 {
89 self.location
90 }
91
92 pub fn scale(&self) -> f64 {
103 self.scale
104 }
105}
106
107impl std::fmt::Display for Cauchy {
108 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
109 write!(f, "Cauchy({}, {})", self.location, self.scale)
110 }
111}
112
113#[cfg(feature = "rand")]
114#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
115impl ::rand::distributions::Distribution<f64> for Cauchy {
116 fn sample<R: ::rand::Rng + ?Sized>(&self, r: &mut R) -> f64 {
117 self.location + self.scale * (f64::consts::PI * (r.gen::<f64>() - 0.5)).tan()
118 }
119}
120
121impl ContinuousCDF<f64, f64> for Cauchy {
122 fn cdf(&self, x: f64) -> f64 {
133 (1.0 / f64::consts::PI) * ((x - self.location) / self.scale).atan() + 0.5
134 }
135
136 fn sf(&self, x: f64) -> f64 {
149 (1.0 / f64::consts::PI) * ((self.location - x) / self.scale).atan() + 0.5
150 }
151
152 fn inverse_cdf(&self, x: f64) -> f64 {
163 if !(0.0..=1.0).contains(&x) {
164 panic!("x must be in [0, 1]");
165 } else {
166 self.location + self.scale * (f64::consts::PI * (x - 0.5)).tan()
167 }
168 }
169}
170
171impl Min<f64> for Cauchy {
172 fn min(&self) -> f64 {
181 f64::NEG_INFINITY
182 }
183}
184
185impl Max<f64> for Cauchy {
186 fn max(&self) -> f64 {
195 f64::INFINITY
196 }
197}
198
199impl Distribution<f64> for Cauchy {
200 fn entropy(&self) -> Option<f64> {
210 Some((4.0 * f64::consts::PI * self.scale).ln())
211 }
212}
213
214impl Median<f64> for Cauchy {
215 fn median(&self) -> f64 {
225 self.location
226 }
227}
228
229impl Mode<Option<f64>> for Cauchy {
230 fn mode(&self) -> Option<f64> {
240 Some(self.location)
241 }
242}
243
244impl Continuous<f64, f64> for Cauchy {
245 fn pdf(&self, x: f64) -> f64 {
256 1.0 / (f64::consts::PI
257 * self.scale
258 * (1.0 + ((x - self.location) / self.scale) * ((x - self.location) / self.scale)))
259 }
260
261 fn ln_pdf(&self, x: f64) -> f64 {
272 -(f64::consts::PI
273 * self.scale
274 * (1.0 + ((x - self.location) / self.scale) * ((x - self.location) / self.scale)))
275 .ln()
276 }
277}
278
279#[rustfmt::skip]
280#[cfg(test)]
281mod tests {
282 use super::*;
283 use crate::distribution::internal::*;
284 use crate::testing_boiler;
285
286 testing_boiler!(location: f64, scale: f64; Cauchy; CauchyError);
287
288 #[test]
289 fn test_create() {
290 create_ok(0.0, 0.1);
291 create_ok(0.0, 1.0);
292 create_ok(0.0, 10.0);
293 create_ok(10.0, 11.0);
294 create_ok(-5.0, 100.0);
295 create_ok(0.0, f64::INFINITY);
296 }
297
298 #[test]
299 fn test_bad_create() {
300 let invalid = [
301 (f64::NAN, 1.0, CauchyError::LocationInvalid),
302 (1.0, f64::NAN, CauchyError::ScaleInvalid),
303 (f64::NAN, f64::NAN, CauchyError::LocationInvalid),
304 (1.0, 0.0, CauchyError::ScaleInvalid),
305 ];
306
307 for (location, scale, err) in invalid {
308 test_create_err(location, scale, err);
309 }
310 }
311
312 #[test]
313 fn test_entropy() {
314 let entropy = |x: Cauchy| x.entropy().unwrap();
315 test_exact(0.0, 2.0, 3.224171427529236102395, entropy);
316 test_exact(0.1, 4.0, 3.917318608089181411812, entropy);
317 test_exact(1.0, 10.0, 4.833609339963336476996, entropy);
318 test_exact(10.0, 11.0, 4.92891951976766133704, entropy);
319 }
320
321 #[test]
322 fn test_mode() {
323 let mode = |x: Cauchy| x.mode().unwrap();
324 test_exact(0.0, 2.0, 0.0, mode);
325 test_exact(0.1, 4.0, 0.1, mode);
326 test_exact(1.0, 10.0, 1.0, mode);
327 test_exact(10.0, 11.0, 10.0, mode);
328 test_exact(0.0, f64::INFINITY, 0.0, mode);
329 }
330
331 #[test]
332 fn test_median() {
333 let median = |x: Cauchy| x.median();
334 test_exact(0.0, 2.0, 0.0, median);
335 test_exact(0.1, 4.0, 0.1, median);
336 test_exact(1.0, 10.0, 1.0, median);
337 test_exact(10.0, 11.0, 10.0, median);
338 test_exact(0.0, f64::INFINITY, 0.0, median);
339 }
340
341 #[test]
342 fn test_min_max() {
343 let min = |x: Cauchy| x.min();
344 let max = |x: Cauchy| x.max();
345 test_exact(0.0, 1.0, f64::NEG_INFINITY, min);
346 test_exact(0.0, 1.0, f64::INFINITY, max);
347 }
348
349 #[test]
350 fn test_pdf() {
351 let pdf = |arg: f64| move |x: Cauchy| x.pdf(arg);
352 test_exact(0.0, 0.1, 0.001272730452554141029739, pdf(-5.0));
353 test_exact(0.0, 0.1, 0.03151583031522679916216, pdf(-1.0));
354 test_absolute(0.0, 0.1, 3.183098861837906715378, 1e-14, pdf(0.0));
355 test_exact(0.0, 0.1, 0.03151583031522679916216, pdf(1.0));
356 test_exact(0.0, 0.1, 0.001272730452554141029739, pdf(5.0));
357 test_absolute(0.0, 1.0, 0.01224268793014579505914, 1e-17, pdf(-5.0));
358 test_exact(0.0, 1.0, 0.1591549430918953357689, pdf(-1.0));
359 test_exact(0.0, 1.0, 0.3183098861837906715378, pdf(0.0));
360 test_exact(0.0, 1.0, 0.1591549430918953357689, pdf(1.0));
361 test_absolute(0.0, 1.0, 0.01224268793014579505914, 1e-17, pdf(5.0));
362 test_exact(0.0, 10.0, 0.02546479089470325372302, pdf(-5.0));
363 test_exact(0.0, 10.0, 0.03151583031522679916216, pdf(-1.0));
364 test_exact(0.0, 10.0, 0.03183098861837906715378, pdf(0.0));
365 test_exact(0.0, 10.0, 0.03151583031522679916216, pdf(1.0));
366 test_exact(0.0, 10.0, 0.02546479089470325372302, pdf(5.0));
367 test_exact(-5.0, 100.0, 0.003183098861837906715378, pdf(-5.0));
368 test_absolute(-5.0, 100.0, 0.003178014039374906864395, 1e-17, pdf(-1.0));
369 test_exact(-5.0, 100.0, 0.003175160959439308444267, pdf(0.0));
370 test_exact(-5.0, 100.0, 0.003171680810918599756255, pdf(1.0));
371 test_absolute(-5.0, 100.0, 0.003151583031522679916216, 1e-17, pdf(5.0));
372 test_exact(0.0, f64::INFINITY, 0.0, pdf(-5.0));
373 test_exact(0.0, f64::INFINITY, 0.0, pdf(-1.0));
374 test_exact(0.0, f64::INFINITY, 0.0, pdf(0.0));
375 test_exact(0.0, f64::INFINITY, 0.0, pdf(1.0));
376 test_exact(0.0, f64::INFINITY, 0.0, pdf(5.0));
377 test_exact(f64::INFINITY, 1.0, 0.0, pdf(-5.0));
378 test_exact(f64::INFINITY, 1.0, 0.0, pdf(-1.0));
379 test_exact(f64::INFINITY, 1.0, 0.0, pdf(0.0));
380 test_exact(f64::INFINITY, 1.0, 0.0, pdf(1.0));
381 test_exact(f64::INFINITY, 1.0, 0.0, pdf(5.0));
382 }
383
384 #[test]
385 fn test_ln_pdf() {
386 let ln_pdf = |arg: f64| move |x: Cauchy| x.ln_pdf(arg);
387 test_exact(0.0, 0.1, -6.666590723732973542744, ln_pdf(-5.0));
388 test_absolute(0.0, 0.1, -3.457265309696613941009, 1e-14, ln_pdf(-1.0));
389 test_exact(0.0, 0.1, 1.157855207144645509875, ln_pdf(0.0));
390 test_absolute(0.0, 0.1, -3.457265309696613941009, 1e-14, ln_pdf(1.0));
391 test_exact(0.0, 0.1, -6.666590723732973542744, ln_pdf(5.0));
392 test_exact(0.0, 1.0, -4.402826423870882219615, ln_pdf(-5.0));
393 test_absolute(0.0, 1.0, -1.837877066409345483561, 1e-15, ln_pdf(-1.0));
394 test_exact(0.0, 1.0, -1.144729885849400174143, ln_pdf(0.0));
395 test_absolute(0.0, 1.0, -1.837877066409345483561, 1e-15, ln_pdf(1.0));
396 test_exact(0.0, 1.0, -4.402826423870882219615, ln_pdf(5.0));
397 test_exact(0.0, 10.0, -3.670458530157655613928, ln_pdf(-5.0));
398 test_absolute(0.0, 10.0, -3.457265309696613941009, 1e-14, ln_pdf(-1.0));
399 test_exact(0.0, 10.0, -3.447314978843445858161, ln_pdf(0.0));
400 test_absolute(0.0, 10.0, -3.457265309696613941009, 1e-14, ln_pdf(1.0));
401 test_exact(0.0, 10.0, -3.670458530157655613928, ln_pdf(5.0));
402 test_exact(-5.0, 100.0, -5.749900071837491542179, ln_pdf(-5.0));
403 test_exact(-5.0, 100.0, -5.751498793201188569872, ln_pdf(-1.0));
404 test_exact(-5.0, 100.0, -5.75239695203607874116, ln_pdf(0.0));
405 test_exact(-5.0, 100.0, -5.75349360734762171285, ln_pdf(1.0));
406 test_exact(-5.0, 100.0, -5.759850402690659625027, ln_pdf(5.0));
407 test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-5.0));
408 test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-1.0));
409 test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.0));
410 test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(1.0));
411 test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(5.0));
412 test_exact(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(-5.0));
413 test_exact(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(-1.0));
414 test_exact(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
415 test_exact(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(1.0));
416 test_exact(f64::INFINITY, 1.0, f64::NEG_INFINITY, ln_pdf(5.0));
417 }
418
419 #[test]
420 fn test_cdf() {
421 let cdf = |arg: f64| move |x: Cauchy| x.cdf(arg);
422 test_absolute(0.0, 0.1, 0.006365349100972796679298, 1e-16, cdf(-5.0));
423 test_absolute(0.0, 0.1, 0.03172551743055356951498, 1e-16, cdf(-1.0));
424 test_exact(0.0, 0.1, 0.5, cdf(0.0));
425 test_exact(0.0, 0.1, 0.968274482569446430485, cdf(1.0));
426 test_exact(0.0, 0.1, 0.9936346508990272033207, cdf(5.0));
427 test_absolute(0.0, 1.0, 0.06283295818900118381375, 1e-16, cdf(-5.0));
428 test_exact(0.0, 1.0, 0.25, cdf(-1.0));
429 test_exact(0.0, 1.0, 0.5, cdf(0.0));
430 test_exact(0.0, 1.0, 0.75, cdf(1.0));
431 test_exact(0.0, 1.0, 0.9371670418109988161863, cdf(5.0));
432 test_exact(0.0, 10.0, 0.3524163823495667258246, cdf(-5.0));
433 test_exact(0.0, 10.0, 0.468274482569446430485, cdf(-1.0));
434 test_exact(0.0, 10.0, 0.5, cdf(0.0));
435 test_exact(0.0, 10.0, 0.531725517430553569515, cdf(1.0));
436 test_exact(0.0, 10.0, 0.6475836176504332741754, cdf(5.0));
437 test_exact(-5.0, 100.0, 0.5, cdf(-5.0));
438 test_exact(-5.0, 100.0, 0.5127256113479918307809, cdf(-1.0));
439 test_exact(-5.0, 100.0, 0.5159022512561763751816, cdf(0.0));
440 test_exact(-5.0, 100.0, 0.5190757242358362337495, cdf(1.0));
441 test_exact(-5.0, 100.0, 0.531725517430553569515, cdf(5.0));
442 test_exact(0.0, f64::INFINITY, 0.5, cdf(-5.0));
443 test_exact(0.0, f64::INFINITY, 0.5, cdf(-1.0));
444 test_exact(0.0, f64::INFINITY, 0.5, cdf(0.0));
445 test_exact(0.0, f64::INFINITY, 0.5, cdf(1.0));
446 test_exact(0.0, f64::INFINITY, 0.5, cdf(5.0));
447 test_exact(f64::INFINITY, 1.0, 0.0, cdf(-5.0));
448 test_exact(f64::INFINITY, 1.0, 0.0, cdf(-1.0));
449 test_exact(f64::INFINITY, 1.0, 0.0, cdf(0.0));
450 test_exact(f64::INFINITY, 1.0, 0.0, cdf(1.0));
451 test_exact(f64::INFINITY, 1.0, 0.0, cdf(5.0));
452 }
453
454 #[test]
455 fn test_sf() {
456 let sf = |arg: f64| move |x: Cauchy| x.sf(arg);
457 test_absolute(0.0, 0.1, 0.9936346508990272, 1e-16, sf(-5.0));
458 test_absolute(0.0, 0.1, 0.9682744825694465, 1e-16, sf(-1.0));
459 test_exact(0.0, 0.1, 0.5, sf(0.0));
460 test_absolute(0.0, 0.1, 0.03172551743055352, 1e-16, sf(1.0));
461 test_exact(0.0, 0.1, 0.006365349100972806, sf(5.0));
462 test_absolute(0.0, 1.0, 0.9371670418109989, 1e-16, sf(-5.0));
463 test_exact(0.0, 1.0, 0.75, sf(-1.0));
464 test_exact(0.0, 1.0, 0.5, sf(0.0));
465 test_exact(0.0, 1.0, 0.25, sf(1.0));
466 test_exact(0.0, 1.0, 0.06283295818900114, sf(5.0));
467 test_exact(0.0, 10.0, 0.6475836176504333, sf(-5.0));
468 test_exact(0.0, 10.0, 0.5317255174305535, sf(-1.0));
469 test_exact(0.0, 10.0, 0.5, sf(0.0));
470 test_exact(0.0, 10.0, 0.4682744825694464, sf(1.0));
471 test_exact(0.0, 10.0, 0.35241638234956674, sf(5.0));
472 test_exact(-5.0, 100.0, 0.5, sf(-5.0));
473 test_exact(-5.0, 100.0, 0.4872743886520082, sf(-1.0));
474 test_exact(-5.0, 100.0, 0.4840977487438236, sf(0.0));
475 test_exact(-5.0, 100.0, 0.48092427576416374, sf(1.0));
476 test_exact(-5.0, 100.0, 0.4682744825694464, sf(5.0));
477 test_exact(0.0, f64::INFINITY, 0.5, sf(-5.0));
478 test_exact(0.0, f64::INFINITY, 0.5, sf(-1.0));
479 test_exact(0.0, f64::INFINITY, 0.5, sf(0.0));
480 test_exact(0.0, f64::INFINITY, 0.5, sf(1.0));
481 test_exact(0.0, f64::INFINITY, 0.5, sf(5.0));
482 test_exact(f64::INFINITY, 1.0, 1.0, sf(-5.0));
483 test_exact(f64::INFINITY, 1.0, 1.0, sf(-1.0));
484 test_exact(f64::INFINITY, 1.0, 1.0, sf(0.0));
485 test_exact(f64::INFINITY, 1.0, 1.0, sf(1.0));
486 test_exact(f64::INFINITY, 1.0, 1.0, sf(5.0));
487 }
488
489 #[test]
490 fn test_inverse_cdf() {
491 let func = |arg: f64| move |x: Cauchy| x.inverse_cdf(x.cdf(arg));
492 test_absolute(0.0, 0.1, -5.0, 1e-10, func(-5.0));
493 test_absolute(0.0, 0.1, -1.0, 1e-14, func(-1.0));
494 test_exact(0.0, 0.1, 0.0, func(0.0));
495 test_absolute(0.0, 0.1, 1.0, 1e-14, func(1.0));
496 test_absolute(0.0, 0.1, 5.0, 1e-10, func(5.0));
497 test_absolute(0.0, 1.0, -5.0, 1e-14, func(-5.0));
498 test_absolute(0.0, 1.0, -1.0, 1e-15, func(-1.0));
499 test_exact(0.0, 1.0, 0.0, func(0.0));
500 test_absolute(0.0, 1.0, 1.0, 1e-15, func(1.0));
501 test_absolute(0.0, 1.0, 5.0, 1e-14, func(5.0));
502 test_absolute(0.0, 10.0, -5.0, 1e-14, func(-5.0));
503 test_absolute(0.0, 10.0, -1.0, 1e-14, func(-1.0));
504 test_exact(0.0, 10.0, 0.0, func(0.0));
505 test_absolute(0.0, 10.0, 1.0, 1e-14, func(1.0));
506 test_absolute(0.0, 10.0, 5.0, 1e-14, func(5.0));
507 test_exact(-5.0, 100.0, -5.0, func(-5.0));
508 test_absolute(-5.0, 100.0, -1.0, 1e-10, func(-1.0));
509 test_absolute(-5.0, 100.0, 0.0, 1e-14, func(0.0));
510 test_absolute(-5.0, 100.0, 1.0, 1e-14, func(1.0));
511 test_absolute(-5.0, 100.0, 5.0, 1e-10, func(5.0));
512 }
513
514 #[test]
515 fn test_continuous() {
516 test::check_continuous_distribution(&create_ok(-1.2, 3.4), -1500.0, 1500.0);
517 test::check_continuous_distribution(&create_ok(-4.5, 6.7), -5000.0, 5000.0);
518 }
519}