From ef38c4508f8850d89d200d9b1ea2ccbb94beab5e Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sun, 25 Nov 2012 14:51:25 +0100
Subject: [PATCH] Distributions: consistent sample methods for discrete
 distributions #32

---
 .../Distributions/Discrete/Bernoulli.cs       | 113 ++++-------
 .../Distributions/Discrete/Binomial.cs        | 130 +++++-------
 .../Distributions/Discrete/Categorical.cs     | 153 ++++++--------
 .../Discrete/ConwayMaxwellPoisson.cs          | 171 ++++++++--------
 .../Distributions/Discrete/DiscreteUniform.cs | 118 ++++-------
 .../Distributions/Discrete/Geometric.cs       | 119 ++++++-----
 .../Distributions/Discrete/Hypergeometric.cs  | 115 ++++++-----
 .../Discrete/NegativeBinomial.cs              | 164 +++++++--------
 .../Distributions/Discrete/Poisson.cs         | 186 +++++++-----------
 src/Numerics/Distributions/Discrete/Zipf.cs   | 142 ++++++-------
 .../Distributions/Multivariate/Multinomial.cs |   4 +-
 11 files changed, 620 insertions(+), 795 deletions(-)

diff --git a/src/Numerics/Distributions/Discrete/Bernoulli.cs b/src/Numerics/Distributions/Discrete/Bernoulli.cs
index 15f86424..941c9f80 100644
--- a/src/Numerics/Distributions/Discrete/Bernoulli.cs
+++ b/src/Numerics/Distributions/Discrete/Bernoulli.cs
@@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The probability of generating a one.
         /// </summary>
-        private double _p;
+        double _p;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the Bernoulli class.
@@ -77,7 +77,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="p">The probability of generating a one.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double p)
+        static bool IsValidParameterSet(double p)
         {
             if (p >= 0.0 && p <= 1.0)
             {
@@ -92,7 +92,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="p">The probability of generating a one.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double p)
+        void SetParameters(double p)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
             {
@@ -107,15 +107,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double P
         {
-            get
-            {
-                return _p;
-            }
+            get { return _p; }
 
-            set
-            {
-                SetParameters(value);
-            }
+            set { SetParameters(value); }
         }
 
         #region IDistribution Members
@@ -125,10 +119,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -146,10 +137,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return _p;
-            }
+            get { return _p; }
         }
 
         /// <summary>
@@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(_p * (1.0 - _p));
-            }
+            get { return Math.Sqrt(_p * (1.0 - _p)); }
         }
 
         /// <summary>
@@ -168,10 +153,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _p * (1.0 - _p);
-            }
+            get { return _p * (1.0 - _p); }
         }
 
         /// <summary>
@@ -179,10 +161,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p));
-            }
+            get { return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p)); }
         }
 
         /// <summary>
@@ -190,10 +169,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p));
-            }
+            get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p)); }
         }
 
         /// <summary>
@@ -201,10 +177,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -212,10 +185,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return 1;
-            }
+            get { return 1; }
         }
 
         /// <summary>
@@ -229,7 +199,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (x < 1.0)
             {
                 return 1.0 - _p;
@@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return _p > 0.5 ? 1 : 0;
-            }
+            get { return _p > 0.5 ? 1 : 0; }
         }
 
         /// <summary>
@@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                throw new NotSupportedException("The median of the Bernoulli distribution is undefined.");
-            }
+            get { throw new NotSupportedException("The median of the Bernoulli distribution is undefined."); }
         }
 
         /// <summary>
@@ -299,13 +263,31 @@ namespace MathNet.Numerics.Distributions
             return k == 1 ? Math.Log(_p) : Double.NegativeInfinity;
         }
 
+        #endregion
+
+        /// <summary>
+        /// Generates one sample from the Bernoulli distribution.
+        /// </summary>
+        /// <param name="rnd">The random source to use.</param>
+        /// <param name="p">The probability of generating a one.</param>
+        /// <returns>A random sample from the Bernoulli distribution.</returns>
+        internal static int SampleUnchecked(Random rnd, double p)
+        {
+            if (rnd.NextDouble() < p)
+            {
+                return 1;
+            }
+
+            return 0;
+        }
+
         /// <summary>
         /// Samples a Bernoulli distributed random variable.
         /// </summary>
         /// <returns>A sample from the Bernoulli distribution.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _p);
+            return SampleUnchecked(RandomSource, _p);
         }
 
         /// <summary>
@@ -316,12 +298,10 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _p);
+                yield return SampleUnchecked(RandomSource, _p);
             }
         }
 
-        #endregion
-
         /// <summary>
         /// Samples a Bernoulli distributed random variable.
         /// </summary>
@@ -335,7 +315,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
             }
 
-            return DoSample(rnd, p);
+            return SampleUnchecked(rnd, p);
         }
 
         /// <summary>
@@ -353,24 +333,9 @@ namespace MathNet.Numerics.Distributions
 
             while (true)
             {
-                yield return DoSample(rnd, p);
+                yield return SampleUnchecked(rnd, p);
             }
         }
 
-        /// <summary>
-        /// Generates one sample from the Bernoulli distribution.
-        /// </summary>
-        /// <param name="rnd">The random source to use.</param>
-        /// <param name="p">The probability of generating a one.</param>
-        /// <returns>A random sample from the Bernoulli distribution.</returns>
-        private static int DoSample(Random rnd, double p)
-        {
-            if (rnd.NextDouble() < p)
-            {
-                return 1;
-            }
-
-            return 0;
-        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Binomial.cs b/src/Numerics/Distributions/Discrete/Binomial.cs
index 5f836591..f37a5881 100644
--- a/src/Numerics/Distributions/Discrete/Binomial.cs
+++ b/src/Numerics/Distributions/Discrete/Binomial.cs
@@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// Stores the normalized binomial probability.
         /// </summary>
-        private double _p;
+        double _p;
 
         /// <summary>
         /// The number of trials.
         /// </summary>
-        private int _n;
+        int _n;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the Binomial class.
@@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="p">The success probability of a trial.</param>
         /// <param name="n">The number of trials.</param>
         /// <returns><c>false</c> <paramref name="p"/> is not in the interval [0.0,1.0] or <paramref name="n"/> is negative, <c>true</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double p, int n)
+        static bool IsValidParameterSet(double p, int n)
         {
             if (p < 0.0 || p > 1.0 || Double.IsNaN(p))
             {
@@ -107,7 +107,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="n">The number of trials.</param>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="p"/> is not in the interval [0.0,1.0].</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="n"/> is negative.</exception>
-        private void SetParameters(double p, int n)
+        void SetParameters(double p, int n)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
             {
@@ -123,15 +123,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double P
         {
-            get
-            {
-                return _p;
-            }
+            get { return _p; }
 
-            set
-            {
-                SetParameters(value, _n);
-            }
+            set { SetParameters(value, _n); }
         }
 
         /// <summary>
@@ -139,15 +133,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int N
         {
-            get
-            {
-                return _n;
-            }
+            get { return _n; }
 
-            set
-            {
-                SetParameters(_p, value);
-            }
+            set { SetParameters(_p, value); }
         }
 
         #region IDistribution Members
@@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return _p * _n;
-            }
+            get { return _p * _n; }
         }
 
         /// <summary>
@@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(_p * (1.0 - _p) * _n);
-            }
+            get { return Math.Sqrt(_p * (1.0 - _p) * _n); }
         }
 
         /// <summary>
@@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _p * (1.0 - _p) * _n;
-            }
+            get { return _p * (1.0 - _p) * _n; }
         }
 
         /// <summary>
@@ -234,10 +210,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p));
-            }
+            get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p)); }
         }
 
         /// <summary>
@@ -245,10 +218,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -256,10 +226,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return _n;
-            }
+            get { return _n; }
         }
 
         /// <summary>
@@ -273,7 +240,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (x > _n)
             {
                 return 1.0;
@@ -303,7 +270,7 @@ namespace MathNet.Numerics.Distributions
                 {
                     return _n;
                 }
-                
+
                 if (_p == 0.0)
                 {
                     return 0;
@@ -318,10 +285,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                return (int)Math.Floor(_p * _n);
-            }
+            get { return (int)Math.Floor(_p * _n); }
         }
 
         /// <summary>
@@ -345,7 +309,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 1.0;
             }
-            
+
             if (_p == 0.0)
             {
                 return 0.0;
@@ -355,7 +319,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 1.0;
             }
-            
+
             if (_p == 1.0)
             {
                 return 0.0;
@@ -385,7 +349,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (_p == 0.0)
             {
                 return Double.NegativeInfinity;
@@ -395,7 +359,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (_p == 1.0)
             {
                 return Double.NegativeInfinity;
@@ -404,13 +368,33 @@ namespace MathNet.Numerics.Distributions
             return SpecialFunctions.BinomialLn(_n, k) + (k * Math.Log(_p)) + ((_n - k) * Math.Log(1.0 - _p));
         }
 
+        #endregion
+
+        /// <summary>
+        /// Generates a sample from the Binomial distribution without doing parameter checking.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="p">The success probability of a trial; must be in the interval [0.0, 1.0].</param>
+        /// <param name="n">The number of trials; must be positive.</param>
+        /// <returns>The number of successful trials.</returns>
+        internal static int SampleUnchecked(Random rnd, double p, int n)
+        {
+            var k = 0;
+            for (var i = 0; i < n; i++)
+            {
+                k += rnd.NextDouble() < p ? 1 : 0;
+            }
+
+            return k;
+        }
+
         /// <summary>
         /// Samples a Binomially distributed random variable.
         /// </summary>
         /// <returns>The number of successes in N trials.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _p, _n);
+            return SampleUnchecked(RandomSource, _p, _n);
         }
 
         /// <summary>
@@ -421,12 +405,10 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _p, _n);
+                yield return SampleUnchecked(RandomSource, _p, _n);
             }
         }
 
-        #endregion
-
         /// <summary>
         /// Samples a binomially distributed random variable.
         /// </summary>
@@ -441,7 +423,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
             }
 
-            return DoSample(rnd, p, n);
+            return SampleUnchecked(rnd, p, n);
         }
 
         /// <summary>
@@ -460,26 +442,8 @@ namespace MathNet.Numerics.Distributions
 
             while (true)
             {
-                yield return DoSample(rnd, p, n);
-            }
-        }
-
-        /// <summary>
-        /// Generates a sample from the Binomial distribution without doing parameter checking.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="p">The success probability of a trial; must be in the interval [0.0, 1.0].</param>
-        /// <param name="n">The number of trials; must be positive.</param>
-        /// <returns>The number of successful trials.</returns>
-        private static int DoSample(Random rnd, double p, int n)
-        {
-            var k = 0;
-            for (var i = 0; i < n; i++)
-            {
-                k += rnd.NextDouble() < p ? 1 : 0;
+                yield return SampleUnchecked(rnd, p, n);
             }
-
-            return k;
         }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Categorical.cs b/src/Numerics/Distributions/Discrete/Categorical.cs
index 7f8674e1..97e51e86 100644
--- a/src/Numerics/Distributions/Discrete/Categorical.cs
+++ b/src/Numerics/Distributions/Discrete/Categorical.cs
@@ -50,12 +50,12 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// Stores the unnormalized categorical probabilities.
         /// </summary>
-        private double[] _p;
+        double[] _p;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the Categorical class.
@@ -110,7 +110,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
         /// as this is often impossible using floating point arithmetic.</param>
         /// <returns>If any of the probabilities are negative returns <c>false</c>, or if the sum of parameters is 0.0; otherwise <c>true</c></returns>
-        private static bool IsValidParameterSet(IEnumerable<double> p)
+        static bool IsValidParameterSet(IEnumerable<double> p)
         {
             var sum = 0.0;
             foreach (double t in p)
@@ -119,7 +119,7 @@ namespace MathNet.Numerics.Distributions
                 {
                     return false;
                 }
-                
+
                 sum += t;
             }
 
@@ -132,7 +132,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
         /// as this is often impossible using floating point arithmetic.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double[] p)
+        void SetParameters(double[] p)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
             {
@@ -163,10 +163,7 @@ namespace MathNet.Numerics.Distributions
                 return p;
             }
 
-            set
-            {
-                SetParameters(value);
-            }
+            set { SetParameters(value); }
         }
 
         #region IDistribution Members
@@ -176,10 +173,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -197,10 +191,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return _p.Mean();
-            }
+            get { return _p.Mean(); }
         }
 
         /// <summary>
@@ -208,10 +199,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return _p.StandardDeviation();
-            }
+            get { return _p.StandardDeviation(); }
         }
 
         /// <summary>
@@ -219,10 +207,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _p.Variance();
-            }
+            get { return _p.Variance(); }
         }
 
         /// <summary>
@@ -230,10 +215,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                return _p.Sum(p => p * Math.Log(p));
-            }
+            get { return _p.Sum(p => p * Math.Log(p)); }
         }
 
         /// <summary>
@@ -242,10 +224,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks>
         public double Skewness
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -253,10 +232,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -264,10 +240,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return _p.Length - 1;
-            }
+            get { return _p.Length - 1; }
         }
 
         /// <summary>
@@ -281,7 +254,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (x >= _p.Length)
             {
                 return 1.0;
@@ -301,10 +274,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks>
         public int Mode
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -312,10 +282,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                return (int)_p.Median();
-            }
+            get { return (int)_p.Median(); }
         }
 
         /// <summary>
@@ -358,6 +325,47 @@ namespace MathNet.Numerics.Distributions
             return Math.Log(_p[k]);
         }
 
+        #endregion
+
+        /// <summary>
+        /// Computes the unnormalized cumulative distribution function. This method performs no
+        /// parameter checking.
+        /// </summary>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <returns>An array representing the unnormalized cumulative distribution function.</returns>
+        internal static double[] UnnormalizedCdf(double[] p)
+        {
+            var cp = (double[])p.Clone();
+
+            for (var i = 1; i < p.Length; i++)
+            {
+                cp[i] += cp[i - 1];
+            }
+
+            return cp;
+        }
+
+        /// <summary>
+        /// Returns one trials from the categorical distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="cdf">The cumulative distribution of the probability distribution.</param>
+        /// <returns>One sample from the categorical distribution implied by <paramref name="cdf"/>.</returns>
+        internal static int SampleUnchecked(Random rnd, double[] cdf)
+        {
+            // TODO : use binary search to speed up this procedure.
+            var u = rnd.NextDouble() * cdf[cdf.Length - 1];
+
+            var idx = 0;
+            while (u > cdf[idx])
+            {
+                idx++;
+            }
+
+            return idx;
+        }
+
         /// <summary>
         /// Samples a Binomially distributed random variable.
         /// </summary>
@@ -376,8 +384,6 @@ namespace MathNet.Numerics.Distributions
             return Samples(RandomSource, _p);
         }
 
-        #endregion
-
         /// <summary>
         /// Samples one categorical distributed random variable; also known as the Discrete distribution.
         /// </summary>
@@ -395,7 +401,7 @@ namespace MathNet.Numerics.Distributions
             // The cumulative density of p.
             var cp = UnnormalizedCdf(p);
 
-            return DoSample(rnd, cp);
+            return SampleUnchecked(rnd, cp);
         }
 
         /// <summary>
@@ -417,47 +423,8 @@ namespace MathNet.Numerics.Distributions
 
             while (true)
             {
-                yield return DoSample(rnd, cp);
+                yield return SampleUnchecked(rnd, cp);
             }
         }
-
-        /// <summary>
-        /// Computes the unnormalized cumulative distribution function. This method performs no
-        /// parameter checking.
-        /// </summary>
-        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
-        /// as this is often impossible using floating point arithmetic.</param>
-        /// <returns>An array representing the unnormalized cumulative distribution function.</returns>
-        internal static double[] UnnormalizedCdf(double[] p)
-        {
-            var cp = (double[])p.Clone();
-
-            for (var i = 1; i < p.Length; i++)
-            {
-                cp[i] += cp[i - 1];
-            }
-
-            return cp;
-        }
-
-        /// <summary>
-        /// Returns one trials from the categorical distribution.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="cdf">The cumulative distribution of the probability distribution.</param>
-        /// <returns>One sample from the categorical distribution implied by <paramref name="cdf"/>.</returns>
-        internal static int DoSample(Random rnd, double[] cdf)
-        {
-            // TODO : use binary search to speed up this procedure.
-            var u = rnd.NextDouble() * cdf[cdf.Length - 1];
-            
-            var idx = 0;
-            while (u > cdf[idx])
-            {
-                idx++;
-            }
-
-            return idx;
-        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs b/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs
index 719d8455..9dcdf47f 100644
--- a/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs
+++ b/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs
@@ -52,37 +52,37 @@ namespace MathNet.Numerics.Distributions
         /// Since many properties of the distribution can only be computed approximately, the tolerance
         /// level specifies how much error we accept.
         /// </summary>
-        private const double Tolerance = 1e-12;
+        const double Tolerance = 1e-12;
 
         /// <summary>
         /// The mean of the distribution.
         /// </summary>
-        private double _mean = double.MinValue;
+        double _mean = double.MinValue;
 
         /// <summary>
         ///  The variance of the distribution.
         /// </summary>
-        private double _variance = double.MinValue;
+        double _variance = double.MinValue;
 
         /// <summary>
         /// Caches the value of the normalization constant.
         /// </summary>
-        private double _z = double.MinValue;
+        double _z = double.MinValue;
 
         /// <summary>
         /// The lambda parameter.
         /// </summary>
-        private double _lambda;
+        double _lambda;
 
         /// <summary>
         /// The nu parameter.
         /// </summary>
-        private double _nu;
+        double _nu;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the <see cref="ConwayMaxwellPoisson"/> class. 
@@ -105,7 +105,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="lambda">The lambda parameter.</param>
         /// <param name="nu">The nu parameter.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double lambda, double nu)
+        void SetParameters(double lambda, double nu)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
             {
@@ -122,7 +122,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="lambda">The lambda parameter.</param>
         /// <param name="nu">The nu parameter.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double lambda, double nu)
+        static bool IsValidParameterSet(double lambda, double nu)
         {
             if (lambda <= 0.0)
             {
@@ -143,15 +143,9 @@ namespace MathNet.Numerics.Distributions
         /// <value>The value of the lambda parameter.</value>
         public double Lambda
         {
-            get
-            {
-                return _lambda;
-            }
+            get { return _lambda; }
 
-            set
-            {
-                SetParameters(value, _nu);
-            }
+            set { SetParameters(value, _nu); }
         }
 
         /// <summary>
@@ -160,15 +154,9 @@ namespace MathNet.Numerics.Distributions
         /// <value>The value of the Nu parameter.</value>
         public double Nu
         {
-            get
-            {
-                return _nu;
-            }
+            get { return _nu; }
 
-            set
-            {
-                SetParameters(_lambda, value);
-            }
+            set { SetParameters(_lambda, value); }
         }
 
         /// <summary>
@@ -189,10 +177,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -333,10 +318,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(Variance);
-            }
+            get { return Math.Sqrt(Variance); }
         }
 
         /// <summary>
@@ -344,10 +326,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -355,10 +334,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -386,10 +362,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -397,10 +370,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -408,10 +378,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -419,10 +386,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -449,35 +413,12 @@ namespace MathNet.Numerics.Distributions
             return Math.Log(Probability(k));
         }
 
-        /// <summary>
-        /// Samples a Conway-Maxwell-Poisson distributed random variable.
-        /// </summary>
-        /// <returns>a sample from the distribution.</returns>
-        public int Sample()
-        {
-            return DoSample(RandomSource, _lambda, _nu, Z);
-        }
-
-        /// <summary>
-        /// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables.
-        /// </summary>
-        /// <returns>
-        /// a sequence of samples from a Conway-Maxwell-Poisson distribution.
-        /// </returns>
-        public IEnumerable<int> Samples()
-        {
-            while (true)
-            {
-                yield return DoSample(RandomSource, _lambda, _nu, Z);
-            }
-        }
-
         #endregion
 
         /// <summary>
         /// Gets the normalization constant of the Conway-Maxwell-Poisson distribution.
         /// </summary>
-        private double Z
+        double Z
         {
             get
             {
@@ -488,7 +429,7 @@ namespace MathNet.Numerics.Distributions
 
                 _z = Normalization(_lambda, _nu);
                 return _z;
-            } 
+            }
         }
 
         /// <summary>
@@ -499,7 +440,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>
         /// an approximate normalization constant for the CMP distribution.
         /// </returns>
-        private static double Normalization(double lambda, double nu)
+        static double Normalization(double lambda, double nu)
         {
             // Initialize Z with the first two terms.
             var z = 1.0 + lambda;
@@ -515,7 +456,7 @@ namespace MathNet.Numerics.Distributions
 
                 // The new term.
                 t = t * e;
-    
+
                 // The updated normalization constant.
                 z = z + t;
 
@@ -542,7 +483,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>
         /// One sample from the distribution implied by <paramref name="lambda"/>, <paramref name="nu"/>, and <paramref name="z"/>.
         /// </returns>
-        private static int DoSample(Random rnd, double lambda, double nu, double z)
+        internal static int SampleUnchecked(Random rnd, double lambda, double nu, double z)
         {
             var u = rnd.NextDouble();
             var p = 1.0 / z;
@@ -558,5 +499,65 @@ namespace MathNet.Numerics.Distributions
 
             return i;
         }
+
+        /// <summary>
+        /// Samples a Conway-Maxwell-Poisson distributed random variable.
+        /// </summary>
+        /// <returns>a sample from the distribution.</returns>
+        public int Sample()
+        {
+            return SampleUnchecked(RandomSource, _lambda, _nu, Z);
+        }
+
+        /// <summary>
+        /// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables.
+        /// </summary>
+        /// <returns>
+        /// a sequence of samples from a Conway-Maxwell-Poisson distribution.
+        /// </returns>
+        public IEnumerable<int> Samples()
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(RandomSource, _lambda, _nu, Z);
+            }
+        }
+
+        /// <summary>
+        /// Samples a random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lambda">The lambda parameter</param>
+        /// <param name="nu">The nu parameter.</param>
+        public static int Sample(Random rnd, double lambda, double nu)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            var z = Normalization(lambda, nu);
+            return SampleUnchecked(rnd, lambda, nu, z);
+        }
+
+        /// <summary>
+        /// Samples a sequence of this random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lambda">The lambda parameter</param>
+        /// <param name="nu">The nu parameter.</param>
+        public static IEnumerable<int> Samples(Random rnd, double lambda, double nu)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            var z = Normalization(lambda, nu);
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, lambda, nu, z);
+            }
+        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/DiscreteUniform.cs b/src/Numerics/Distributions/Discrete/DiscreteUniform.cs
index e096f2d9..c4d20347 100644
--- a/src/Numerics/Distributions/Discrete/DiscreteUniform.cs
+++ b/src/Numerics/Distributions/Discrete/DiscreteUniform.cs
@@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The distribution's lower bound.
         /// </summary>
-        private int _lower;
+        int _lower;
 
         /// <summary>
         /// The distribution's upper bound.
         /// </summary>
-        private int _upper;
+        int _upper;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the DiscreteUniform class.
@@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="lower">Lower bound.</param>
         /// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(int lower, int upper)
+        static bool IsValidParameterSet(int lower, int upper)
         {
             if (lower <= upper)
             {
@@ -101,7 +101,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="lower">Lower bound.</param>
         /// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(int lower, int upper)
+        void SetParameters(int lower, int upper)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(lower, upper))
             {
@@ -117,15 +117,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int LowerBound
         {
-            get
-            {
-                return _lower;
-            }
+            get { return _lower; }
 
-            set
-            {
-                SetParameters(value, _upper);
-            }
+            set { SetParameters(value, _upper); }
         }
 
         /// <summary>
@@ -133,15 +127,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int UpperBound
         {
-            get
-            {
-                return _upper;
-            }
+            get { return _upper; }
 
-            set
-            {
-                SetParameters(_lower, value);
-            }
+            set { SetParameters(_lower, value); }
         }
 
         #region IDistribution Members
@@ -151,10 +139,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -172,10 +157,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return (_lower + _upper) / 2.0;
-            }
+            get { return (_lower + _upper) / 2.0; }
         }
 
         /// <summary>
@@ -183,10 +165,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0);
-            }
+            get { return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0); }
         }
 
         /// <summary>
@@ -194,10 +173,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0;
-            }
+            get { return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0; }
         }
 
         /// <summary>
@@ -205,10 +181,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                return Math.Log(_upper - _lower + 1.0);
-            }
+            get { return Math.Log(_upper - _lower + 1.0); }
         }
 
         /// <summary>
@@ -216,10 +189,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return 0.0;
-            }
+            get { return 0.0; }
         }
 
         /// <summary>
@@ -227,10 +197,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return _lower;
-            }
+            get { return _lower; }
         }
 
         /// <summary>
@@ -238,10 +205,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return _upper;
-            }
+            get { return _upper; }
         }
 
         /// <summary>
@@ -255,7 +219,7 @@ namespace MathNet.Numerics.Distributions
             {
                 return 0.0;
             }
-            
+
             if (x >= _upper)
             {
                 return 1.0;
@@ -273,10 +237,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return (int)Math.Floor((_lower + _upper) / 2.0);
-            }
+            get { return (int)Math.Floor((_lower + _upper) / 2.0); }
         }
 
         /// <summary>
@@ -284,10 +245,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                return (int)Math.Floor((_lower + _upper) / 2.0);
-            }
+            get { return (int)Math.Floor((_lower + _upper) / 2.0); }
         }
 
         /// <summary>
@@ -324,13 +282,27 @@ namespace MathNet.Numerics.Distributions
             return Double.NegativeInfinity;
         }
 
+        #endregion
+
+        /// <summary>
+        /// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking.
+        /// </summary>
+        /// <param name="rnd">The random source to use.</param>
+        /// <param name="lower">The lower bound of the uniform random variable.</param>
+        /// <param name="upper">The upper bound of the uniform random variable.</param>
+        /// <returns>A random sample from the discrete uniform distribution.</returns>
+        internal static int SampleUnchecked(Random rnd, int lower, int upper)
+        {
+            return (rnd.Next() % (upper - lower + 1)) + lower;
+        }
+
         /// <summary>
         /// Draws a random sample from the distribution.
         /// </summary>
         /// <returns>a sample from the distribution.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _lower, _upper);
+            return SampleUnchecked(RandomSource, _lower, _upper);
         }
 
         /// <summary>
@@ -341,12 +313,10 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _lower, _upper);
+                yield return SampleUnchecked(RandomSource, _lower, _upper);
             }
         }
 
-        #endregion
-
         /// <summary>
         /// Samples a uniformly distributed random variable.
         /// </summary>
@@ -361,7 +331,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
             }
 
-            return DoSample(rnd, lower, upper);
+            return SampleUnchecked(rnd, lower, upper);
         }
 
         /// <summary>
@@ -380,20 +350,8 @@ namespace MathNet.Numerics.Distributions
 
             while (true)
             {
-                yield return DoSample(rnd, lower, upper);
+                yield return SampleUnchecked(rnd, lower, upper);
             }
         }
-
-        /// <summary>
-        /// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking.
-        /// </summary>
-        /// <param name="rnd">The random source to use.</param>
-        /// <param name="lower">The lower bound of the uniform random variable.</param>
-        /// <param name="upper">The upper bound of the uniform random variable.</param>
-        /// <returns>A random sample from the discrete uniform distribution.</returns>
-        private static int DoSample(Random rnd, int lower, int upper)
-        {
-            return (rnd.Next() % (upper - lower + 1)) + lower;
-        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Geometric.cs b/src/Numerics/Distributions/Discrete/Geometric.cs
index 93b67cac..a0f01b5a 100644
--- a/src/Numerics/Distributions/Discrete/Geometric.cs
+++ b/src/Numerics/Distributions/Discrete/Geometric.cs
@@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The geometric distribution parameter. 
         /// </summary>
-        private double _p;
+        double _p;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the Geometric class.
@@ -68,7 +68,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="p">The probability of generating a one.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double p)
+        void SetParameters(double p)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
             {
@@ -83,7 +83,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="p">The probability of generating a one.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double p)
+        static bool IsValidParameterSet(double p)
         {
             if (p >= 0.0 && p <= 1.0)
             {
@@ -98,15 +98,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double P
         {
-            get
-            {
-                return _p;
-            }
+            get { return _p; }
 
-            set
-            {
-                SetParameters(value);
-            }
+            set { SetParameters(value); }
         }
 
         /// <summary>
@@ -127,10 +121,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -148,10 +139,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return 1.0 / _p;
-            }
+            get { return 1.0 / _p; }
         }
 
         /// <summary>
@@ -159,10 +147,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return (1.0 - _p) / (_p * _p);
-            }
+            get { return (1.0 - _p) / (_p * _p); }
         }
 
         /// <summary>
@@ -170,10 +155,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(1.0 - _p) / _p;
-            }
+            get { return Math.Sqrt(1.0 - _p) / _p; }
         }
 
         /// <summary>
@@ -181,10 +163,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p;
-            }
+            get { return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p; }
         }
 
         /// <summary>
@@ -193,10 +172,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>Throws a not supported exception.</remarks>
         public double Skewness
         {
-            get
-            {
-                return (2.0 - _p) / Math.Sqrt(1.0 - _p);
-            }
+            get { return (2.0 - _p) / Math.Sqrt(1.0 - _p); }
         }
 
         /// <summary>
@@ -218,10 +194,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return 1;
-            }
+            get { return 1; }
         }
 
         /// <summary>
@@ -229,10 +202,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p));
-            }
+            get { return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p)); }
         }
 
         /// <summary>
@@ -240,10 +210,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 1;
-            }
+            get { return 1; }
         }
 
         /// <summary>
@@ -251,10 +218,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return int.MaxValue;
-            }
+            get { return int.MaxValue; }
         }
 
         /// <summary>
@@ -291,13 +255,28 @@ namespace MathNet.Numerics.Distributions
             return ((k - 1) * Math.Log(1.0 - _p)) + Math.Log(_p);
         }
 
+        #endregion
+
+        /// <summary>
+        /// Returns one sample from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="p">The p parameter</param>
+        /// <returns>
+        /// One sample from the distribution implied by <paramref name="p"/>.
+        /// </returns>
+        internal static int SampleUnchecked(Random rnd, double p)
+        {
+            return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p));
+        }
+
         /// <summary>
         /// Samples a Geometric distributed random variable.
         /// </summary>
         /// <returns>A sample from the Geometric distribution.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _p);
+            return SampleUnchecked(RandomSource, _p);
         }
 
         /// <summary>
@@ -308,23 +287,41 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _p);
+                yield return SampleUnchecked(RandomSource, _p);
             }
         }
 
-        #endregion
+        /// <summary>
+        /// Samples a random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="p">The p parameter</param>
+        public static int Sample(Random rnd, double p)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            return SampleUnchecked(rnd, p);
+        }
 
         /// <summary>
-        /// Returns one sample from the distribution.
+        /// Samples a sequence of this random variable.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
         /// <param name="p">The p parameter</param>
-        /// <returns>
-        /// One sample from the distribution implied by <paramref name="p"/>.
-        /// </returns>
-        private static int DoSample(Random rnd, double p)
+        public static IEnumerable<int> Samples(Random rnd, double p)
         {
-            return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p));
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, p);
+            }
         }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Hypergeometric.cs b/src/Numerics/Distributions/Discrete/Hypergeometric.cs
index 68ad0679..4901135b 100644
--- a/src/Numerics/Distributions/Discrete/Hypergeometric.cs
+++ b/src/Numerics/Distributions/Discrete/Hypergeometric.cs
@@ -52,22 +52,22 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The size of the population.
         /// </summary>
-        private int _populationSize;
+        int _populationSize;
 
         /// <summary>
         /// The m parameter of the distribution.
         /// </summary>
-        private int _m;
+        int _m;
 
         /// <summary>
         /// The n parameter (number to draw) of the distribution.
         /// </summary>
-        private int _n;
+        int _n;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the Hypergeometric class.
@@ -87,7 +87,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="total">The Total parameter of the distribution.</param>
         /// <param name="m">The m parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
-        private void SetParameters(int total, int m, int n)
+        void SetParameters(int total, int m, int n)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(total, m, n))
             {
@@ -106,7 +106,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="m">The m parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(int total, int m, int n)
+        static bool IsValidParameterSet(int total, int m, int n)
         {
             if (total < 0 || m < 0 || n < 0)
             {
@@ -166,10 +166,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -187,10 +184,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return (double)_m * _n / _populationSize;
-            }
+            get { return (double)_m * _n / _populationSize; }
         }
 
         /// <summary>
@@ -198,10 +192,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0));
-            }
+            get { return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0)); }
         }
 
         /// <summary>
@@ -225,10 +216,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0));
-            }
+            get { return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0)); }
         }
 
         /// <summary>
@@ -320,27 +308,6 @@ namespace MathNet.Numerics.Distributions
             return Math.Log(Probability(k));
         }
 
-        /// <summary>
-        /// Samples a Hypergeometric distributed random variable.
-        /// </summary>
-        /// <returns>The number of successes in n trials.</returns>
-        public int Sample()
-        {
-            return DoSample(RandomSource, _populationSize, _m, _n);
-        }
-
-        /// <summary>
-        /// Samples an array of Hypergeometric distributed random variables.
-        /// </summary>
-        /// <returns>a sequence of successes in n trials.</returns>
-        public IEnumerable<int> Samples()
-        {
-            while (true)
-            {
-                yield return DoSample(RandomSource, _populationSize, _m, _n);
-            }
-        }
-
         #endregion
 
         /// <summary>
@@ -351,7 +318,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="m">The m parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
         /// <returns>a random number from the Hypergeometric distribution.</returns>
-        private static int DoSample(Random rnd, int size, int m, int n)
+        internal static int SampleUnchecked(Random rnd, int size, int m, int n)
         {
             var x = 0;
 
@@ -367,10 +334,68 @@ namespace MathNet.Numerics.Distributions
 
                 size--;
                 n--;
-            } 
+            }
             while (0 < n);
 
             return x;
         }
+
+        /// <summary>
+        /// Samples a Hypergeometric distributed random variable.
+        /// </summary>
+        /// <returns>The number of successes in n trials.</returns>
+        public int Sample()
+        {
+            return SampleUnchecked(RandomSource, _populationSize, _m, _n);
+        }
+
+        /// <summary>
+        /// Samples an array of Hypergeometric distributed random variables.
+        /// </summary>
+        /// <returns>a sequence of successes in n trials.</returns>
+        public IEnumerable<int> Samples()
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(RandomSource, _populationSize, _m, _n);
+            }
+        }
+
+        /// <summary>
+        /// Samples a random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="populationSize">The population size.</param>
+        /// <param name="m">The m parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        public static int Sample(Random rnd, int populationSize, int m, int n)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            return SampleUnchecked(rnd, populationSize, m, n);
+        }
+
+        /// <summary>
+        /// Samples a sequence of this random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="populationSize">The population size.</param>
+        /// <param name="m">The m parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        public static IEnumerable<int> Samples(Random rnd, int populationSize, int m, int n)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, populationSize, m, n);
+            }
+        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/NegativeBinomial.cs b/src/Numerics/Distributions/Discrete/NegativeBinomial.cs
index e9bfdf29..c619a974 100644
--- a/src/Numerics/Distributions/Discrete/NegativeBinomial.cs
+++ b/src/Numerics/Distributions/Discrete/NegativeBinomial.cs
@@ -46,32 +46,26 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The r parameter of the distribution.
         /// </summary>
-        private double _r;
+        double _r;
 
         /// <summary>
         /// The p parameter of the distribution.
         /// </summary>
-        private double _p;
+        double _p;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Gets or sets the number of trials.
         /// </summary>
         public double R
         {
-            get
-            {
-                return _r;
-            }
+            get { return _r; }
 
-            set
-            {
-                SetParameters(value, _p);
-            }
+            set { SetParameters(value, _p); }
         }
 
         /// <summary>
@@ -79,15 +73,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double P
         {
-            get
-            {
-                return _p;
-            }
+            get { return _p; }
 
-            set
-            {
-                SetParameters(_r, value);
-            }
+            set { SetParameters(_r, value); }
         }
 
         /// <summary>
@@ -107,7 +95,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="r">The number of trials.</param>
         /// <param name="p">The probability of a trial resulting in success.</param>
         /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double r, double p)
+        void SetParameters(double r, double p)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
             {
@@ -124,7 +112,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="r">The number of trials.</param>
         /// <param name="p">The probability of a trial resulting in success.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>        
-        private static bool IsValidParameterSet(double r, double p)
+        static bool IsValidParameterSet(double r, double p)
         {
             if (r < 0.0 || Double.IsNaN(r))
             {
@@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return _r * (1.0 - _p) / _p;
-            }
+            get { return _r * (1.0 - _p) / _p; }
         }
 
         /// <summary>
@@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _r * (1.0 - _p) / (_p * _p);
-            }
+            get { return _r * (1.0 - _p) / (_p * _p); }
         }
 
         /// <summary>
@@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(_r * (1.0 - _p)) / _p;
-            }
+            get { return Math.Sqrt(_r * (1.0 - _p)) / _p; }
         }
 
         /// <summary>
@@ -211,10 +187,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Entropy
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -222,10 +195,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p));
-            }
+            get { return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p)); }
         }
 
         /// <summary>
@@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0;
-            }
+            get { return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0; }
         }
 
         /// <summary>
@@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -269,10 +233,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -280,10 +241,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return int.MaxValue;
-            }
+            get { return int.MaxValue; }
         }
 
         /// <summary>
@@ -296,10 +254,10 @@ namespace MathNet.Numerics.Distributions
         public double Probability(int k)
         {
             var ln = SpecialFunctions.GammaLn(_r + k)
-                     - SpecialFunctions.GammaLn(_r)
-                     - SpecialFunctions.GammaLn(k + 1.0)
-                     + (_r * Math.Log(_p))
-                     + (k * Math.Log(1.0 - _p));
+                - SpecialFunctions.GammaLn(_r)
+                - SpecialFunctions.GammaLn(k + 1.0)
+                + (_r * Math.Log(_p))
+                + (k * Math.Log(1.0 - _p));
             return Math.Exp(ln);
         }
 
@@ -313,20 +271,44 @@ namespace MathNet.Numerics.Distributions
         public double ProbabilityLn(int k)
         {
             var ln = SpecialFunctions.GammaLn(_r + k)
-                     - SpecialFunctions.GammaLn(_r)
-                     - SpecialFunctions.GammaLn(k + 1.0)
-                     + (_r * Math.Log(_p))
-                     + (k * Math.Log(1.0 - _p));
+                - SpecialFunctions.GammaLn(_r)
+                - SpecialFunctions.GammaLn(k + 1.0)
+                + (_r * Math.Log(_p))
+                + (k * Math.Log(1.0 - _p));
             return ln;
         }
 
+        #endregion
+
+        /// <summary>
+        /// Samples a negative binomial distributed random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="r">The r parameter.</param>
+        /// <param name="p">The p parameter.</param>
+        /// <returns>a sample from the distribution.</returns>
+        internal static int SampleUnchecked(Random rnd, double r, double p)
+        {
+            var lambda = Gamma.SampleUnchecked(rnd, r, p);
+            var c = Math.Exp(-lambda);
+            var p1 = 1.0;
+            var k = 0;
+            do
+            {
+                k = k + 1;
+                p1 = p1 * rnd.NextDouble();
+            }
+            while (p1 >= c);
+            return k - 1;
+        }
+
         /// <summary>
         /// Samples a <c>NegativeBinomial</c> distributed random variable.
         /// </summary>
         /// <returns>a sample from the distribution.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _r, _p);
+            return SampleUnchecked(RandomSource, _r, _p);
         }
 
         /// <summary>
@@ -337,43 +319,43 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _r, _p);
+                yield return SampleUnchecked(RandomSource, _r, _p);
             }
         }
 
-        #endregion
-
         /// <summary>
-        /// Samples a negative binomial distributed random variable.
+        /// Samples a random variable.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
         /// <param name="r">The r parameter.</param>
         /// <param name="p">The p parameter.</param>
-        /// <returns>a sample from the distribution.</returns>
-        private static int DoSample(Random rnd, double r, double p)
+        public static int Sample(Random rnd, double r, double p)
         {
-            var lambda = Gamma.Sample(rnd, r, p);
-            return DoPoissonSample(rnd, lambda);
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            return SampleUnchecked(rnd, r, p);
         }
 
         /// <summary>
-        /// Use Knuth's method to generate Poisson distributed random variables.
+        /// Samples a sequence of this random variable.
         /// </summary>
-        /// <param name="rnd">The random number generator.</param>
-        /// <param name="lambda">The lambda value.</param>
-        /// <returns>a Poisson distributed random number.</returns>
-        private static int DoPoissonSample(Random rnd, double lambda)
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="r">The r parameter.</param>
+        /// <param name="p">The p parameter.</param>
+        public static IEnumerable<int> Samples(Random rnd, double r, double p)
         {
-            var c = Math.Exp(-lambda);
-            var p = 1.0;
-            var k = 0;
-            do
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
             {
-                k = k + 1;
-                p = p * rnd.NextDouble();
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, r, p);
             }
-            while (p >= c);
-            return k - 1;
         }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Poisson.cs b/src/Numerics/Distributions/Discrete/Poisson.cs
index f3fc4b22..39b3d432 100644
--- a/src/Numerics/Distributions/Discrete/Poisson.cs
+++ b/src/Numerics/Distributions/Discrete/Poisson.cs
@@ -43,27 +43,21 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The Poisson distribution parameter λ.
         /// </summary>
-        private double _lambda;
+        double _lambda;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Gets or sets the Poisson distribution parameter λ.
         /// </summary>
         public double Lambda
         {
-            get
-            {
-                return _lambda;
-            }
+            get { return _lambda; }
 
-            set
-            {
-                SetParameters(value);
-            }
+            set { SetParameters(value); }
         }
 
         /// <summary>
@@ -82,7 +76,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="lambda">The mean (λ) of the distribution.</param>
         /// <exception cref="System.ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
-        private void SetParameters(double lambda)
+        void SetParameters(double lambda)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
             {
@@ -97,7 +91,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="lambda">The mean (λ) of the distribution.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double lambda)
+        static bool IsValidParameterSet(double lambda)
         {
             return lambda > 0.00;
         }
@@ -120,10 +114,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -141,10 +132,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return _lambda;
-            }
+            get { return _lambda; }
         }
 
         /// <summary>
@@ -152,10 +140,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Variance
         {
-            get
-            {
-                return _lambda;
-            }
+            get { return _lambda; }
         }
 
         /// <summary>
@@ -163,10 +148,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(_lambda);
-            }
+            get { return Math.Sqrt(_lambda); }
         }
 
         /// <summary>
@@ -175,10 +157,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
         public double Entropy
         {
-            get
-            {
-                return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda));
-            }
+            get { return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda)); }
         }
 
         /// <summary>
@@ -186,10 +165,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get
-            {
-                return 1.0 / Math.Sqrt(_lambda);
-            }
+            get { return 1.0 / Math.Sqrt(_lambda); }
         }
 
         /// <summary>
@@ -197,10 +173,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 0;
-            }
+            get { return 0; }
         }
 
         /// <summary>
@@ -208,10 +181,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return int.MaxValue;
-            }
+            get { return int.MaxValue; }
         }
 
         /// <summary>
@@ -233,10 +203,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return (int)Math.Floor(_lambda);
-            }
+            get { return (int)Math.Floor(_lambda); }
         }
 
         /// <summary>
@@ -245,10 +212,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
         public int Median
         {
-            get
-            {
-                return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda));
-            }
+            get { return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda)); }
         }
 
         /// <summary>
@@ -271,71 +235,15 @@ namespace MathNet.Numerics.Distributions
             return -_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k);
         }
 
-        /// <summary>
-        /// Samples a Poisson distributed random variable.
-        /// </summary>
-        /// <returns>A sample from the Poisson distribution.</returns>
-        public int Sample()
-        {
-            return DoSample(RandomSource, _lambda);
-        }
-
-        /// <summary>
-        /// Samples an array of Poisson distributed random variables.
-        /// </summary>
-        /// <returns>a sequence of successes in N trials.</returns>
-        public IEnumerable<int> Samples()
-        {
-            while (true)
-            {
-                yield return DoSample(RandomSource, _lambda);
-            }
-        }
-
         #endregion
 
-        /// <summary>
-        /// Samples a Poisson distributed random variable.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="lambda">The Poisson distribution parameter λ.</param>
-        /// <returns>A sample from the Poisson distribution.</returns>
-        public static int Sample(Random rnd, double lambda)
-        {
-            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
-            {
-                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
-            }
-
-            return DoSample(rnd, lambda);
-        }
-
-        /// <summary>
-        /// Samples a sequence of Poisson distributed random variables.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="lambda">The Poisson distribution parameter λ.</param>
-        /// <returns>a sequence of samples from the distribution.</returns>
-        public static IEnumerable<int> Samples(Random rnd, double lambda)
-        {
-            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
-            {
-                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
-            }
-
-            while (true)
-            {
-                yield return DoSample(rnd, lambda);
-            }
-        }
-
         /// <summary>
         /// Generates one sample from the Poisson distribution.
         /// </summary>
         /// <param name="rnd">The random source to use.</param>
         /// <param name="lambda">The Poisson distribution parameter λ.</param>
         /// <returns>A random sample from the Poisson distribution.</returns>
-        private static int DoSample(Random rnd, double lambda)
+        internal static int SampleUnchecked(Random rnd, double lambda)
         {
             return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda);
         }
@@ -346,7 +254,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="rnd">The random source to use.</param>
         /// <param name="lambda">The Poisson distribution parameter λ.</param>
         /// <returns>A random sample from the Poisson distribution.</returns>
-        private static int DoSampleShort(Random rnd, double lambda)
+        static int DoSampleShort(Random rnd, double lambda)
         {
             var limit = Math.Exp(-lambda);
             var count = 0;
@@ -367,7 +275,7 @@ namespace MathNet.Numerics.Distributions
         /// <remarks>"Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson,
         /// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
         /// The article is on pages 29-35. The algorithm given here is on page 32. </remarks>
-        private static int DoSampleLarge(Random rnd, double lambda)
+        static int DoSampleLarge(Random rnd, double lambda)
         {
             var c = 0.767 - (3.36 / lambda);
             var beta = Math.PI / Math.Sqrt(3.0 * lambda);
@@ -395,5 +303,61 @@ namespace MathNet.Numerics.Distributions
                 }
             }
         }
+
+        /// <summary>
+        /// Samples a Poisson distributed random variable.
+        /// </summary>
+        /// <returns>A sample from the Poisson distribution.</returns>
+        public int Sample()
+        {
+            return SampleUnchecked(RandomSource, _lambda);
+        }
+
+        /// <summary>
+        /// Samples an array of Poisson distributed random variables.
+        /// </summary>
+        /// <returns>a sequence of successes in N trials.</returns>
+        public IEnumerable<int> Samples()
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(RandomSource, _lambda);
+            }
+        }
+
+        /// <summary>
+        /// Samples a Poisson distributed random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lambda">The Poisson distribution parameter λ.</param>
+        /// <returns>A sample from the Poisson distribution.</returns>
+        public static int Sample(Random rnd, double lambda)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            return SampleUnchecked(rnd, lambda);
+        }
+
+        /// <summary>
+        /// Samples a sequence of Poisson distributed random variables.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lambda">The Poisson distribution parameter λ.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static IEnumerable<int> Samples(Random rnd, double lambda)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, lambda);
+            }
+        }
     }
 }
diff --git a/src/Numerics/Distributions/Discrete/Zipf.cs b/src/Numerics/Distributions/Discrete/Zipf.cs
index 34849ea7..8bd0bc88 100644
--- a/src/Numerics/Distributions/Discrete/Zipf.cs
+++ b/src/Numerics/Distributions/Discrete/Zipf.cs
@@ -47,17 +47,17 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// The s parameter of the distribution.
         /// </summary>
-        private double _s;
+        double _s;
 
         /// <summary>
         /// The n parameter of the distribution.
         /// </summary>
-        private int _n;
+        int _n;
 
         /// <summary>
         /// The distribution's random number generator.
         /// </summary>
-        private Random _random;
+        Random _random;
 
         /// <summary>
         /// Initializes a new instance of the <see cref="Zipf"/> class. 
@@ -79,7 +79,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="s">The s parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
-        private void SetParameters(double s, int n)
+        void SetParameters(double s, int n)
         {
             if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
             {
@@ -96,7 +96,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="s">The s parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
         /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
-        private static bool IsValidParameterSet(double s, int n)
+        static bool IsValidParameterSet(double s, int n)
         {
             if (n <= 0 || s <= 0 || Double.IsNaN(s))
             {
@@ -111,15 +111,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double S
         {
-            get
-            {
-                return _s;
-            }
+            get { return _s; }
 
-            set
-            {
-                SetParameters(value, _n);
-            }
+            set { SetParameters(value, _n); }
         }
 
         /// <summary>
@@ -127,15 +121,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int N
         {
-            get
-            {
-                return _n;
-            }
+            get { return _n; }
 
-            set
-            {
-                SetParameters(_s, value);
-            }
+            set { SetParameters(_s, value); }
         }
 
         /// <summary>
@@ -154,10 +142,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public Random RandomSource
         {
-            get
-            {
-                return _random;
-            }
+            get { return _random; }
 
             set
             {
@@ -175,10 +160,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Mean
         {
-            get
-            {
-                return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s);
-            }
+            get { return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s); }
         }
 
         /// <summary>
@@ -203,10 +185,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get
-            {
-                return Math.Sqrt(Variance);
-            }
+            get { return Math.Sqrt(Variance); }
         }
 
         /// <summary>
@@ -266,10 +245,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Mode
         {
-            get
-            {
-                return 1;
-            }
+            get { return 1; }
         }
 
         /// <summary>
@@ -277,10 +253,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Median
         {
-            get
-            {
-                throw new NotSupportedException();
-            }
+            get { throw new NotSupportedException(); }
         }
 
         /// <summary>
@@ -288,10 +261,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Minimum
         {
-            get
-            {
-                return 1;
-            }
+            get { return 1; }
         }
 
         /// <summary>
@@ -299,12 +269,9 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Maximum
         {
-            get
-            {
-                return _n;
-            }
+            get { return _n; }
         }
-        
+
         /// <summary>
         /// Computes values of the probability mass function.
         /// </summary>
@@ -329,13 +296,45 @@ namespace MathNet.Numerics.Distributions
             return Math.Log(Probability(k));
         }
 
+        #endregion
+
+        /// <summary>
+        /// Generates a sample from the Zipf distribution without doing parameter checking.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="s">The s parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        /// <returns>a random number from the Zipf distribution.</returns>
+        internal static int SampleUnchecked(Random rnd, double s, int n)
+        {
+            var r = 0.0;
+            while (r == 0.0)
+            {
+                r = rnd.NextDouble();
+            }
+
+            var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s);
+            int i;
+            var sum = 0.0;
+            for (i = 1; i <= n; i++)
+            {
+                sum += p / Math.Pow(i, s);
+                if (sum >= r)
+                {
+                    break;
+                }
+            }
+
+            return i;
+        }
+
         /// <summary>
         /// Draws a random sample from the distribution.
         /// </summary>
         /// <returns>a sample from the distribution.</returns>
         public int Sample()
         {
-            return DoSample(RandomSource, _s, _n);
+            return SampleUnchecked(RandomSource, _s, _n);
         }
 
         /// <summary>
@@ -346,40 +345,43 @@ namespace MathNet.Numerics.Distributions
         {
             while (true)
             {
-                yield return DoSample(RandomSource, _s, _n);
+                yield return SampleUnchecked(RandomSource, _s, _n);
             }
         }
 
-        #endregion
-
         /// <summary>
-        /// Generates a sample from the Zipf distribution without doing parameter checking.
+        /// Samples a random variable.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
         /// <param name="s">The s parameter of the distribution.</param>
         /// <param name="n">The n parameter of the distribution.</param>
-        /// <returns>a random number from the Zipf distribution.</returns>
-        private static int DoSample(Random rnd, double s, int n)
+        public static int Sample(Random rnd, double s, int n)
         {
-            var r = 0.0;
-            while (r == 0.0)
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
             {
-                r = rnd.NextDouble();
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
             }
 
-            var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s);
-            int i;
-            var sum = 0.0;
-            for (i = 1; i <= n; i++)
+            return SampleUnchecked(rnd, s, n);
+        }
+
+        /// <summary>
+        /// Samples a sequence of this random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="s">The s parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        public static IEnumerable<int> Samples(Random rnd, double s, int n)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
             {
-                sum += p / Math.Pow(i, s);
-                if (sum >= r)
-                {
-                    break;
-                }
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
             }
 
-            return i;
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, s, n);
+            }
         }
     }
 }
diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs
index ba424745..6816902d 100644
--- a/src/Numerics/Distributions/Multivariate/Multinomial.cs
+++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs
@@ -368,7 +368,7 @@ namespace MathNet.Numerics.Distributions
 
             for (var i = 0; i < n; i++)
             {
-                ret[Categorical.DoSample(rnd, cp)]++;
+                ret[Categorical.SampleUnchecked(rnd, cp)]++;
             }
 
             return ret;
@@ -399,7 +399,7 @@ namespace MathNet.Numerics.Distributions
 
                 for (var i = 0; i < n; i++)
                 {
-                    ret[Categorical.DoSample(rnd, cp)]++;
+                    ret[Categorical.SampleUnchecked(rnd, cp)]++;
                 }
 
                 yield return ret;