Implement sampling from the unit circle
This commit is contained in:
committed by
Diggory Hardy
parent
c7a76b487e
commit
96f43ed79c
@@ -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;
|
||||
|
||||
@@ -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<f64>,
|
||||
}
|
||||
|
||||
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<R: Rng + ?Sized>(&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]);
|
||||
}
|
||||
}
|
||||
@@ -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<f64>,
|
||||
|
||||
Reference in New Issue
Block a user