Browse Source

Distributions: array-sampling 2

pull/222/head
Christoph Ruegg 12 years ago
parent
commit
4b7d61e920
  1. 72
      src/Numerics/Distributions/Exponential.cs
  2. 1
      src/Numerics/Distributions/LogNormal.cs
  3. 75
      src/Numerics/Distributions/Pareto.cs
  4. 63
      src/Numerics/Distributions/Rayleigh.cs
  5. 87
      src/Numerics/Distributions/Triangular.cs
  6. 95
      src/Numerics/Distributions/Weibull.cs

72
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);
}
}
}

1
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;

75
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);
}
}
}

63
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);
}
}
}

87
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);
}
}
}

95
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);
}
}
}

Loading…
Cancel
Save