diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs index 621516e3..082d2227 100644 --- a/src/distributions/mod.rs +++ b/src/distributions/mod.rs @@ -101,6 +101,7 @@ //! - Multivariate probability distributions //! - [`Dirichlet`] distribution //! - [`UnitSphereSurface`] distribution +//! - [`UnitCircle`] distribution //! //! # Examples //! @@ -171,6 +172,7 @@ //! [`Uniform::new`]: struct.Uniform.html#method.new //! [`Uniform::new_inclusive`]: struct.Uniform.html#method.new_inclusive //! [`UnitSphereSurface`]: struct.UnitSphereSurface.html +//! [`UnitCircle`]: struct.UnitCircle.html //! [`WeightedIndex`]: struct.WeightedIndex.html use Rng; @@ -181,6 +183,7 @@ pub use self::float::{OpenClosed01, Open01}; pub use self::bernoulli::Bernoulli; #[cfg(feature="alloc")] pub use self::weighted::{WeightedIndex, WeightedError}; #[cfg(feature="std")] pub use self::unit_sphere::UnitSphereSurface; +#[cfg(feature="std")] pub use self::unit_circle::UnitCircle; #[cfg(feature="std")] pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT}; #[cfg(feature="std")] pub use self::normal::{Normal, LogNormal, StandardNormal}; #[cfg(feature="std")] pub use self::exponential::{Exp, Exp1}; @@ -194,6 +197,7 @@ pub mod uniform; mod bernoulli; #[cfg(feature="alloc")] mod weighted; #[cfg(feature="std")] mod unit_sphere; +#[cfg(feature="std")] mod unit_circle; #[cfg(feature="std")] mod gamma; #[cfg(feature="std")] mod normal; #[cfg(feature="std")] mod exponential; diff --git a/src/distributions/unit_circle.rs b/src/distributions/unit_circle.rs new file mode 100644 index 00000000..a1cf66ae --- /dev/null +++ b/src/distributions/unit_circle.rs @@ -0,0 +1,94 @@ +use Rng; +use distributions::{Distribution, Uniform}; + +/// Samples uniformly from the edge of the unit circle in two dimensions. +/// +/// Implemented via a method by von Neumann[^1]. +/// +/// +/// # Example +/// +/// ``` +/// use rand::distributions::{UnitCircle, Distribution}; +/// +/// let circle = UnitCircle::new(); +/// let v = circle.sample(&mut rand::thread_rng()); +/// println!("{:?} is from the unit circle.", v) +/// ``` +/// +/// [^1]: von Neumann, J. (1951) [*Various Techniques Used in Connection with +/// Random Digits.*](https://mcnp.lanl.gov/pdf_files/nbs_vonneumann.pdf) +/// NBS Appl. Math. Ser., No. 12. Washington, DC: U.S. Government Printing +/// Office, pp. 36-38. +#[derive(Clone, Copy, Debug)] +pub struct UnitCircle { + uniform: Uniform, +} + +impl UnitCircle { + /// Construct a new `UnitCircle` distribution. + #[inline] + pub fn new() -> UnitCircle { + UnitCircle { uniform: Uniform::new(-1., 1.) } + } +} + +impl Distribution<[f64; 2]> for UnitCircle { + #[inline] + fn sample(&self, rng: &mut R) -> [f64; 2] { + let mut x1; + let mut x2; + let mut sum; + loop { + x1 = self.uniform.sample(rng); + x2 = self.uniform.sample(rng); + sum = x1*x1 + x2*x2; + if sum < 1. { + break; + } + } + let diff = x1*x1 - x2*x2; + [diff / sum, 2.*x1*x2 / sum] + } +} + +#[cfg(test)] +mod tests { + use distributions::Distribution; + use super::UnitCircle; + + /// Assert that two numbers are almost equal to each other. + /// + /// On panic, this macro will print the values of the expressions with their + /// debug representations. + macro_rules! assert_almost_eq { + ($a:expr, $b:expr, $prec:expr) => ( + let diff = ($a - $b).abs(); + if diff > $prec { + panic!(format!( + "assertion failed: `abs(left - right) = {:.1e} < {:e}`, \ + (left: `{}`, right: `{}`)", + diff, $prec, $a, $b)); + } + ); + } + + #[test] + fn norm() { + let mut rng = ::test::rng(1); + let dist = UnitCircle::new(); + for _ in 0..1000 { + let x = dist.sample(&mut rng); + assert_almost_eq!(x[0]*x[0] + x[1]*x[1], 1., 1e-15); + } + } + + #[test] + fn value_stability() { + let mut rng = ::test::rng(2); + let dist = UnitCircle::new(); + assert_eq!(dist.sample(&mut rng), [-0.8150602311723979, 0.5793762331690843]); + assert_eq!(dist.sample(&mut rng), [-0.056204569973983196, 0.998419273809375]); + assert_eq!(dist.sample(&mut rng), [0.7761923749562624, -0.630496151502733]); + } +} diff --git a/src/distributions/unit_sphere.rs b/src/distributions/unit_sphere.rs index e0138e1d..039e1f6c 100644 --- a/src/distributions/unit_sphere.rs +++ b/src/distributions/unit_sphere.rs @@ -18,7 +18,7 @@ use distributions::{Distribution, Uniform}; /// /// [^1]: Marsaglia, George (1972). [*Choosing a Point from the Surface of a /// Sphere.*](https://doi.org/10.1214/aoms/1177692644) -/// Ann. Math. Statist. 43 (1972), no. 2, 645--646. +/// Ann. Math. Statist. 43, no. 2, 645--646. #[derive(Clone, Copy, Debug)] pub struct UnitSphereSurface { uniform: Uniform,