From 6f50ce3d9483ebb5910769870b3bddfeb5ab2c21 Mon Sep 17 00:00:00 2001 From: Tom Conway Date: Thu, 21 May 2026 11:26:25 +1000 Subject: [PATCH 1/3] Add the Cauchy distribution. --- src/distr/cauchy.rs | 215 ++++++++++++++++++++++++++++++++++++++++++++ src/distr/mod.rs | 2 + 2 files changed, 217 insertions(+) create mode 100644 src/distr/cauchy.rs diff --git a/src/distr/cauchy.rs b/src/distr/cauchy.rs new file mode 100644 index 0000000000..9fc974bbe5 --- /dev/null +++ b/src/distr/cauchy.rs @@ -0,0 +1,215 @@ +// Copyright 2018 Developers of the Rand project. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +//! The Cauchy distribution `Cauchy(x_0, γ)`. + +use crate::distr::Distribution; +use crate::{Rng, RngExt}; +use core::fmt; + +#[cfg(feature = "serde")] +use serde::{Deserialize, Serialize}; + +/// The [Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution) `Cauchy(x_0, γ)`. +/// +/// The Cauchy distribution is a continuous probability distribution with +/// location parameter `x_0` (median) and scale parameter `γ > 0` (half-width +/// at half-maximum). It has no finite mean or variance. +/// +/// # Plot +/// +/// The following plot shows the Cauchy distribution with `x_0 = 0` and +/// `γ = 1` (the standard Cauchy distribution). +/// +/// # Example +/// +/// ```rust +/// use rand::distr::{Cauchy, Distribution}; +/// +/// let d = Cauchy::new(0.0, 1.0).unwrap(); +/// let v = d.sample(&mut rand::rng()); +/// println!("{} is from a Cauchy distribution", v); +/// ``` +#[derive(Clone, Copy, Debug, PartialEq)] +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +pub struct Cauchy { + median: f64, + scale: f64, +} + +/// Error type returned from [`Cauchy::new`]. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum CauchyError { + /// `scale <= 0` or `scale` is NaN. + InvalidScale, + /// `median` is infinite or NaN. + InvalidMedian, +} + +impl fmt::Display for CauchyError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + f.write_str(match self { + CauchyError::InvalidScale => "scale is not positive and finite in Cauchy distribution", + CauchyError::InvalidMedian => "median is not finite in Cauchy distribution", + }) + } +} + +impl core::error::Error for CauchyError {} + +impl Cauchy { + /// Construct a new `Cauchy` with the given `median` (`x_0`) and `scale` (`γ`). + /// + /// # Errors + /// + /// Returns [`CauchyError::InvalidMedian`] if `median` is infinite or NaN, + /// or [`CauchyError::InvalidScale`] if `scale` is not strictly positive and finite. + #[inline] + pub fn new(median: f64, scale: f64) -> Result { + if !median.is_finite() { + return Err(CauchyError::InvalidMedian); + } + if !(scale > 0.0) || !scale.is_finite() { + return Err(CauchyError::InvalidScale); + } + Ok(Cauchy { median, scale }) + } + + /// Returns the median (`x_0`) of the distribution. + #[inline] + pub fn median(&self) -> f64 { + self.median + } + + /// Returns the scale (`γ`) of the distribution. + #[inline] + pub fn scale(&self) -> f64 { + self.scale + } +} + +impl Distribution for Cauchy { + fn sample(&self, rng: &mut R) -> f64 { + loop { + let u: f64 = 2.0 * rng.random::() - 1.0; + let v: f64 = 2.0 * rng.random::() - 1.0; + if u * u + v * v > 1.0 || v == 0.0 { + continue; + } + return self.median + self.scale * u / v; + } + } +} + +#[cfg(test)] +mod test { + use std::{println, vec::Vec}; + + use super::*; + use crate::RngExt; + + #[test] + fn test_new_valid() { + assert!(Cauchy::new(0.0, 1.0).is_ok()); + assert!(Cauchy::new(-5.0, 0.1).is_ok()); + } + + #[test] + fn test_new_invalid_scale() { + assert_eq!(Cauchy::new(0.0, 0.0), Err(CauchyError::InvalidScale)); + assert_eq!(Cauchy::new(0.0, -1.0), Err(CauchyError::InvalidScale)); + assert_eq!(Cauchy::new(0.0, f64::NAN), Err(CauchyError::InvalidScale)); + assert_eq!( + Cauchy::new(0.0, f64::INFINITY), + Err(CauchyError::InvalidScale) + ); + } + + #[test] + fn test_new_invalid_median() { + assert_eq!(Cauchy::new(f64::NAN, 1.0), Err(CauchyError::InvalidMedian)); + assert_eq!( + Cauchy::new(f64::INFINITY, 1.0), + Err(CauchyError::InvalidMedian) + ); + assert_eq!( + Cauchy::new(f64::NEG_INFINITY, 1.0), + Err(CauchyError::InvalidMedian) + ); + } + + #[test] + fn test_accessors() { + let d = Cauchy::new(2.0, 3.0).unwrap(); + assert_eq!(d.median(), 2.0); + assert_eq!(d.scale(), 3.0); + } + + #[test] + #[cfg(feature = "serde")] + fn test_serializing_deserializing_cauchy() { + let distr = Cauchy::new(0.0, 1.0).unwrap(); + let de: Cauchy = postcard::from_bytes(&postcard::to_allocvec(&distr).unwrap()).unwrap(); + assert_eq!(distr, de); + } + + #[test] + fn test_distr() { + let mut rng = crate::test::rng(1); + let distr = Cauchy::new(0.0, 1.0).unwrap(); + + const N: usize = 100_000; + let mut values: Vec = Vec::new(); + for _ in 0..N { + values.push(rng.sample(distr)); + } + values.sort_by(|a, b| a.partial_cmp(b).unwrap()); + + // Calculate the Kolmogorov-Smirnoff divergence between the + // empirical and expected distributions, and check it's small. + let mut d0 = 0.0; + for i in 0..N { + let p = (i as f64) / (N as f64); + let x = values[i]; + let q = x.atan() / std::f64::consts::PI + 0.5; + let d = (p - q).abs(); + if d > d0 { + d0 = d; + } + } + assert!(d0 < 0.005); + } + + #[test] + fn value_stability() { + let mut rng = crate::test::rng(1); + let distr = Cauchy::new(0.0, 1.0).unwrap(); + let mut buf = [0.0f64; 10]; + for x in &mut buf { + *x = rng.sample(distr); + } + // Fill in expected values once sample() is implemented. + let expected: [f64; 10] = [ + -0.18291029299250308, + -101.0791569450233, + 0.491799631323629, + 1.5187358196195335, + -0.4600057220965164, + 0.4052753943925359, + 1.7684714852626602, + 5.167567490120088, + 37.10022546727608, + -3.533827480127306, + ]; + assert_eq!(buf, expected); + for (got, exp) in buf.iter().zip(expected.iter()) { + println!("got = {}, exp = {}", got, exp); + assert!((got - exp).abs() < 1e-6); + } + } +} diff --git a/src/distr/mod.rs b/src/distr/mod.rs index 70ecf2410a..2351f639f4 100644 --- a/src/distr/mod.rs +++ b/src/distr/mod.rs @@ -87,6 +87,7 @@ //! [`statrs`]: https://crates.io/crates/statrs mod bernoulli; +mod cauchy; mod distribution; mod float; mod integer; @@ -103,6 +104,7 @@ pub mod uniform; pub mod weighted; pub use self::bernoulli::{Bernoulli, BernoulliError}; +pub use self::cauchy::{Cauchy, CauchyError}; #[cfg(feature = "alloc")] pub use self::distribution::SampleString; pub use self::distribution::{Distribution, Iter, Map}; From 14454332b96804b15724298d1143bff79603167f Mon Sep 17 00:00:00 2001 From: Tom Conway Date: Thu, 21 May 2026 13:40:14 +1000 Subject: [PATCH 2/3] convert a range into an iterator --- src/distr/cauchy.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/distr/cauchy.rs b/src/distr/cauchy.rs index 9fc974bbe5..469e7e733c 100644 --- a/src/distr/cauchy.rs +++ b/src/distr/cauchy.rs @@ -173,9 +173,8 @@ mod test { // Calculate the Kolmogorov-Smirnoff divergence between the // empirical and expected distributions, and check it's small. let mut d0 = 0.0; - for i in 0..N { + for (i, &x) in values.iter().enumerate() { let p = (i as f64) / (N as f64); - let x = values[i]; let q = x.atan() / std::f64::consts::PI + 0.5; let d = (p - q).abs(); if d > d0 { From 2843e30938d90af4c04284c20bce912becfb4042 Mon Sep 17 00:00:00 2001 From: Tom Conway Date: Thu, 21 May 2026 13:43:51 +1000 Subject: [PATCH 3/3] remove uses of std --- src/distr/cauchy.rs | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/distr/cauchy.rs b/src/distr/cauchy.rs index 469e7e733c..27b51d64ae 100644 --- a/src/distr/cauchy.rs +++ b/src/distr/cauchy.rs @@ -108,8 +108,6 @@ impl Distribution for Cauchy { #[cfg(test)] mod test { - use std::{println, vec::Vec}; - use super::*; use crate::RngExt; @@ -164,9 +162,9 @@ mod test { let distr = Cauchy::new(0.0, 1.0).unwrap(); const N: usize = 100_000; - let mut values: Vec = Vec::new(); - for _ in 0..N { - values.push(rng.sample(distr)); + let mut values: [f64; 100000] = [0.0f64; N]; + for x in values.iter_mut() { + *x = rng.sample(distr); } values.sort_by(|a, b| a.partial_cmp(b).unwrap()); @@ -175,7 +173,7 @@ mod test { let mut d0 = 0.0; for (i, &x) in values.iter().enumerate() { let p = (i as f64) / (N as f64); - let q = x.atan() / std::f64::consts::PI + 0.5; + let q = x.atan() / core::f64::consts::PI + 0.5; let d = (p - q).abs(); if d > d0 { d0 = d; @@ -207,7 +205,6 @@ mod test { ]; assert_eq!(buf, expected); for (got, exp) in buf.iter().zip(expected.iter()) { - println!("got = {}, exp = {}", got, exp); assert!((got - exp).abs() < 1e-6); } }