Distributions: array-sampling (wip)

src/FSharp/Distributions.fs

@@ -44,117 +44,117 @@ module Sample =
     let mapSeq3 f dist1 dist2 dist3 : System.Random -> 'T seq = fun rng -> Seq.zip3 (dist1 rng) (dist2 rng) (dist3 rng) |> Seq.map (fun (d1, d2, d3) -> f d1 d2 d3)
 
     /// Bernoulli with probability (p).
-    let bernoulli p rng = Bernoulli.Sample(rng, p)
-    let bernoulliSeq p rng = Bernoulli.Samples(rng, p)
+    let bernoulli p (rng:System.Random) = Bernoulli.Sample(rng, p)
+    let bernoulliSeq p (rng:System.Random) = Bernoulli.Samples(rng, p)
 
     /// Beta with α and β shape parameters.
-    let beta a b rng = Beta.Sample(rng, a, b)
-    let betaSeq a b rng = Beta.Samples(rng, a, b)
+    let beta a b (rng:System.Random) = Beta.Sample(rng, a, b)
+    let betaSeq a b (rng:System.Random) = Beta.Samples(rng, a, b)
 
     /// Binomial with success probability (p) in each trial and number of trials (n).
-    let binomial p n rng = Binomial.Sample(rng, p, n)
-    let binomialSeq p n rng = Binomial.Samples(rng, p, n)
+    let binomial p n (rng:System.Random) = Binomial.Sample(rng, p, n)
+    let binomialSeq p n (rng:System.Random) = Binomial.Samples(rng, p, n)
 
     /// Categorical with an array of nonnegative ratios defining the relative probability mass (unnormalized).
-    let categorical probabilityMass rng = Categorical.Sample(rng, probabilityMass)
-    let categoricalSeq probabilityMass rng = Categorical.Samples(rng, probabilityMass)
+    let categorical probabilityMass (rng:System.Random) = Categorical.Sample(rng, probabilityMass)
+    let categoricalSeq probabilityMass (rng:System.Random) = Categorical.Samples(rng, probabilityMass)
 
     /// Cauchy with location (x0) and scale (γ).
-    let cauchy location scale rng = Cauchy.Sample(rng, location, scale)
-    let cauchySeq location scale rng = Cauchy.Samples(rng, location, scale)
+    let cauchy location scale (rng:System.Random) = Cauchy.Sample(rng, location, scale)
+    let cauchySeq location scale (rng:System.Random) = Cauchy.Samples(rng, location, scale)
 
     /// Chi with degrees of freedom (k).
-    let chi freedom rng = Chi.Sample(rng, freedom)
-    let chiSeq freedom rng = Chi.Samples(rng, freedom)
+    let chi freedom (rng:System.Random) = Chi.Sample(rng, freedom)
+    let chiSeq freedom (rng:System.Random) = Chi.Samples(rng, freedom)
 
     /// Chi-Squared with degrees of freedom (k).
-    let chiSquared freedom rng = ChiSquared.Sample(rng, freedom)
-    let chiSquaredSeq freedom rng = ChiSquared.Samples(rng, freedom)
+    let chiSquared freedom (rng:System.Random) = ChiSquared.Sample(rng, freedom)
+    let chiSquaredSeq freedom (rng:System.Random) = ChiSquared.Samples(rng, freedom)
 
     /// Continuous-Uniform with lower and upper bounds.
-    let continuousUniform lower upper rng = ContinuousUniform.Sample(rng, lower, upper)
-    let continuousUniformSeq lower upper rng = ContinuousUniform.Samples(rng, lower, upper)
+    let continuousUniform lower upper (rng:System.Random) = ContinuousUniform.Sample(rng, lower, upper)
+    let continuousUniformSeq lower upper (rng:System.Random) = ContinuousUniform.Samples(rng, lower, upper)
 
     /// Conway-Maxwell-Poisson with lambda (λ) and rate of decay (ν).
-    let conwayMaxwellPoisson lambda nu rng = ConwayMaxwellPoisson.Sample(rng, lambda, nu)
-    let conwayMaxwellPoissonSeq lambda nu rng = ConwayMaxwellPoisson.Samples(rng, lambda, nu)
+    let conwayMaxwellPoisson lambda nu (rng:System.Random) = ConwayMaxwellPoisson.Sample(rng, lambda, nu)
+    let conwayMaxwellPoissonSeq lambda nu (rng:System.Random) = ConwayMaxwellPoisson.Samples(rng, lambda, nu)
 
     /// Discrete-Uniform with lower and upper bounds (both inclusive).
-    let discreteUniform lower upper rng = DiscreteUniform.Sample(rng, lower, upper)
-    let discreteUniformSeq lower upper rng = DiscreteUniform.Samples(rng, lower, upper)
+    let discreteUniform lower upper (rng:System.Random) = DiscreteUniform.Sample(rng, lower, upper)
+    let discreteUniformSeq lower upper (rng:System.Random) = DiscreteUniform.Samples(rng, lower, upper)
 
     /// Erlang with shape (k) and rate or inverse scale (λ).
-    let erlang shape rate rng = Erlang.Sample(rng, shape, rate)
-    let erlangSeq shape rate rng = Erlang.Samples(rng, shape, rate)
+    let erlang shape rate (rng:System.Random) = Erlang.Sample(rng, shape, rate)
+    let erlangSeq shape rate (rng:System.Random) = Erlang.Samples(rng, shape, rate)
 
     /// Exponential with rate (λ).
-    let exponential rate rng = Exponential.Sample(rng, rate)
-    let exponentialSeq rate rng = Exponential.Samples(rng, rate)
+    let exponential rate (rng:System.Random) = Exponential.Sample(rng, rate)
+    let exponentialSeq rate (rng:System.Random) = Exponential.Samples(rng, rate)
 
     /// Fisher-Snedecor (F-Distribution) with first (d1) and second (d2) degree of freedom.
-    let fisherSnedecor d1 d2 rng = FisherSnedecor.Sample(rng, d1, d2)
-    let fisherSnedecorSeq d1 d2 rng = FisherSnedecor.Samples(rng, d1, d2)
+    let fisherSnedecor d1 d2 (rng:System.Random) = FisherSnedecor.Sample(rng, d1, d2)
+    let fisherSnedecorSeq d1 d2 (rng:System.Random) = FisherSnedecor.Samples(rng, d1, d2)
 
     /// Gamma with shape (k, α) and rate or inverse scale (β).
-    let gamma shape rate rng = Gamma.Sample(rng, shape, rate)
-    let gammaSeq shape rate rng = Gamma.Sample(rng, shape, rate)
+    let gamma shape rate (rng:System.Random) = Gamma.Sample(rng, shape, rate)
+    let gammaSeq shape rate (rng:System.Random) = Gamma.Sample(rng, shape, rate)
 
     /// Geometric with probability (p) of generating one.
-    let geometric p rng = Geometric.Sample(rng, p)
-    let geometricSeq p rng = Geometric.Samples(rng, p)
+    let geometric p (rng:System.Random) = Geometric.Sample(rng, p)
+    let geometricSeq p (rng:System.Random) = Geometric.Samples(rng, p)
 
     /// Hypergeometric with size of the population (N), number successes within the population (K, M) and number of draws without replacement (n).
-    let hypergeometric population success draws rng = Hypergeometric.Sample(rng, population, success, draws)
-    let hypergeometricSeq population success draws rng = Hypergeometric.Samples(rng, population, success, draws)
+    let hypergeometric population success draws (rng:System.Random) = Hypergeometric.Sample(rng, population, success, draws)
+    let hypergeometricSeq population success draws (rng:System.Random) = Hypergeometric.Samples(rng, population, success, draws)
 
     /// Inverse-Gamma with shape (α) and scale (β)
-    let inverseGamma shape scale rng = InverseGamma.Sample(rng, shape, scale)
-    let inverseGammaSeq shape scale rng = InverseGamma.Samples(rng, shape, scale)
+    let inverseGamma shape scale (rng:System.Random) = InverseGamma.Sample(rng, shape, scale)
+    let inverseGammaSeq shape scale (rng:System.Random) = InverseGamma.Samples(rng, shape, scale)
 
     /// Laplace with location (μ) and scale (b).
-    let laplace location scale rng = Laplace.Sample(rng, location, scale)
-    let laplaceSeq location scale rng = Laplace.Samples(rng, location, scale)
+    let laplace location scale (rng:System.Random) = Laplace.Sample(rng, location, scale)
+    let laplaceSeq location scale (rng:System.Random) = Laplace.Samples(rng, location, scale)
 
     /// Log-Normal with log-scale (μ) and shape (σ).
-    let logNormal mu sigma rng = LogNormal.Sample(rng, mu, sigma)
-    let logNormalSeq mu sigma rng = LogNormal.Samples(rng, mu, sigma)
+    let logNormal mu sigma (rng:System.Random) = LogNormal.Sample(rng, mu, sigma)
+    let logNormalSeq mu sigma (rng:System.Random) = LogNormal.Samples(rng, mu, sigma)
 
     /// Negative-Binomial with number of failures (r) until the experiment stopped and probability (p) of a trial resulting in success.
-    let negativeBinomial r p rng = NegativeBinomial.Sample(rng, r, p)
-    let negativeBinomialSeq r p rng = NegativeBinomial.Samples(rng, r, p)
+    let negativeBinomial r p (rng:System.Random) = NegativeBinomial.Sample(rng, r, p)
+    let negativeBinomialSeq r p (rng:System.Random) = NegativeBinomial.Samples(rng, r, p)
 
     /// Normal with mean (μ) and standard deviation (σ).
-    let normal mean stddev rng = Normal.Sample(rng, mean, stddev)
-    let normalSeq mean stddev rng = Normal.Samples(rng, mean, stddev)
+    let normal mean stddev (rng:System.Random) = Normal.Sample(rng, mean, stddev)
+    let normalSeq mean stddev (rng:System.Random) = Normal.Samples(rng, mean, stddev)
 
     /// Standard Gaussian.
-    let standard rng = Normal.Sample(rng, 0.0, 1.0)
-    let standardSeq rng = Normal.Samples(rng, 0.0, 1.0)
+    let standard (rng:System.Random) = Normal.Sample(rng, 0.0, 1.0)
+    let standardSeq (rng:System.Random) = Normal.Samples(rng, 0.0, 1.0)
 
     /// Pareto with scale (xm) and shape (α).
-    let pareto scale shape rng = Pareto.Sample(rng, scale, shape)
-    let paretoSeq scale shape rng = Pareto.Samples(rng, scale, shape)
+    let pareto scale shape (rng:System.Random) = Pareto.Sample(rng, scale, shape)
+    let paretoSeq scale shape (rng:System.Random) = Pareto.Samples(rng, scale, shape)
 
     /// Poisson with lambda (λ).
-    let poisson lambda rng = Poisson.Sample(rng, lambda)
-    let poissonSeq lambda rng = Poisson.Samples(rng, lambda)
+    let poisson lambda (rng:System.Random) = Poisson.Sample(rng, lambda)
+    let poissonSeq lambda (rng:System.Random) = Poisson.Samples(rng, lambda)
 
     /// Rayleigh with scale (σ).
-    let rayleigh scale rng = Rayleigh.Sample(rng, scale)
-    let rayleighSeq scale rng = Rayleigh.Sample(rng, scale)
+    let rayleigh scale (rng:System.Random) = Rayleigh.Sample(rng, scale)
+    let rayleighSeq scale (rng:System.Random) = Rayleigh.Sample(rng, scale)
 
     /// Stable with stability (α), skewness (β), scale (c) and location (μ).
-    let stable alpha beta scale location rng = Stable.Sample(rng, alpha, beta, scale, location)
-    let stableSeq alpha beta scale location rng = Stable.Samples(rng, alpha, beta, scale, location)
+    let stable alpha beta scale location (rng:System.Random) = Stable.Sample(rng, alpha, beta, scale, location)
+    let stableSeq alpha beta scale location (rng:System.Random) = Stable.Samples(rng, alpha, beta, scale, location)
 
     /// Student-T with location (μ), scale (σ) and degrees of freedom (ν).
-    let studentT location scale freedom rng = StudentT.Sample(rng, location, scale, freedom)
-    let studentTSeq location scale freedom rng = StudentT.Samples(rng, location, scale, freedom)
+    let studentT location scale freedom (rng:System.Random) = StudentT.Sample(rng, location, scale, freedom)
+    let studentTSeq location scale freedom (rng:System.Random) = StudentT.Samples(rng, location, scale, freedom)
 
     /// Weibull with shape (k) and scale (λ).
-    let weibull shape scale rng = Weibull.Sample(rng, shape, scale)
-    let weibullSeq shape scale rng = Weibull.Samples(rng, shape, scale)
+    let weibull shape scale (rng:System.Random) = Weibull.Sample(rng, shape, scale)
+    let weibullSeq shape scale (rng:System.Random) = Weibull.Samples(rng, shape, scale)
 
     /// Zipf with s and n parameters.
-    let zipf s n rng = Zipf.Sample(rng, s, n)
-    let zipfSeq s n rng = Zipf.Samples(rng, s, n)
+    let zipf s n (rng:System.Random) = Zipf.Sample(rng, s, n)
+    let zipfSeq s n (rng:System.Random) = Zipf.Samples(rng, s, n)

src/Numerics/Distributions/Cauchy.cs

@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -254,7 +255,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public double Sample()
         {
-            return _location + _scale*Math.Tan(Constants.Pi*(_random.NextDouble() - 0.5));
+            return SampleUnchecked(_random, _location, _scale);
         }
 
         /// <summary>
@@ -262,13 +263,35 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
+        {
+            return SamplesUnchecked(_random, _location, _scale);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double location, double scale)
+        {
+            return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
         {
             while (true)
             {
-                yield return _location + _scale*Math.Tan(Constants.Pi*(_random.NextDouble() - 0.5));
+                yield return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
             }
         }
 
+        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale)
+        {
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = location + scale*Math.Tan(Constants.Pi*(values[i] - 0.5));
+                }
+            });
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
@@ -342,7 +365,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
+            return SampleUnchecked(SystemRandomSource.Default, location, scale);
         }
 
         /// <summary>
@@ -356,10 +379,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
-            }
+            return SamplesUnchecked(rnd, location, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="location">The location (x0) of the distribution.</param>
+        /// <param name="scale">The scale (γ) of the distribution. Range: γ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double location, double scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, location, scale);
         }
 
         /// <summary>
@@ -370,7 +404,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double location, double scale)
         {
-            return Sample(SystemRandomSource.Default, location, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, location, scale);
         }
 
         /// <summary>
@@ -381,7 +417,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double location, double scale)
         {
-            return Samples(SystemRandomSource.Default, location, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, location, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="location">The location (x0) of the distribution.</param>
+        /// <param name="scale">The scale (γ) of the distribution. Range: γ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double location, double scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, location, scale);
         }
     }
 }

src/Numerics/Distributions/ContinuousUniform.cs

@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -258,7 +259,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return _lower + _random.NextDouble()*(_upper - _lower);
+            return SampleUnchecked(_random, _lower, _upper);
         }
 
         /// <summary>
@@ -267,12 +268,36 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
+            return SamplesUnchecked(_random, _lower, _upper);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double lower, double upper)
+        {
+            return lower + rnd.NextDouble()*(upper - lower);
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double lower, double upper)
+        {
+            double difference = upper - lower;
             while (true)
             {
-                yield return _lower + _random.NextDouble()*(_upper - _lower);
+                yield return lower + rnd.NextDouble()*difference;
             }
         }
 
+        static void SamplesUnchecked(System.Random rnd, double[] values, double lower, double upper)
+        {
+            rnd.NextDoubles(values);
+            var difference = upper - lower;
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = lower + values[i]*difference;
+                }
+            });
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
@@ -345,7 +370,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return lower + rnd.NextDouble()*(upper - lower);
+            return SampleUnchecked(rnd, lower, upper);
         }
 
         /// <summary>
@@ -359,10 +384,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return lower + rnd.NextDouble()*(upper - lower);
-            }
+            return SamplesUnchecked(rnd, lower, upper);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double lower, double upper)
+        {
+            if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, lower, upper);
         }
 
         /// <summary>
@@ -373,7 +410,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a uniformly distributed sample.</returns>
         public static double Sample(double lower, double upper)
         {
-            return Sample(SystemRandomSource.Default, lower, upper);
+            if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, lower, upper);
         }
 
         /// <summary>
@@ -384,7 +423,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of uniformly distributed samples.</returns>
         public static IEnumerable<double> Samples(double lower, double upper)
         {
-            return Samples(SystemRandomSource.Default, lower, upper);
+            if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, lower, upper);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double lower, double upper)
+        {
+            if (upper < lower) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, lower, upper);
         }
     }
 }

src/Numerics/Distributions/Laplace.cs

@@ -255,10 +255,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _location, _scale);
-            }
+            return SamplesUnchecked(_random, _location, _scale);
         }
 
         /// <summary>
@@ -271,7 +268,26 @@ namespace MathNet.Numerics.Distributions
         static double SampleUnchecked(System.Random rnd, double location, double scale)
         {
             var u = rnd.NextDouble() - 0.5;
-            return location - (scale * Math.Sign(u) * Math.Log(1.0 - (2.0 * Math.Abs(u))));
+            return location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
+        {
+            while (true)
+            {
+                var u = rnd.NextDouble() - 0.5;
+                yield return location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
+            }
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale)
+        {
+            rnd.NextDoubles(values);
+            for (int i = 0; i < values.Length; i++)
+            {
+                var u = values[i] - 0.5;
+                values[i] = location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
+            }
         }
 
         /// <summary>
@@ -344,10 +360,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, location, scale);
-            }
+            return SamplesUnchecked(rnd, location, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <param name="scale">The scale (b) of the distribution. Range: b > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double location, double scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, location, scale);
         }
 
         /// <summary>
@@ -358,7 +386,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double location, double scale)
         {
-            return Sample(SystemRandomSource.Default, location, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, location, scale);
         }
 
         /// <summary>
@@ -369,7 +399,24 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double location, double scale)
         {
-            return Samples(SystemRandomSource.Default, location, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, location, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <param name="scale">The scale (b) of the distribution. Range: b > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double location, double scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, location, scale);
         }
     }
 }

src/Numerics/Distributions/LogNormal.cs

@@ -31,9 +31,11 @@
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Reflection.Emit;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.Statistics;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -98,7 +100,7 @@ namespace MathNet.Numerics.Distributions
         public static LogNormal WithMeanVariance(double mean, double var, System.Random randomSource = null)
         {
             var sigma2 = Math.Log(var/(mean*mean) + 1.0);
-            return new LogNormal(Math.Log(mean) - sigma2 / 2.0, Math.Sqrt(sigma2), randomSource);
+            return new LogNormal(Math.Log(mean) - sigma2/2.0, Math.Sqrt(sigma2), randomSource);
         }
 
         /// <summary>
@@ -264,7 +266,7 @@ namespace MathNet.Numerics.Distributions
                 return 0.0;
             }
 
-            var a = (Math.Log(x) - _mu) / _sigma;
+            var a = (Math.Log(x) - _mu)/_sigma;
             return Math.Exp(-0.5*a*a)/(x*_sigma*Constants.Sqrt2Pi);
         }
 
@@ -281,7 +283,7 @@ namespace MathNet.Numerics.Distributions
                 return Double.NegativeInfinity;
             }
 
-            var a = (Math.Log(x) - _mu) / _sigma;
+            var a = (Math.Log(x) - _mu)/_sigma;
             return (-0.5*a*a) - Math.Log(x*_sigma) - Constants.LogSqrt2Pi;
         }
 
@@ -316,7 +318,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return Math.Exp(Normal.Sample(_random, _mu, _sigma));
+            return SampleUnchecked(_random, _mu, _sigma);
         }
 
         /// <summary>
@@ -325,7 +327,29 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            return Normal.Samples(_random, _mu, _sigma).Select(Math.Exp);
+            return SamplesUnchecked(_random, _mu, _sigma);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double mu, double sigma)
+        {
+            return Math.Exp(Normal.SampleUnchecked(rnd, mu, sigma));
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double mu, double sigma)
+        {
+            return Normal.SamplesUnchecked(rnd, mu, sigma).Select(Math.Exp);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double mu, double sigma)
+        {
+            Normal.SamplesUnchecked(rnd, values, mu, sigma);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = Math.Exp(values[i]);
+                }
+            });
         }
 
         /// <summary>
@@ -346,7 +370,7 @@ namespace MathNet.Numerics.Distributions
                 return 0.0;
             }
 
-            var a = (Math.Log(x) - mu) / sigma;
+            var a = (Math.Log(x) - mu)/sigma;
             return Math.Exp(-0.5*a*a)/(x*sigma*Constants.Sqrt2Pi);
         }
 
@@ -367,7 +391,7 @@ namespace MathNet.Numerics.Distributions
                 return Double.NegativeInfinity;
             }
 
-            var a = (Math.Log(x) - mu) / sigma;
+            var a = (Math.Log(x) - mu)/sigma;
             return (-0.5*a*a) - Math.Log(x*sigma) - Constants.LogSqrt2Pi;
         }
 
@@ -415,7 +439,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(System.Random rnd, double mu, double sigma)
         {
-            return Math.Exp(Normal.Sample(rnd, mu, sigma));
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(rnd, mu, sigma);
         }
 
         /// <summary>
@@ -427,7 +453,24 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(System.Random rnd, double mu, double sigma)
         {
-            return Normal.Samples(rnd, mu, sigma).Select(Math.Exp);
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(rnd, mu, sigma);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="mu">The log-scale (μ) of the distribution.</param>
+        /// <param name="sigma">The shape (σ) of the distribution. Range: σ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double mu, double sigma)
+        {
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, mu, sigma);
         }
 
         /// <summary>
@@ -438,7 +481,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double mu, double sigma)
         {
-            return Sample(SystemRandomSource.Default, mu, sigma);
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, mu, sigma);
         }
 
         /// <summary>
@@ -449,7 +494,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double mu, double sigma)
         {
-            return Samples(SystemRandomSource.Default, mu, sigma);
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, mu, sigma);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="mu">The log-scale (μ) of the distribution.</param>
+        /// <param name="sigma">The shape (σ) of the distribution. Range: σ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double mu, double sigma)
+        {
+            if (sigma < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, mu, sigma);
         }
     }
 }

src/Numerics/Distributions/Normal.cs

@@ -324,7 +324,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return _mean + (_stdDev*SampleStandardBoxMuller(_random).Item1);
+            return SampleUnchecked(_random, _mean, _stdDev);
         }
 
         /// <summary>
@@ -333,33 +333,93 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
+            return SamplesUnchecked(_random, _mean, _stdDev);
+        }
+
+        internal static double SampleUnchecked(System.Random rnd, double mean, double stddev)
+        {
+            double x, y;
+            while (!PolarTransform(rnd.NextDouble(), rnd.NextDouble(), out x, out y))
+            {
+            }
+            return mean + (stddev*x);
+        }
+
+        internal static IEnumerable<double> SamplesUnchecked(System.Random rnd, double mean, double stddev)
+        {
+            double x, y;
             while (true)
             {
-                var sample = SampleStandardBoxMuller(_random);
-                yield return _mean + (_stdDev*sample.Item1);
-                yield return _mean + (_stdDev*sample.Item2);
+                if (!PolarTransform(rnd.NextDouble(), rnd.NextDouble(), out x, out y))
+                {
+                    continue;
+                }
+                yield return mean + (stddev*x);
+                yield return mean + (stddev*y);
             }
         }
 
-        /// <summary>
-        /// Samples a pair of standard normal distributed random variables using the <i>Box-Muller</i> algorithm.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <returns>a pair of random numbers from the standard normal distribution.</returns>
-        static Tuple<double, double> SampleStandardBoxMuller(System.Random rnd)
+        internal static void SamplesUnchecked(System.Random rnd, double[] values, double mean, double stddev)
         {
-            var v1 = (2.0 * rnd.NextDouble()) - 1.0;
-            var v2 = (2.0 * rnd.NextDouble()) - 1.0;
-            var r = (v1 * v1) + (v2 * v2);
-            while (r >= 1.0 || r == 0.0)
+            if (values.Length == 0)
+            {
+                return;
+            }
+
+            // Since we only accept points within the unit circle
+            // we need to generate roughly 4/pi=1.27 times the numbers needed.
+            int n = (int)Math.Ceiling(values.Length*4*Constants.InvPi);
+            if (n.IsOdd())
+            {
+                n++;
+            }
+            var uniform = rnd.NextDoubles(n);
+
+            // Polar transform
+            double x, y;
+            int index = 0;
+            for (int i = 0; i < uniform.Length && index < values.Length; i += 2)
             {
-                v1 = (2.0 * rnd.NextDouble()) - 1.0;
-                v2 = (2.0 * rnd.NextDouble()) - 1.0;
-                r = (v1 * v1) + (v2 * v2);
+                if (!PolarTransform(uniform[i], uniform[i + 1], out x, out y))
+                {
+                    continue;
+                }
+                values[index++] = mean + stddev*x;
+                if (index == values.Length) return;
+                values[index++] = mean + stddev*y;
+                if (index == values.Length) return;
             }
 
-            var fac = Math.Sqrt(-2.0 * Math.Log(r) / r);
-            return new Tuple<double, double>(v1 * fac, v2 * fac);
+            // remaining, if any
+            while (index < values.Length)
+            {
+                if (!PolarTransform(rnd.NextDouble(), rnd.NextDouble(), out x, out y))
+                {
+                    continue;
+                }
+                values[index++] = mean + stddev*x;
+                if (index == values.Length) return;
+                values[index++] = mean + stddev*y;
+                if (index == values.Length) return;
+            }
+        }
+
+        static bool PolarTransform(double a, double b, out double x, out double y)
+        {
+            var v1 = (2.0*a) - 1.0;
+            var v2 = (2.0*b) - 1.0;
+            var r = (v1*v1) + (v2*v2);
+            if (r >= 1.0 || r == 0.0)
+            {
+                x = 0;
+                y = 0;
+                return false;
+            }
+
+            var fac = Math.Sqrt(-2.0*Math.Log(r)/r);
+            x = v1*fac;
+            y = v2*fac;
+            return true;
         }
 
         /// <summary>
@@ -439,7 +499,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return mean + (stddev*SampleStandardBoxMuller(rnd).Item1);
+            return SampleUnchecked(rnd, mean, stddev);
         }
 
         /// <summary>
@@ -453,12 +513,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                var sample = SampleStandardBoxMuller(rnd);
-                yield return mean + (stddev*sample.Item1);
-                yield return mean + (stddev*sample.Item2);
-            }
+            return SamplesUnchecked(rnd, mean, stddev);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="mean">The mean (μ) of the normal distribution.</param>
+        /// <param name="stddev">The standard deviation (σ) of the normal distribution. Range: σ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double mean, double stddev)
+        {
+            if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, mean, stddev);
         }
 
         /// <summary>
@@ -469,7 +539,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double mean, double stddev)
         {
-            return Sample(SystemRandomSource.Default, mean, stddev);
+            if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, mean, stddev);
         }
 
         /// <summary>
@@ -480,7 +552,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double mean, double stddev)
         {
-            return Samples(SystemRandomSource.Default, mean, stddev);
+            if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, mean, stddev);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="mean">The mean (μ) of the normal distribution.</param>
+        /// <param name="stddev">The standard deviation (σ) of the normal distribution. Range: σ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double mean, double stddev)
+        {
+            if (stddev < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, mean, stddev);
         }
     }
 }

src/Numerics/Distributions/Rayleigh.cs

@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -252,6 +253,18 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        static void SampleUnchecked(System.Random rnd, double[] values, double scale)
+        {
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = scale*Math.Sqrt(-2.0*Math.Log(values[i]));
+                }
+            });
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>

src/Numerics/Distributions/Triangular.cs

@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -295,6 +296,23 @@ namespace MathNet.Numerics.Distributions
             return Samples(_random, _lower, _upper, _mode);
         }
 
+        static void SampleUnchecked(System.Random rnd, double[] values, double lower, double upper, double mode)
+        {
+            double ml = mode - lower;
+            double ul = upper - lower;
+            double um = upper - mode;
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = values[i] < ml/ul
+                        ? lower + Math.Sqrt(values[i]*ul*ml)
+                        : upper - Math.Sqrt((1 - values[i])*ul*um);
+                }
+            });
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>

src/Numerics/Distributions/Weibull.cs

@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -289,6 +290,18 @@ namespace MathNet.Numerics.Distributions
             return scale*Math.Pow(-Math.Log(x), 1.0/shape);
         }
 
+        static void SampleUnchecked(System.Random rnd, double[] values, double shape, double scale)
+        {
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = scale*Math.Pow(-Math.Log(values[i]), 1.0/shape);
+                }
+            });
+        }
+
         /// <summary>
         /// Generates a sample from the Weibull distribution.
         /// </summary>

src/Numerics/Random/RandomExtensions.cs

@@ -46,7 +46,7 @@ namespace MathNet.Numerics.Random
         /// <param name="values">The array to fill with random values.</param>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static void NextDoubles(this System.Random rnd, double[] values)
         {
@@ -70,7 +70,7 @@ namespace MathNet.Numerics.Random
         /// <param name="count">The size of the array to fill.</param>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static double[] NextDoubles(this System.Random rnd, int count)
         {
@@ -84,7 +84,7 @@ namespace MathNet.Numerics.Random
         /// </summary>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static IEnumerable<double> NextDoubleSequence(this System.Random rnd)
         {
@@ -112,7 +112,7 @@ namespace MathNet.Numerics.Random
         /// <param name="count">The size of the array to fill.</param>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static byte[] NextBytes(this System.Random rnd, int count)
         {
@@ -132,7 +132,7 @@ namespace MathNet.Numerics.Random
         /// <seealso cref="NextFullRangeInt64"/>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static long NextInt64(this System.Random rnd)
         {
@@ -160,7 +160,7 @@ namespace MathNet.Numerics.Random
         /// <seealso cref="System.Random.Next()"/>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static int NextFullRangeInt32(this System.Random rnd)
         {
@@ -180,7 +180,7 @@ namespace MathNet.Numerics.Random
         /// <seealso cref="NextInt64"/>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static long NextFullRangeInt64(this System.Random rnd)
         {
@@ -199,7 +199,7 @@ namespace MathNet.Numerics.Random
         /// </returns>
         /// <remarks>
         /// This extension is thread-safe if and only if called on an random number
-        /// generator provided by Math.NET Nummerics or derived from the RandomSource class.
+        /// generator provided by Math.NET Numerics or derived from the RandomSource class.
         /// </remarks>
         public static decimal NextDecimal(this System.Random rnd)
         {