diff --git a/Cargo.toml b/Cargo.toml index 4b3d2a4b290..0a9e28cfe40 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -55,6 +55,10 @@ fuchsia-zircon = { version = "0.3.2", optional = true } stdweb = { version = "0.4", optional = true } wasm-bindgen = { version = "0.2.12", optional = true } +[dev-dependencies] +# This has a histogram implementation used for testing uniformity. +average = "0.9.2" + [build-dependencies] rustc_version = "0.2" diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs index 03d2d58e7bc..082d22278f4 100644 --- a/src/distributions/mod.rs +++ b/src/distributions/mod.rs @@ -100,6 +100,8 @@ //! - [`FisherF`] distribution //! - Multivariate probability distributions //! - [`Dirichlet`] distribution +//! - [`UnitSphereSurface`] distribution +//! - [`UnitCircle`] distribution //! //! # Examples //! @@ -169,6 +171,8 @@ //! [`Uniform`]: struct.Uniform.html //! [`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; @@ -178,6 +182,8 @@ pub use self::other::Alphanumeric; 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}; @@ -190,6 +196,8 @@ pub use self::bernoulli::Bernoulli; 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; @@ -578,7 +586,7 @@ mod tests { Weighted { weight: x, item: 2 }, Weighted { weight: 1, item: 3 }]); } - + #[cfg(feature="std")] #[test] fn test_distributions_iter() { diff --git a/src/distributions/unit_circle.rs b/src/distributions/unit_circle.rs new file mode 100644 index 00000000000..a1cf66ae6c0 --- /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 new file mode 100644 index 00000000000..039e1f6c7b2 --- /dev/null +++ b/src/distributions/unit_sphere.rs @@ -0,0 +1,92 @@ +use Rng; +use distributions::{Distribution, Uniform}; + +/// Samples uniformly from the surface of the unit sphere in three dimensions. +/// +/// Implemented via a method by Marsaglia[^1]. +/// +/// +/// # Example +/// +/// ``` +/// use rand::distributions::{UnitSphereSurface, Distribution}; +/// +/// let sphere = UnitSphereSurface::new(); +/// let v = sphere.sample(&mut rand::thread_rng()); +/// println!("{:?} is from the unit sphere surface.", v) +/// ``` +/// +/// [^1]: Marsaglia, George (1972). [*Choosing a Point from the Surface of a +/// Sphere.*](https://doi.org/10.1214/aoms/1177692644) +/// Ann. Math. Statist. 43, no. 2, 645--646. +#[derive(Clone, Copy, Debug)] +pub struct UnitSphereSurface { + uniform: Uniform, +} + +impl UnitSphereSurface { + /// Construct a new `UnitSphereSurface` distribution. + #[inline] + pub fn new() -> UnitSphereSurface { + UnitSphereSurface { uniform: Uniform::new(-1., 1.) } + } +} + +impl Distribution<[f64; 3]> for UnitSphereSurface { + #[inline] + fn sample(&self, rng: &mut R) -> [f64; 3] { + loop { + let (x1, x2) = (self.uniform.sample(rng), self.uniform.sample(rng)); + let sum = x1*x1 + x2*x2; + if sum >= 1. { + continue; + } + let factor = 2. * (1.0_f64 - sum).sqrt(); + return [x1 * factor, x2 * factor, 1. - 2.*sum]; + } + } +} + +#[cfg(test)] +mod tests { + use distributions::Distribution; + use super::UnitSphereSurface; + + /// 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 = UnitSphereSurface::new(); + for _ in 0..1000 { + let x = dist.sample(&mut rng); + assert_almost_eq!(x[0]*x[0] + x[1]*x[1] + x[2]*x[2], 1., 1e-15); + } + } + + #[test] + fn value_stability() { + let mut rng = ::test::rng(2); + let dist = UnitSphereSurface::new(); + assert_eq!(dist.sample(&mut rng), + [0.1660727273690104, 0.5202698793511903, 0.8376986939610902]); + assert_eq!(dist.sample(&mut rng), + [0.40052443371799834, 0.4237054996596154, -0.812436968356975]); + assert_eq!(dist.sample(&mut rng), + [0.9209910250970449, -0.32692477745072107, 0.21188610520628948]); + } +} diff --git a/tests/uniformity.rs b/tests/uniformity.rs new file mode 100644 index 00000000000..a39026bc253 --- /dev/null +++ b/tests/uniformity.rs @@ -0,0 +1,59 @@ +#![cfg(feature = "std")] + +#[macro_use] +extern crate average; +extern crate rand; + +use std as core; +use rand::FromEntropy; +use rand::distributions::Distribution; +use average::Histogram; + +const N_BINS: usize = 100; +const N_SAMPLES: u32 = 1_000_000; +const TOL: f64 = 1e-3; +define_histogram!(hist, 100); +use hist::Histogram as Histogram100; + +#[test] +fn unit_sphere() { + const N_DIM: usize = 3; + let h = Histogram100::with_const_width(-1., 1.); + let mut histograms = [h.clone(), h.clone(), h]; + let dist = rand::distributions::UnitSphereSurface::new(); + let mut rng = rand::rngs::SmallRng::from_entropy(); + for _ in 0..N_SAMPLES { + let v = dist.sample(&mut rng); + for i in 0..N_DIM { + histograms[i].add(v[i]).map_err( + |e| { println!("v: {}", v[i]); e } + ).unwrap(); + } + } + for h in &histograms { + let sum: u64 = h.bins().iter().sum(); + println!("{:?}", h); + for &b in h.bins() { + let p = (b as f64) / (sum as f64); + assert!((p - 1.0 / (N_BINS as f64)).abs() < TOL, "{}", p); + } + } +} + +#[test] +fn unit_circle() { + use ::std::f64::consts::PI; + let mut h = Histogram100::with_const_width(-PI, PI); + let dist = rand::distributions::UnitCircle::new(); + let mut rng = rand::rngs::SmallRng::from_entropy(); + for _ in 0..N_SAMPLES { + let v = dist.sample(&mut rng); + h.add(v[0].atan2(v[1])).unwrap(); + } + let sum: u64 = h.bins().iter().sum(); + println!("{:?}", h); + for &b in h.bins() { + let p = (b as f64) / (sum as f64); + assert!((p - 1.0 / (N_BINS as f64)).abs() < TOL, "{}", p); + } +}