From 4b7d61e920751eda6baaa1fb83eb6846ca2d5d0f Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 9 Apr 2014 22:43:47 +0200
Subject: [PATCH] Distributions: array-sampling 2

---
 src/Numerics/Distributions/Exponential.cs | 72 ++++++++++++++---
 src/Numerics/Distributions/LogNormal.cs   |  1 -
 src/Numerics/Distributions/Pareto.cs      | 75 +++++++++++++++---
 src/Numerics/Distributions/Rayleigh.cs    | 63 +++++++++++----
 src/Numerics/Distributions/Triangular.cs  | 87 +++++++++++++++------
 src/Numerics/Distributions/Weibull.cs     | 95 ++++++++++++++---------
 6 files changed, 298 insertions(+), 95 deletions(-)

diff --git a/src/Numerics/Distributions/Exponential.cs b/src/Numerics/Distributions/Exponential.cs
index e6fc1194..c6937557 100644
--- a/src/Numerics/Distributions/Exponential.cs
+++ b/src/Numerics/Distributions/Exponential.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -241,10 +243,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _rate);
-            }
+            return SamplesUnchecked(_random, _rate);
         }
 
         /// <summary>
@@ -261,7 +260,30 @@ namespace MathNet.Numerics.Distributions
                 r = rnd.NextDouble();
             }
 
-            return -Math.Log(r) / rate;
+            return -Math.Log(r)/rate;
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double rate)
+        {
+            return rnd.NextDoubleSequence().Where(r => r != 0.0).Select(r => -Math.Log(r)/rate);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double rate)
+        {
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    // this happens very rarely
+                    var r = values[i];
+                    while (r == 0.0)
+                    {
+                        r = rnd.NextDouble();
+                    }
+                    values[i] = -Math.Log(r)/rate;
+                }
+            });
         }
 
         /// <summary>
@@ -344,10 +366,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, rate);
-            }
+            return SamplesUnchecked(rnd, rate);
+        }
+
+        /// <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="rate">The rate (λ) parameter 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 rate)
+        {
+            if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, rate);
         }
 
         /// <summary>
@@ -357,7 +390,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public static double Sample(double rate)
         {
-            return Sample(SystemRandomSource.Default, rate);
+            if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, rate);
         }
 
         /// <summary>
@@ -367,7 +402,22 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double rate)
         {
-            return Samples(SystemRandomSource.Default, rate);
+            if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, rate);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double rate)
+        {
+            if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, rate);
         }
     }
 }
diff --git a/src/Numerics/Distributions/LogNormal.cs b/src/Numerics/Distributions/LogNormal.cs
index 73b71cd7..cbf8f243 100644
--- a/src/Numerics/Distributions/LogNormal.cs
+++ b/src/Numerics/Distributions/LogNormal.cs
@@ -31,7 +31,6 @@
 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;
diff --git a/src/Numerics/Distributions/Pareto.cs b/src/Numerics/Distributions/Pareto.cs
index ba962f8b..3f4d1e98 100644
--- a/src/Numerics/Distributions/Pareto.cs
+++ b/src/Numerics/Distributions/Pareto.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -266,7 +268,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public double Sample()
         {
-            return _scale*Math.Pow(_random.NextDouble(), -1.0/_shape);
+            return SampleUnchecked(_random, _scale, _shape);
         }
 
         /// <summary>
@@ -275,11 +277,31 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            var power = -1.0/_shape;
-            while (true)
+            return SamplesUnchecked(_random, _scale, _shape);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double scale, double shape)
+        {
+            return scale*Math.Pow(rnd.NextDouble(), -1.0/shape);
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double scale, double shape)
+        {
+            var power = -1.0/shape;
+            return rnd.NextDoubleSequence().Select(x => scale*Math.Pow(x, power));
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double scale, double shape)
+        {
+            var power = -1.0/shape;
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
             {
-                yield return _scale*Math.Pow(_random.NextDouble(), power);
-            }
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = scale*Math.Pow(values[i], power);
+                }
+            });
         }
 
         /// <summary>
@@ -368,11 +390,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0 || shape <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            var power = -1.0 / shape;
-            while (true)
-            {
-                yield return scale*Math.Pow(rnd.NextDouble(), power);
-            }
+            return SamplesUnchecked(rnd, scale, shape);
+        }
+
+        /// <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="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">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 scale, double shape)
+        {
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, scale, shape);
         }
 
         /// <summary>
@@ -383,7 +416,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double scale, double shape)
         {
-            return Sample(SystemRandomSource.Default, scale, shape);
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, scale, shape);
         }
 
         /// <summary>
@@ -394,7 +429,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double scale, double shape)
         {
-            return Samples(SystemRandomSource.Default, scale, shape);
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, scale, shape);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double scale, double shape)
+        {
+            if (scale <= 0.0 || shape <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, scale, shape);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Rayleigh.cs b/src/Numerics/Distributions/Rayleigh.cs
index e254eb80..d6e0c205 100644
--- a/src/Numerics/Distributions/Rayleigh.cs
+++ b/src/Numerics/Distributions/Rayleigh.cs
@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.Threading;
@@ -238,7 +239,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public double Sample()
         {
-            return _scale*Math.Sqrt(-2.0*Math.Log(_random.NextDouble()));
+            return SampleUnchecked(_random, _scale);
         }
 
         /// <summary>
@@ -247,13 +248,20 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return _scale*Math.Sqrt(-2.0*Math.Log(_random.NextDouble()));
-            }
+            return SamplesUnchecked(_random, _scale);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double scale)
+        {
+            return scale*Math.Sqrt(-2.0*Math.Log(rnd.NextDouble()));
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double scale)
+        {
+            return rnd.NextDoubleSequence().Select(x => scale*Math.Sqrt(-2.0*Math.Log(x)));
         }
 
-        static void SampleUnchecked(System.Random rnd, double[] values, double scale)
+        static void SamplesUnchecked(System.Random rnd, double[] values, double scale)
         {
             rnd.NextDoubles(values);
             CommonParallel.For(0, values.Length, 4096, (a, b) =>
@@ -332,7 +340,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return scale*Math.Sqrt(-2.0*Math.Log(rnd.NextDouble()));
+            return SampleUnchecked(rnd, scale);
         }
 
         /// <summary>
@@ -345,10 +353,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return scale*Math.Sqrt(-2.0*Math.Log(rnd.NextDouble()));
-            }
+            return SamplesUnchecked(rnd, 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="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 scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, scale);
         }
 
         /// <summary>
@@ -358,7 +377,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double scale)
         {
-            return Sample(SystemRandomSource.Default, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, scale);
         }
 
         /// <summary>
@@ -368,7 +389,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double scale)
         {
-            return Samples(SystemRandomSource.Default, scale);
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, 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="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 scale)
+        {
+            if (scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, scale);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Triangular.cs b/src/Numerics/Distributions/Triangular.cs
index 991ea0d4..cb070b52 100644
--- a/src/Numerics/Distributions/Triangular.cs
+++ b/src/Numerics/Distributions/Triangular.cs
@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.Threading;
@@ -284,7 +285,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return Sample(_random, _lower, _upper, _mode);
+            return SampleUnchecked(_random, _lower, _upper, _mode);
         }
 
         /// <summary>
@@ -293,22 +294,38 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            return Samples(_random, _lower, _upper, _mode);
+            return SamplesUnchecked(_random, _lower, _upper, _mode);
         }
 
-        static void SampleUnchecked(System.Random rnd, double[] values, double lower, double upper, double mode)
+        static double SampleUnchecked(System.Random rnd, double lower, double upper, double mode)
         {
-            double ml = mode - lower;
-            double ul = upper - lower;
-            double um = upper - mode;
+            var u = rnd.NextDouble();
+            return u < (mode - lower)/(upper - lower)
+                ? lower + Math.Sqrt(u*(upper - lower)*(mode - lower))
+                : upper - Math.Sqrt((1 - u)*(upper - lower)*(upper - mode));
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double lower, double upper, double mode)
+        {
+            double ml = mode - lower, ul = upper - lower, um = upper - mode;
+            double u = ml/ul, v = ul*ml, w = ul*um;
+
+            return rnd.NextDoubleSequence().Select(x => x < u ? lower + Math.Sqrt(x*v) : upper - Math.Sqrt((1 - x)*w));
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double lower, double upper, double mode)
+        {
+            double ml = mode - lower, ul = upper - lower, um = upper - mode;
+            double u = ml/ul, v = ul*ml, w = ul*um;
+
             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);
+                    values[i] = values[i] < u
+                        ? lower + Math.Sqrt(values[i]*v)
+                        : upper - Math.Sqrt((1 - values[i])*w);
                 }
             });
         }
@@ -409,14 +426,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            var a = lower;
-            var b = upper;
-            var c = mode;
-            var u = rnd.NextDouble();
-
-            return u < (c - a)/(b - a)
-                ? a + Math.Sqrt(u*(b - a)*(c - a))
-                : b - Math.Sqrt((1 - u)*(b - a)*(b - c));
+            return SampleUnchecked(rnd, lower, upper, mode);
         }
 
         /// <summary>
@@ -431,10 +441,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return Sample(rnd, lower, upper, mode);
-            }
+            return SamplesUnchecked(rnd, lower, upper, mode);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
+        /// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
+        /// <param name="mode">Mode (most frequent value).  Range: lower ≤ mode ≤ upper</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double lower, double upper, double mode)
+        {
+            if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, lower, upper, mode);
         }
 
         /// <summary>
@@ -446,7 +468,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double lower, double upper, double mode)
         {
-            return Sample(SystemRandomSource.Default, lower, upper, mode);
+            if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, lower, upper, mode);
         }
 
         /// <summary>
@@ -458,7 +482,24 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double lower, double upper, double mode)
         {
-            return Samples(SystemRandomSource.Default, lower, upper, mode);
+            if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, lower, upper, mode);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
+        /// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
+        /// <param name="mode">Mode (most frequent value).  Range: lower ≤ mode ≤ upper</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double lower, double upper, double mode)
+        {
+            if (!(upper >= mode && mode >= lower)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, lower, upper, mode);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Weibull.cs b/src/Numerics/Distributions/Weibull.cs
index 9330789c..cbcfae15 100644
--- a/src/Numerics/Distributions/Weibull.cs
+++ b/src/Numerics/Distributions/Weibull.cs
@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.Threading;
@@ -276,32 +277,6 @@ namespace MathNet.Numerics.Distributions
             return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, _shape)*_scalePowShapeInv);
         }
 
-        /// <summary>
-        /// Generates one sample from the Weibull distribution. This method doesn't perform
-        /// any parameter checks.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="shape">The shape (k) of the Weibull distribution. Range: k > 0.</param>
-        /// <param name="scale">The scale (λ) of the Weibull distribution. Range: λ > 0.</param>
-        /// <returns>A sample from a Weibull distributed random variable.</returns>
-        static double SampleUnchecked(System.Random rnd, double shape, double scale)
-        {
-            var x = rnd.NextDouble();
-            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>
@@ -317,10 +292,32 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
+            return SamplesUnchecked(_random, _shape, _scale);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double shape, double scale)
+        {
+            var x = rnd.NextDouble();
+            return scale*Math.Pow(-Math.Log(x), 1.0/shape);
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double shape, double scale)
+        {
+            var exponent = 1.0/shape;
+            return rnd.NextDoubleSequence().Select(x => scale*Math.Pow(-Math.Log(x), exponent));
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double shape, double scale)
+        {
+            var exponent = 1.0/shape;
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
             {
-                yield return SampleUnchecked(_random, _shape, _scale);
-            }
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = scale*Math.Pow(-Math.Log(values[i]), exponent);
+                }
+            });
         }
 
         /// <summary>
@@ -421,10 +418,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, shape, scale);
-            }
+            return SamplesUnchecked(rnd, shape, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="shape">The shape (k) of the Weibull distribution. Range: k > 0.</param>
+        /// <param name="scale">The scale (λ) of the Weibull distribution. Range: λ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double shape, double scale)
+        {
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, shape, scale);
         }
 
         /// <summary>
@@ -435,7 +443,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double shape, double scale)
         {
-            return Sample(SystemRandomSource.Default, shape, scale);
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, shape, scale);
         }
 
         /// <summary>
@@ -446,7 +456,22 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double shape, double scale)
         {
-            return Samples(SystemRandomSource.Default, shape, scale);
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, shape, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="shape">The shape (k) of the Weibull distribution. Range: k > 0.</param>
+        /// <param name="scale">The scale (λ) of the Weibull distribution. Range: λ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double shape, double scale)
+        {
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, shape, scale);
         }
     }
 }