Support for Dirichlet distribution

src/distributions/dirichlet.rs

@@ -0,0 +1,139 @@
+// Copyright 2013 The Rust Project Developers. See the COPYRIGHT
+// file at the top-level directory of this distribution and at
+// https://rust-lang.org/COPYRIGHT.
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+//! The dirichlet distribution.
+
+use Rng;
+use distributions::Distribution;
+use distributions::gamma::Gamma;
+
+/// The dirichelet distribution `Dirichlet(alpha)`.
+///
+/// The Dirichlet distribution } is a family of continuous multivariate probability distributions parameterized by
+/// a vector alpha of positive reals
+/// It is a multivariate generalization of the beta distribution.
+///
+/// # Example
+///
+/// ```
+/// use rand::prelude::*;
+/// use rand::distributions::Dirichlet;
+///
+/// let dirichlet = Dirichlet::new(&vec![1.0, 2.0, 3.0]);
+/// let samples = dirichlet.sample(&mut rand::thread_rng());
+/// println!("{:?} is from a Dirichlet([1.0, 2.0, 3.0]) distribution", samples);
+/// ```
+
+#[derive(Clone, Debug)]
+pub struct Dirichlet {
+    /// Concentration parameters (alpha)
+    alpha: Vec<f64>,
+}
+
+impl Dirichlet {
+    /// Construct a new `Dirichlet` with the given alpha parameter
+    /// `alpha`. Panics if `alpha.len() < 2`.
+    #[inline]
+    pub fn new(alpha: &[f64]) -> Dirichlet {
+        assert!(
+            alpha.len() > 1,
+            "Dirichlet::new called with `alpha` with length <  2"
+        );
+        for i in 0..alpha.len() {
+            assert!(
+                alpha[i] > 0.0,
+                "Dirichlet::new called with `alpha`  <=  0.0"
+            );
+        }
+
+        Dirichlet {
+            alpha: alpha.to_vec(),
+        }
+    }
+
+    /// Construct a new `Dirichlet` with the given shape parameter and size
+    /// `alpha`. Panics if `alpha <= 0.0`.
+    /// `size` . Panic if `size < 2`
+    #[inline]
+    pub fn new_with_param(alpha: f64, size: usize) -> Dirichlet {
+        assert!(alpha > 0.0, "Dirichlet::new called with `alpha`  <=  0.0");
+        assert!(size > 1, "Dirichlet::new called with `size`  <=  1");
+        Dirichlet {
+            alpha: vec![alpha; size],
+        }
+    }
+}
+
+impl Distribution<Vec<f64>> for Dirichlet {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec<f64> {
+        let n = self.alpha.len();
+        let mut samples = vec![0.0f64; n];
+        let mut sum = 0.0f64;
+
+        for i in 0..n {
+            let g = Gamma::new(self.alpha[i], 1.0);
+            samples[i] = g.sample(rng);
+            sum += samples[i];
+        }
+        let invacc = 1.0 / sum;
+        for i in 0..n {
+            samples[i] *= invacc;
+        }
+        samples
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::Dirichlet;
+    use distributions::Distribution;
+
+    #[test]
+    fn test_dirichlet() {
+        let d = Dirichlet::new(&vec![1.0, 2.0, 3.0]);
+        let mut rng = ::test::rng(221);
+        let samples = d.sample(&mut rng);
+        let _: Vec<f64> = samples
+            .into_iter()
+            .map(|x| {
+                assert!(x > 0.0);
+                x
+            })
+            .collect();
+    }
+
+    #[test]
+    fn test_dirichlet_with_param() {
+        let alpha = 0.5f64;
+        let size = 2;
+        let d = Dirichlet::new_with_param(alpha, size);
+        let mut rng = ::test::rng(221);
+        let samples = d.sample(&mut rng);
+        let _: Vec<f64> = samples
+            .into_iter()
+            .map(|x| {
+                assert!(x > 0.0);
+                x
+            })
+            .collect();
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_dirichlet_invalid_length() {
+        Dirichlet::new_with_param(0.5f64, 1);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_dirichlet_invalid_alpha() {
+        Dirichlet::new_with_param(0.0f64, 2);
+    }
+}

src/distributions/mod.rs

@@ -95,7 +95,7 @@
 //!   - [`ChiSquared`] distribution
 //!   - [`StudentT`] distribution
 //!   - [`FisherF`] distribution
-//!
+//! - Dirichlet distribution
 //!
 //! # Examples
 //!
@@ -148,6 +148,7 @@
 //! [`Bernoulli`]: struct.Bernoulli.html
 //! [`Binomial`]: struct.Binomial.html
 //! [`ChiSquared`]: struct.ChiSquared.html
+//! [`Dirichlet`]: struct.Dirichlet.html
 //! [`Exp`]: struct.Exp.html
 //! [`Exp1`]: struct.Exp1.html
 //! [`FisherF`]: struct.FisherF.html
@@ -180,6 +181,8 @@ pub use self::uniform::Uniform as Range;
 #[cfg(feature = "std")]
 #[doc(inline)] pub use self::binomial::Binomial;
 #[doc(inline)] pub use self::bernoulli::Bernoulli;
+#[cfg(feature = "std")]
+#[doc(inline)] pub use self::dirichlet::Dirichlet;
 
 pub mod uniform;
 #[cfg(feature="std")]
@@ -193,6 +196,8 @@ pub mod uniform;
 #[cfg(feature = "std")]
 #[doc(hidden)] pub mod binomial;
 #[doc(hidden)] pub mod bernoulli;
+#[cfg(feature = "std")]
+#[doc(hidden)] pub mod dirichlet;
 
 mod float;
 mod integer;
@@ -204,6 +209,7 @@ mod ziggurat_tables;
 #[cfg(feature="std")]
 use distributions::float::IntoFloat;
 
+
 /// Types that can be used to create a random instance of `Support`.
 #[deprecated(since="0.5.0", note="use Distribution instead")]
 pub trait Sample<Support> {