Browse Source

Distributions: consistent sample methods for discrete distributions #32

la-knuth
Christoph Ruegg 14 years ago
parent
commit
ef38c4508f
  1. 113
      src/Numerics/Distributions/Discrete/Bernoulli.cs
  2. 130
      src/Numerics/Distributions/Discrete/Binomial.cs
  3. 153
      src/Numerics/Distributions/Discrete/Categorical.cs
  4. 171
      src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs
  5. 118
      src/Numerics/Distributions/Discrete/DiscreteUniform.cs
  6. 119
      src/Numerics/Distributions/Discrete/Geometric.cs
  7. 115
      src/Numerics/Distributions/Discrete/Hypergeometric.cs
  8. 164
      src/Numerics/Distributions/Discrete/NegativeBinomial.cs
  9. 186
      src/Numerics/Distributions/Discrete/Poisson.cs
  10. 142
      src/Numerics/Distributions/Discrete/Zipf.cs
  11. 4
      src/Numerics/Distributions/Multivariate/Multinomial.cs

113
src/Numerics/Distributions/Discrete/Bernoulli.cs

@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The probability of generating a one. /// The probability of generating a one.
/// </summary> /// </summary>
private double _p; double _p;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the Bernoulli class. /// Initializes a new instance of the Bernoulli class.
@ -77,7 +77,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="p">The probability of generating a one.</param> /// <param name="p">The probability of generating a one.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double p) static bool IsValidParameterSet(double p)
{ {
if (p >= 0.0 && p <= 1.0) if (p >= 0.0 && p <= 1.0)
{ {
@ -92,7 +92,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="p">The probability of generating a one.</param> /// <param name="p">The probability of generating a one.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double p) void SetParameters(double p)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{ {
@ -107,15 +107,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double P public double P
{ {
get get { return _p; }
{
return _p;
}
set set { SetParameters(value); }
{
SetParameters(value);
}
} }
#region IDistribution Members #region IDistribution Members
@ -125,10 +119,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -146,10 +137,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return _p; }
{
return _p;
}
} }
/// <summary> /// <summary>
@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(_p * (1.0 - _p)); }
{
return Math.Sqrt(_p * (1.0 - _p));
}
} }
/// <summary> /// <summary>
@ -168,10 +153,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _p * (1.0 - _p); }
{
return _p * (1.0 - _p);
}
} }
/// <summary> /// <summary>
@ -179,10 +161,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p)); }
{
return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p));
}
} }
/// <summary> /// <summary>
@ -190,10 +169,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p)); }
{
return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p));
}
} }
/// <summary> /// <summary>
@ -201,10 +177,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -212,10 +185,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return 1; }
{
return 1;
}
} }
/// <summary> /// <summary>
@ -229,7 +199,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (x < 1.0) if (x < 1.0)
{ {
return 1.0 - _p; return 1.0 - _p;
@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return _p > 0.5 ? 1 : 0; }
{
return _p > 0.5 ? 1 : 0;
}
} }
/// <summary> /// <summary>
@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { throw new NotSupportedException("The median of the Bernoulli distribution is undefined."); }
{
throw new NotSupportedException("The median of the Bernoulli distribution is undefined.");
}
} }
/// <summary> /// <summary>
@ -299,13 +263,31 @@ namespace MathNet.Numerics.Distributions
return k == 1 ? Math.Log(_p) : Double.NegativeInfinity; return k == 1 ? Math.Log(_p) : Double.NegativeInfinity;
} }
#endregion
/// <summary>
/// Generates one sample from the Bernoulli distribution.
/// </summary>
/// <param name="rnd">The random source to use.</param>
/// <param name="p">The probability of generating a one.</param>
/// <returns>A random sample from the Bernoulli distribution.</returns>
internal static int SampleUnchecked(Random rnd, double p)
{
if (rnd.NextDouble() < p)
{
return 1;
}
return 0;
}
/// <summary> /// <summary>
/// Samples a Bernoulli distributed random variable. /// Samples a Bernoulli distributed random variable.
/// </summary> /// </summary>
/// <returns>A sample from the Bernoulli distribution.</returns> /// <returns>A sample from the Bernoulli distribution.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _p); return SampleUnchecked(RandomSource, _p);
} }
/// <summary> /// <summary>
@ -316,12 +298,10 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _p); yield return SampleUnchecked(RandomSource, _p);
} }
} }
#endregion
/// <summary> /// <summary>
/// Samples a Bernoulli distributed random variable. /// Samples a Bernoulli distributed random variable.
/// </summary> /// </summary>
@ -335,7 +315,7 @@ namespace MathNet.Numerics.Distributions
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
} }
return DoSample(rnd, p); return SampleUnchecked(rnd, p);
} }
/// <summary> /// <summary>
@ -353,24 +333,9 @@ namespace MathNet.Numerics.Distributions
while (true) while (true)
{ {
yield return DoSample(rnd, p); yield return SampleUnchecked(rnd, p);
} }
} }
/// <summary>
/// Generates one sample from the Bernoulli distribution.
/// </summary>
/// <param name="rnd">The random source to use.</param>
/// <param name="p">The probability of generating a one.</param>
/// <returns>A random sample from the Bernoulli distribution.</returns>
private static int DoSample(Random rnd, double p)
{
if (rnd.NextDouble() < p)
{
return 1;
}
return 0;
}
} }
} }

130
src/Numerics/Distributions/Discrete/Binomial.cs

@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// Stores the normalized binomial probability. /// Stores the normalized binomial probability.
/// </summary> /// </summary>
private double _p; double _p;
/// <summary> /// <summary>
/// The number of trials. /// The number of trials.
/// </summary> /// </summary>
private int _n; int _n;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the Binomial class. /// Initializes a new instance of the Binomial class.
@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="p">The success probability of a trial.</param> /// <param name="p">The success probability of a trial.</param>
/// <param name="n">The number of trials.</param> /// <param name="n">The number of trials.</param>
/// <returns><c>false</c> <paramref name="p"/> is not in the interval [0.0,1.0] or <paramref name="n"/> is negative, <c>true</c> otherwise.</returns> /// <returns><c>false</c> <paramref name="p"/> is not in the interval [0.0,1.0] or <paramref name="n"/> is negative, <c>true</c> otherwise.</returns>
private static bool IsValidParameterSet(double p, int n) static bool IsValidParameterSet(double p, int n)
{ {
if (p < 0.0 || p > 1.0 || Double.IsNaN(p)) if (p < 0.0 || p > 1.0 || Double.IsNaN(p))
{ {
@ -107,7 +107,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="n">The number of trials.</param> /// <param name="n">The number of trials.</param>
/// <exception cref="ArgumentOutOfRangeException">If <paramref name="p"/> is not in the interval [0.0,1.0].</exception> /// <exception cref="ArgumentOutOfRangeException">If <paramref name="p"/> is not in the interval [0.0,1.0].</exception>
/// <exception cref="ArgumentOutOfRangeException">If <paramref name="n"/> is negative.</exception> /// <exception cref="ArgumentOutOfRangeException">If <paramref name="n"/> is negative.</exception>
private void SetParameters(double p, int n) void SetParameters(double p, int n)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n)) if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{ {
@ -123,15 +123,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double P public double P
{ {
get get { return _p; }
{
return _p;
}
set set { SetParameters(value, _n); }
{
SetParameters(value, _n);
}
} }
/// <summary> /// <summary>
@ -139,15 +133,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int N public int N
{ {
get get { return _n; }
{
return _n;
}
set set { SetParameters(_p, value); }
{
SetParameters(_p, value);
}
} }
#region IDistribution Members #region IDistribution Members
@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return _p * _n; }
{
return _p * _n;
}
} }
/// <summary> /// <summary>
@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(_p * (1.0 - _p) * _n); }
{
return Math.Sqrt(_p * (1.0 - _p) * _n);
}
} }
/// <summary> /// <summary>
@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _p * (1.0 - _p) * _n; }
{
return _p * (1.0 - _p) * _n;
}
} }
/// <summary> /// <summary>
@ -234,10 +210,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p)); }
{
return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p));
}
} }
/// <summary> /// <summary>
@ -245,10 +218,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -256,10 +226,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return _n; }
{
return _n;
}
} }
/// <summary> /// <summary>
@ -273,7 +240,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (x > _n) if (x > _n)
{ {
return 1.0; return 1.0;
@ -303,7 +270,7 @@ namespace MathNet.Numerics.Distributions
{ {
return _n; return _n;
} }
if (_p == 0.0) if (_p == 0.0)
{ {
return 0; return 0;
@ -318,10 +285,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { return (int)Math.Floor(_p * _n); }
{
return (int)Math.Floor(_p * _n);
}
} }
/// <summary> /// <summary>
@ -345,7 +309,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 1.0; return 1.0;
} }
if (_p == 0.0) if (_p == 0.0)
{ {
return 0.0; return 0.0;
@ -355,7 +319,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 1.0; return 1.0;
} }
if (_p == 1.0) if (_p == 1.0)
{ {
return 0.0; return 0.0;
@ -385,7 +349,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (_p == 0.0) if (_p == 0.0)
{ {
return Double.NegativeInfinity; return Double.NegativeInfinity;
@ -395,7 +359,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (_p == 1.0) if (_p == 1.0)
{ {
return Double.NegativeInfinity; return Double.NegativeInfinity;
@ -404,13 +368,33 @@ namespace MathNet.Numerics.Distributions
return SpecialFunctions.BinomialLn(_n, k) + (k * Math.Log(_p)) + ((_n - k) * Math.Log(1.0 - _p)); return SpecialFunctions.BinomialLn(_n, k) + (k * Math.Log(_p)) + ((_n - k) * Math.Log(1.0 - _p));
} }
#endregion
/// <summary>
/// Generates a sample from the Binomial distribution without doing parameter checking.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The success probability of a trial; must be in the interval [0.0, 1.0].</param>
/// <param name="n">The number of trials; must be positive.</param>
/// <returns>The number of successful trials.</returns>
internal static int SampleUnchecked(Random rnd, double p, int n)
{
var k = 0;
for (var i = 0; i < n; i++)
{
k += rnd.NextDouble() < p ? 1 : 0;
}
return k;
}
/// <summary> /// <summary>
/// Samples a Binomially distributed random variable. /// Samples a Binomially distributed random variable.
/// </summary> /// </summary>
/// <returns>The number of successes in N trials.</returns> /// <returns>The number of successes in N trials.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _p, _n); return SampleUnchecked(RandomSource, _p, _n);
} }
/// <summary> /// <summary>
@ -421,12 +405,10 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _p, _n); yield return SampleUnchecked(RandomSource, _p, _n);
} }
} }
#endregion
/// <summary> /// <summary>
/// Samples a binomially distributed random variable. /// Samples a binomially distributed random variable.
/// </summary> /// </summary>
@ -441,7 +423,7 @@ namespace MathNet.Numerics.Distributions
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
} }
return DoSample(rnd, p, n); return SampleUnchecked(rnd, p, n);
} }
/// <summary> /// <summary>
@ -460,26 +442,8 @@ namespace MathNet.Numerics.Distributions
while (true) while (true)
{ {
yield return DoSample(rnd, p, n); yield return SampleUnchecked(rnd, p, n);
}
}
/// <summary>
/// Generates a sample from the Binomial distribution without doing parameter checking.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The success probability of a trial; must be in the interval [0.0, 1.0].</param>
/// <param name="n">The number of trials; must be positive.</param>
/// <returns>The number of successful trials.</returns>
private static int DoSample(Random rnd, double p, int n)
{
var k = 0;
for (var i = 0; i < n; i++)
{
k += rnd.NextDouble() < p ? 1 : 0;
} }
return k;
} }
} }
} }

153
src/Numerics/Distributions/Discrete/Categorical.cs

@ -50,12 +50,12 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// Stores the unnormalized categorical probabilities. /// Stores the unnormalized categorical probabilities.
/// </summary> /// </summary>
private double[] _p; double[] _p;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the Categorical class. /// Initializes a new instance of the Categorical class.
@ -110,7 +110,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="p">An array of nonnegative ratios: this array does not need to be normalized /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized
/// as this is often impossible using floating point arithmetic.</param> /// as this is often impossible using floating point arithmetic.</param>
/// <returns>If any of the probabilities are negative returns <c>false</c>, or if the sum of parameters is 0.0; otherwise <c>true</c></returns> /// <returns>If any of the probabilities are negative returns <c>false</c>, or if the sum of parameters is 0.0; otherwise <c>true</c></returns>
private static bool IsValidParameterSet(IEnumerable<double> p) static bool IsValidParameterSet(IEnumerable<double> p)
{ {
var sum = 0.0; var sum = 0.0;
foreach (double t in p) foreach (double t in p)
@ -119,7 +119,7 @@ namespace MathNet.Numerics.Distributions
{ {
return false; return false;
} }
sum += t; sum += t;
} }
@ -132,7 +132,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="p">An array of nonnegative ratios: this array does not need to be normalized /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized
/// as this is often impossible using floating point arithmetic.</param> /// as this is often impossible using floating point arithmetic.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double[] p) void SetParameters(double[] p)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{ {
@ -163,10 +163,7 @@ namespace MathNet.Numerics.Distributions
return p; return p;
} }
set set { SetParameters(value); }
{
SetParameters(value);
}
} }
#region IDistribution Members #region IDistribution Members
@ -176,10 +173,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -197,10 +191,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return _p.Mean(); }
{
return _p.Mean();
}
} }
/// <summary> /// <summary>
@ -208,10 +199,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return _p.StandardDeviation(); }
{
return _p.StandardDeviation();
}
} }
/// <summary> /// <summary>
@ -219,10 +207,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _p.Variance(); }
{
return _p.Variance();
}
} }
/// <summary> /// <summary>
@ -230,10 +215,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { return _p.Sum(p => p * Math.Log(p)); }
{
return _p.Sum(p => p * Math.Log(p));
}
} }
/// <summary> /// <summary>
@ -242,10 +224,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks> /// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks>
public double Skewness public double Skewness
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -253,10 +232,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -264,10 +240,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return _p.Length - 1; }
{
return _p.Length - 1;
}
} }
/// <summary> /// <summary>
@ -281,7 +254,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (x >= _p.Length) if (x >= _p.Length)
{ {
return 1.0; return 1.0;
@ -301,10 +274,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks> /// <remarks>Throws a <see cref="NotSupportedException"/>.</remarks>
public int Mode public int Mode
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -312,10 +282,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { return (int)_p.Median(); }
{
return (int)_p.Median();
}
} }
/// <summary> /// <summary>
@ -358,6 +325,47 @@ namespace MathNet.Numerics.Distributions
return Math.Log(_p[k]); return Math.Log(_p[k]);
} }
#endregion
/// <summary>
/// Computes the unnormalized cumulative distribution function. This method performs no
/// parameter checking.
/// </summary>
/// <param name="p">An array of nonnegative ratios: this array does not need to be normalized
/// as this is often impossible using floating point arithmetic.</param>
/// <returns>An array representing the unnormalized cumulative distribution function.</returns>
internal static double[] UnnormalizedCdf(double[] p)
{
var cp = (double[])p.Clone();
for (var i = 1; i < p.Length; i++)
{
cp[i] += cp[i - 1];
}
return cp;
}
/// <summary>
/// Returns one trials from the categorical distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="cdf">The cumulative distribution of the probability distribution.</param>
/// <returns>One sample from the categorical distribution implied by <paramref name="cdf"/>.</returns>
internal static int SampleUnchecked(Random rnd, double[] cdf)
{
// TODO : use binary search to speed up this procedure.
var u = rnd.NextDouble() * cdf[cdf.Length - 1];
var idx = 0;
while (u > cdf[idx])
{
idx++;
}
return idx;
}
/// <summary> /// <summary>
/// Samples a Binomially distributed random variable. /// Samples a Binomially distributed random variable.
/// </summary> /// </summary>
@ -376,8 +384,6 @@ namespace MathNet.Numerics.Distributions
return Samples(RandomSource, _p); return Samples(RandomSource, _p);
} }
#endregion
/// <summary> /// <summary>
/// Samples one categorical distributed random variable; also known as the Discrete distribution. /// Samples one categorical distributed random variable; also known as the Discrete distribution.
/// </summary> /// </summary>
@ -395,7 +401,7 @@ namespace MathNet.Numerics.Distributions
// The cumulative density of p. // The cumulative density of p.
var cp = UnnormalizedCdf(p); var cp = UnnormalizedCdf(p);
return DoSample(rnd, cp); return SampleUnchecked(rnd, cp);
} }
/// <summary> /// <summary>
@ -417,47 +423,8 @@ namespace MathNet.Numerics.Distributions
while (true) while (true)
{ {
yield return DoSample(rnd, cp); yield return SampleUnchecked(rnd, cp);
} }
} }
/// <summary>
/// Computes the unnormalized cumulative distribution function. This method performs no
/// parameter checking.
/// </summary>
/// <param name="p">An array of nonnegative ratios: this array does not need to be normalized
/// as this is often impossible using floating point arithmetic.</param>
/// <returns>An array representing the unnormalized cumulative distribution function.</returns>
internal static double[] UnnormalizedCdf(double[] p)
{
var cp = (double[])p.Clone();
for (var i = 1; i < p.Length; i++)
{
cp[i] += cp[i - 1];
}
return cp;
}
/// <summary>
/// Returns one trials from the categorical distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="cdf">The cumulative distribution of the probability distribution.</param>
/// <returns>One sample from the categorical distribution implied by <paramref name="cdf"/>.</returns>
internal static int DoSample(Random rnd, double[] cdf)
{
// TODO : use binary search to speed up this procedure.
var u = rnd.NextDouble() * cdf[cdf.Length - 1];
var idx = 0;
while (u > cdf[idx])
{
idx++;
}
return idx;
}
} }
} }

171
src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs

@ -52,37 +52,37 @@ namespace MathNet.Numerics.Distributions
/// Since many properties of the distribution can only be computed approximately, the tolerance /// Since many properties of the distribution can only be computed approximately, the tolerance
/// level specifies how much error we accept. /// level specifies how much error we accept.
/// </summary> /// </summary>
private const double Tolerance = 1e-12; const double Tolerance = 1e-12;
/// <summary> /// <summary>
/// The mean of the distribution. /// The mean of the distribution.
/// </summary> /// </summary>
private double _mean = double.MinValue; double _mean = double.MinValue;
/// <summary> /// <summary>
/// The variance of the distribution. /// The variance of the distribution.
/// </summary> /// </summary>
private double _variance = double.MinValue; double _variance = double.MinValue;
/// <summary> /// <summary>
/// Caches the value of the normalization constant. /// Caches the value of the normalization constant.
/// </summary> /// </summary>
private double _z = double.MinValue; double _z = double.MinValue;
/// <summary> /// <summary>
/// The lambda parameter. /// The lambda parameter.
/// </summary> /// </summary>
private double _lambda; double _lambda;
/// <summary> /// <summary>
/// The nu parameter. /// The nu parameter.
/// </summary> /// </summary>
private double _nu; double _nu;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the <see cref="ConwayMaxwellPoisson"/> class. /// Initializes a new instance of the <see cref="ConwayMaxwellPoisson"/> class.
@ -105,7 +105,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="lambda">The lambda parameter.</param> /// <param name="lambda">The lambda parameter.</param>
/// <param name="nu">The nu parameter.</param> /// <param name="nu">The nu parameter.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double lambda, double nu) void SetParameters(double lambda, double nu)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu)) if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
{ {
@ -122,7 +122,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="lambda">The lambda parameter.</param> /// <param name="lambda">The lambda parameter.</param>
/// <param name="nu">The nu parameter.</param> /// <param name="nu">The nu parameter.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double lambda, double nu) static bool IsValidParameterSet(double lambda, double nu)
{ {
if (lambda <= 0.0) if (lambda <= 0.0)
{ {
@ -143,15 +143,9 @@ namespace MathNet.Numerics.Distributions
/// <value>The value of the lambda parameter.</value> /// <value>The value of the lambda parameter.</value>
public double Lambda public double Lambda
{ {
get get { return _lambda; }
{
return _lambda;
}
set set { SetParameters(value, _nu); }
{
SetParameters(value, _nu);
}
} }
/// <summary> /// <summary>
@ -160,15 +154,9 @@ namespace MathNet.Numerics.Distributions
/// <value>The value of the Nu parameter.</value> /// <value>The value of the Nu parameter.</value>
public double Nu public double Nu
{ {
get get { return _nu; }
{
return _nu;
}
set set { SetParameters(_lambda, value); }
{
SetParameters(_lambda, value);
}
} }
/// <summary> /// <summary>
@ -189,10 +177,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -333,10 +318,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(Variance); }
{
return Math.Sqrt(Variance);
}
} }
/// <summary> /// <summary>
@ -344,10 +326,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -355,10 +334,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -386,10 +362,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -397,10 +370,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -408,10 +378,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -419,10 +386,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -449,35 +413,12 @@ namespace MathNet.Numerics.Distributions
return Math.Log(Probability(k)); return Math.Log(Probability(k));
} }
/// <summary>
/// Samples a Conway-Maxwell-Poisson distributed random variable.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public int Sample()
{
return DoSample(RandomSource, _lambda, _nu, Z);
}
/// <summary>
/// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables.
/// </summary>
/// <returns>
/// a sequence of samples from a Conway-Maxwell-Poisson distribution.
/// </returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return DoSample(RandomSource, _lambda, _nu, Z);
}
}
#endregion #endregion
/// <summary> /// <summary>
/// Gets the normalization constant of the Conway-Maxwell-Poisson distribution. /// Gets the normalization constant of the Conway-Maxwell-Poisson distribution.
/// </summary> /// </summary>
private double Z double Z
{ {
get get
{ {
@ -488,7 +429,7 @@ namespace MathNet.Numerics.Distributions
_z = Normalization(_lambda, _nu); _z = Normalization(_lambda, _nu);
return _z; return _z;
} }
} }
/// <summary> /// <summary>
@ -499,7 +440,7 @@ namespace MathNet.Numerics.Distributions
/// <returns> /// <returns>
/// an approximate normalization constant for the CMP distribution. /// an approximate normalization constant for the CMP distribution.
/// </returns> /// </returns>
private static double Normalization(double lambda, double nu) static double Normalization(double lambda, double nu)
{ {
// Initialize Z with the first two terms. // Initialize Z with the first two terms.
var z = 1.0 + lambda; var z = 1.0 + lambda;
@ -515,7 +456,7 @@ namespace MathNet.Numerics.Distributions
// The new term. // The new term.
t = t * e; t = t * e;
// The updated normalization constant. // The updated normalization constant.
z = z + t; z = z + t;
@ -542,7 +483,7 @@ namespace MathNet.Numerics.Distributions
/// <returns> /// <returns>
/// One sample from the distribution implied by <paramref name="lambda"/>, <paramref name="nu"/>, and <paramref name="z"/>. /// One sample from the distribution implied by <paramref name="lambda"/>, <paramref name="nu"/>, and <paramref name="z"/>.
/// </returns> /// </returns>
private static int DoSample(Random rnd, double lambda, double nu, double z) internal static int SampleUnchecked(Random rnd, double lambda, double nu, double z)
{ {
var u = rnd.NextDouble(); var u = rnd.NextDouble();
var p = 1.0 / z; var p = 1.0 / z;
@ -558,5 +499,65 @@ namespace MathNet.Numerics.Distributions
return i; return i;
} }
/// <summary>
/// Samples a Conway-Maxwell-Poisson distributed random variable.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public int Sample()
{
return SampleUnchecked(RandomSource, _lambda, _nu, Z);
}
/// <summary>
/// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables.
/// </summary>
/// <returns>
/// a sequence of samples from a Conway-Maxwell-Poisson distribution.
/// </returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return SampleUnchecked(RandomSource, _lambda, _nu, Z);
}
}
/// <summary>
/// Samples a random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The lambda parameter</param>
/// <param name="nu">The nu parameter.</param>
public static int Sample(Random rnd, double lambda, double nu)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
var z = Normalization(lambda, nu);
return SampleUnchecked(rnd, lambda, nu, z);
}
/// <summary>
/// Samples a sequence of this random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The lambda parameter</param>
/// <param name="nu">The nu parameter.</param>
public static IEnumerable<int> Samples(Random rnd, double lambda, double nu)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
var z = Normalization(lambda, nu);
while (true)
{
yield return SampleUnchecked(rnd, lambda, nu, z);
}
}
} }
} }

118
src/Numerics/Distributions/Discrete/DiscreteUniform.cs

@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The distribution's lower bound. /// The distribution's lower bound.
/// </summary> /// </summary>
private int _lower; int _lower;
/// <summary> /// <summary>
/// The distribution's upper bound. /// The distribution's upper bound.
/// </summary> /// </summary>
private int _upper; int _upper;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the DiscreteUniform class. /// Initializes a new instance of the DiscreteUniform class.
@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="lower">Lower bound.</param> /// <param name="lower">Lower bound.</param>
/// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param> /// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(int lower, int upper) static bool IsValidParameterSet(int lower, int upper)
{ {
if (lower <= upper) if (lower <= upper)
{ {
@ -101,7 +101,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="lower">Lower bound.</param> /// <param name="lower">Lower bound.</param>
/// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param> /// <param name="upper">Upper bound; must be at least as large as <paramref name="lower"/>.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(int lower, int upper) void SetParameters(int lower, int upper)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(lower, upper)) if (Control.CheckDistributionParameters && !IsValidParameterSet(lower, upper))
{ {
@ -117,15 +117,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int LowerBound public int LowerBound
{ {
get get { return _lower; }
{
return _lower;
}
set set { SetParameters(value, _upper); }
{
SetParameters(value, _upper);
}
} }
/// <summary> /// <summary>
@ -133,15 +127,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int UpperBound public int UpperBound
{ {
get get { return _upper; }
{
return _upper;
}
set set { SetParameters(_lower, value); }
{
SetParameters(_lower, value);
}
} }
#region IDistribution Members #region IDistribution Members
@ -151,10 +139,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -172,10 +157,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return (_lower + _upper) / 2.0; }
{
return (_lower + _upper) / 2.0;
}
} }
/// <summary> /// <summary>
@ -183,10 +165,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0); }
{
return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0);
}
} }
/// <summary> /// <summary>
@ -194,10 +173,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0; }
{
return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0;
}
} }
/// <summary> /// <summary>
@ -205,10 +181,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { return Math.Log(_upper - _lower + 1.0); }
{
return Math.Log(_upper - _lower + 1.0);
}
} }
/// <summary> /// <summary>
@ -216,10 +189,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return 0.0; }
{
return 0.0;
}
} }
/// <summary> /// <summary>
@ -227,10 +197,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return _lower; }
{
return _lower;
}
} }
/// <summary> /// <summary>
@ -238,10 +205,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return _upper; }
{
return _upper;
}
} }
/// <summary> /// <summary>
@ -255,7 +219,7 @@ namespace MathNet.Numerics.Distributions
{ {
return 0.0; return 0.0;
} }
if (x >= _upper) if (x >= _upper)
{ {
return 1.0; return 1.0;
@ -273,10 +237,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return (int)Math.Floor((_lower + _upper) / 2.0); }
{
return (int)Math.Floor((_lower + _upper) / 2.0);
}
} }
/// <summary> /// <summary>
@ -284,10 +245,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { return (int)Math.Floor((_lower + _upper) / 2.0); }
{
return (int)Math.Floor((_lower + _upper) / 2.0);
}
} }
/// <summary> /// <summary>
@ -324,13 +282,27 @@ namespace MathNet.Numerics.Distributions
return Double.NegativeInfinity; return Double.NegativeInfinity;
} }
#endregion
/// <summary>
/// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking.
/// </summary>
/// <param name="rnd">The random source to use.</param>
/// <param name="lower">The lower bound of the uniform random variable.</param>
/// <param name="upper">The upper bound of the uniform random variable.</param>
/// <returns>A random sample from the discrete uniform distribution.</returns>
internal static int SampleUnchecked(Random rnd, int lower, int upper)
{
return (rnd.Next() % (upper - lower + 1)) + lower;
}
/// <summary> /// <summary>
/// Draws a random sample from the distribution. /// Draws a random sample from the distribution.
/// </summary> /// </summary>
/// <returns>a sample from the distribution.</returns> /// <returns>a sample from the distribution.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _lower, _upper); return SampleUnchecked(RandomSource, _lower, _upper);
} }
/// <summary> /// <summary>
@ -341,12 +313,10 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _lower, _upper); yield return SampleUnchecked(RandomSource, _lower, _upper);
} }
} }
#endregion
/// <summary> /// <summary>
/// Samples a uniformly distributed random variable. /// Samples a uniformly distributed random variable.
/// </summary> /// </summary>
@ -361,7 +331,7 @@ namespace MathNet.Numerics.Distributions
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
} }
return DoSample(rnd, lower, upper); return SampleUnchecked(rnd, lower, upper);
} }
/// <summary> /// <summary>
@ -380,20 +350,8 @@ namespace MathNet.Numerics.Distributions
while (true) while (true)
{ {
yield return DoSample(rnd, lower, upper); yield return SampleUnchecked(rnd, lower, upper);
} }
} }
/// <summary>
/// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking.
/// </summary>
/// <param name="rnd">The random source to use.</param>
/// <param name="lower">The lower bound of the uniform random variable.</param>
/// <param name="upper">The upper bound of the uniform random variable.</param>
/// <returns>A random sample from the discrete uniform distribution.</returns>
private static int DoSample(Random rnd, int lower, int upper)
{
return (rnd.Next() % (upper - lower + 1)) + lower;
}
} }
} }

119
src/Numerics/Distributions/Discrete/Geometric.cs

@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The geometric distribution parameter. /// The geometric distribution parameter.
/// </summary> /// </summary>
private double _p; double _p;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the Geometric class. /// Initializes a new instance of the Geometric class.
@ -68,7 +68,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="p">The probability of generating a one.</param> /// <param name="p">The probability of generating a one.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double p) void SetParameters(double p)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{ {
@ -83,7 +83,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="p">The probability of generating a one.</param> /// <param name="p">The probability of generating a one.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double p) static bool IsValidParameterSet(double p)
{ {
if (p >= 0.0 && p <= 1.0) if (p >= 0.0 && p <= 1.0)
{ {
@ -98,15 +98,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double P public double P
{ {
get get { return _p; }
{
return _p;
}
set set { SetParameters(value); }
{
SetParameters(value);
}
} }
/// <summary> /// <summary>
@ -127,10 +121,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -148,10 +139,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return 1.0 / _p; }
{
return 1.0 / _p;
}
} }
/// <summary> /// <summary>
@ -159,10 +147,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return (1.0 - _p) / (_p * _p); }
{
return (1.0 - _p) / (_p * _p);
}
} }
/// <summary> /// <summary>
@ -170,10 +155,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(1.0 - _p) / _p; }
{
return Math.Sqrt(1.0 - _p) / _p;
}
} }
/// <summary> /// <summary>
@ -181,10 +163,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p; }
{
return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p;
}
} }
/// <summary> /// <summary>
@ -193,10 +172,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>Throws a not supported exception.</remarks> /// <remarks>Throws a not supported exception.</remarks>
public double Skewness public double Skewness
{ {
get get { return (2.0 - _p) / Math.Sqrt(1.0 - _p); }
{
return (2.0 - _p) / Math.Sqrt(1.0 - _p);
}
} }
/// <summary> /// <summary>
@ -218,10 +194,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return 1; }
{
return 1;
}
} }
/// <summary> /// <summary>
@ -229,10 +202,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p)); }
{
return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p));
}
} }
/// <summary> /// <summary>
@ -240,10 +210,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 1; }
{
return 1;
}
} }
/// <summary> /// <summary>
@ -251,10 +218,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return int.MaxValue; }
{
return int.MaxValue;
}
} }
/// <summary> /// <summary>
@ -291,13 +255,28 @@ namespace MathNet.Numerics.Distributions
return ((k - 1) * Math.Log(1.0 - _p)) + Math.Log(_p); return ((k - 1) * Math.Log(1.0 - _p)) + Math.Log(_p);
} }
#endregion
/// <summary>
/// Returns one sample from the distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The p parameter</param>
/// <returns>
/// One sample from the distribution implied by <paramref name="p"/>.
/// </returns>
internal static int SampleUnchecked(Random rnd, double p)
{
return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p));
}
/// <summary> /// <summary>
/// Samples a Geometric distributed random variable. /// Samples a Geometric distributed random variable.
/// </summary> /// </summary>
/// <returns>A sample from the Geometric distribution.</returns> /// <returns>A sample from the Geometric distribution.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _p); return SampleUnchecked(RandomSource, _p);
} }
/// <summary> /// <summary>
@ -308,23 +287,41 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _p); yield return SampleUnchecked(RandomSource, _p);
} }
} }
#endregion /// <summary>
/// Samples a random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The p parameter</param>
public static int Sample(Random rnd, double p)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, p);
}
/// <summary> /// <summary>
/// Returns one sample from the distribution. /// Samples a sequence of this random variable.
/// </summary> /// </summary>
/// <param name="rnd">The random number generator to use.</param> /// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The p parameter</param> /// <param name="p">The p parameter</param>
/// <returns> public static IEnumerable<int> Samples(Random rnd, double p)
/// One sample from the distribution implied by <paramref name="p"/>.
/// </returns>
private static int DoSample(Random rnd, double p)
{ {
return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p)); if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
while (true)
{
yield return SampleUnchecked(rnd, p);
}
} }
} }
} }

115
src/Numerics/Distributions/Discrete/Hypergeometric.cs

@ -52,22 +52,22 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The size of the population. /// The size of the population.
/// </summary> /// </summary>
private int _populationSize; int _populationSize;
/// <summary> /// <summary>
/// The m parameter of the distribution. /// The m parameter of the distribution.
/// </summary> /// </summary>
private int _m; int _m;
/// <summary> /// <summary>
/// The n parameter (number to draw) of the distribution. /// The n parameter (number to draw) of the distribution.
/// </summary> /// </summary>
private int _n; int _n;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the Hypergeometric class. /// Initializes a new instance of the Hypergeometric class.
@ -87,7 +87,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="total">The Total parameter of the distribution.</param> /// <param name="total">The Total parameter of the distribution.</param>
/// <param name="m">The m parameter of the distribution.</param> /// <param name="m">The m parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
private void SetParameters(int total, int m, int n) void SetParameters(int total, int m, int n)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(total, m, n)) if (Control.CheckDistributionParameters && !IsValidParameterSet(total, m, n))
{ {
@ -106,7 +106,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="m">The m parameter of the distribution.</param> /// <param name="m">The m parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(int total, int m, int n) static bool IsValidParameterSet(int total, int m, int n)
{ {
if (total < 0 || m < 0 || n < 0) if (total < 0 || m < 0 || n < 0)
{ {
@ -166,10 +166,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -187,10 +184,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return (double)_m * _n / _populationSize; }
{
return (double)_m * _n / _populationSize;
}
} }
/// <summary> /// <summary>
@ -198,10 +192,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0)); }
{
return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0));
}
} }
/// <summary> /// <summary>
@ -225,10 +216,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0)); }
{
return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0));
}
} }
/// <summary> /// <summary>
@ -320,27 +308,6 @@ namespace MathNet.Numerics.Distributions
return Math.Log(Probability(k)); return Math.Log(Probability(k));
} }
/// <summary>
/// Samples a Hypergeometric distributed random variable.
/// </summary>
/// <returns>The number of successes in n trials.</returns>
public int Sample()
{
return DoSample(RandomSource, _populationSize, _m, _n);
}
/// <summary>
/// Samples an array of Hypergeometric distributed random variables.
/// </summary>
/// <returns>a sequence of successes in n trials.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return DoSample(RandomSource, _populationSize, _m, _n);
}
}
#endregion #endregion
/// <summary> /// <summary>
@ -351,7 +318,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="m">The m parameter of the distribution.</param> /// <param name="m">The m parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
/// <returns>a random number from the Hypergeometric distribution.</returns> /// <returns>a random number from the Hypergeometric distribution.</returns>
private static int DoSample(Random rnd, int size, int m, int n) internal static int SampleUnchecked(Random rnd, int size, int m, int n)
{ {
var x = 0; var x = 0;
@ -367,10 +334,68 @@ namespace MathNet.Numerics.Distributions
size--; size--;
n--; n--;
} }
while (0 < n); while (0 < n);
return x; return x;
} }
/// <summary>
/// Samples a Hypergeometric distributed random variable.
/// </summary>
/// <returns>The number of successes in n trials.</returns>
public int Sample()
{
return SampleUnchecked(RandomSource, _populationSize, _m, _n);
}
/// <summary>
/// Samples an array of Hypergeometric distributed random variables.
/// </summary>
/// <returns>a sequence of successes in n trials.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return SampleUnchecked(RandomSource, _populationSize, _m, _n);
}
}
/// <summary>
/// Samples a random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="populationSize">The population size.</param>
/// <param name="m">The m parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param>
public static int Sample(Random rnd, int populationSize, int m, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, populationSize, m, n);
}
/// <summary>
/// Samples a sequence of this random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="populationSize">The population size.</param>
/// <param name="m">The m parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param>
public static IEnumerable<int> Samples(Random rnd, int populationSize, int m, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
while (true)
{
yield return SampleUnchecked(rnd, populationSize, m, n);
}
}
} }
} }

164
src/Numerics/Distributions/Discrete/NegativeBinomial.cs

@ -46,32 +46,26 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The r parameter of the distribution. /// The r parameter of the distribution.
/// </summary> /// </summary>
private double _r; double _r;
/// <summary> /// <summary>
/// The p parameter of the distribution. /// The p parameter of the distribution.
/// </summary> /// </summary>
private double _p; double _p;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Gets or sets the number of trials. /// Gets or sets the number of trials.
/// </summary> /// </summary>
public double R public double R
{ {
get get { return _r; }
{
return _r;
}
set set { SetParameters(value, _p); }
{
SetParameters(value, _p);
}
} }
/// <summary> /// <summary>
@ -79,15 +73,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double P public double P
{ {
get get { return _p; }
{
return _p;
}
set set { SetParameters(_r, value); }
{
SetParameters(_r, value);
}
} }
/// <summary> /// <summary>
@ -107,7 +95,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="r">The number of trials.</param> /// <param name="r">The number of trials.</param>
/// <param name="p">The probability of a trial resulting in success.</param> /// <param name="p">The probability of a trial resulting in success.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double r, double p) void SetParameters(double r, double p)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p)) if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
{ {
@ -124,7 +112,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="r">The number of trials.</param> /// <param name="r">The number of trials.</param>
/// <param name="p">The probability of a trial resulting in success.</param> /// <param name="p">The probability of a trial resulting in success.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double r, double p) static bool IsValidParameterSet(double r, double p)
{ {
if (r < 0.0 || Double.IsNaN(r)) if (r < 0.0 || Double.IsNaN(r))
{ {
@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return _r * (1.0 - _p) / _p; }
{
return _r * (1.0 - _p) / _p;
}
} }
/// <summary> /// <summary>
@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _r * (1.0 - _p) / (_p * _p); }
{
return _r * (1.0 - _p) / (_p * _p);
}
} }
/// <summary> /// <summary>
@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(_r * (1.0 - _p)) / _p; }
{
return Math.Sqrt(_r * (1.0 - _p)) / _p;
}
} }
/// <summary> /// <summary>
@ -211,10 +187,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Entropy public double Entropy
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -222,10 +195,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p)); }
{
return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p));
}
} }
/// <summary> /// <summary>
@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0; }
{
return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0;
}
} }
/// <summary> /// <summary>
@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -269,10 +233,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -280,10 +241,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return int.MaxValue; }
{
return int.MaxValue;
}
} }
/// <summary> /// <summary>
@ -296,10 +254,10 @@ namespace MathNet.Numerics.Distributions
public double Probability(int k) public double Probability(int k)
{ {
var ln = SpecialFunctions.GammaLn(_r + k) var ln = SpecialFunctions.GammaLn(_r + k)
- SpecialFunctions.GammaLn(_r) - SpecialFunctions.GammaLn(_r)
- SpecialFunctions.GammaLn(k + 1.0) - SpecialFunctions.GammaLn(k + 1.0)
+ (_r * Math.Log(_p)) + (_r * Math.Log(_p))
+ (k * Math.Log(1.0 - _p)); + (k * Math.Log(1.0 - _p));
return Math.Exp(ln); return Math.Exp(ln);
} }
@ -313,20 +271,44 @@ namespace MathNet.Numerics.Distributions
public double ProbabilityLn(int k) public double ProbabilityLn(int k)
{ {
var ln = SpecialFunctions.GammaLn(_r + k) var ln = SpecialFunctions.GammaLn(_r + k)
- SpecialFunctions.GammaLn(_r) - SpecialFunctions.GammaLn(_r)
- SpecialFunctions.GammaLn(k + 1.0) - SpecialFunctions.GammaLn(k + 1.0)
+ (_r * Math.Log(_p)) + (_r * Math.Log(_p))
+ (k * Math.Log(1.0 - _p)); + (k * Math.Log(1.0 - _p));
return ln; return ln;
} }
#endregion
/// <summary>
/// Samples a negative binomial distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="r">The r parameter.</param>
/// <param name="p">The p parameter.</param>
/// <returns>a sample from the distribution.</returns>
internal static int SampleUnchecked(Random rnd, double r, double p)
{
var lambda = Gamma.SampleUnchecked(rnd, r, p);
var c = Math.Exp(-lambda);
var p1 = 1.0;
var k = 0;
do
{
k = k + 1;
p1 = p1 * rnd.NextDouble();
}
while (p1 >= c);
return k - 1;
}
/// <summary> /// <summary>
/// Samples a <c>NegativeBinomial</c> distributed random variable. /// Samples a <c>NegativeBinomial</c> distributed random variable.
/// </summary> /// </summary>
/// <returns>a sample from the distribution.</returns> /// <returns>a sample from the distribution.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _r, _p); return SampleUnchecked(RandomSource, _r, _p);
} }
/// <summary> /// <summary>
@ -337,43 +319,43 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _r, _p); yield return SampleUnchecked(RandomSource, _r, _p);
} }
} }
#endregion
/// <summary> /// <summary>
/// Samples a negative binomial distributed random variable. /// Samples a random variable.
/// </summary> /// </summary>
/// <param name="rnd">The random number generator to use.</param> /// <param name="rnd">The random number generator to use.</param>
/// <param name="r">The r parameter.</param> /// <param name="r">The r parameter.</param>
/// <param name="p">The p parameter.</param> /// <param name="p">The p parameter.</param>
/// <returns>a sample from the distribution.</returns> public static int Sample(Random rnd, double r, double p)
private static int DoSample(Random rnd, double r, double p)
{ {
var lambda = Gamma.Sample(rnd, r, p); if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
return DoPoissonSample(rnd, lambda); {
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, r, p);
} }
/// <summary> /// <summary>
/// Use Knuth's method to generate Poisson distributed random variables. /// Samples a sequence of this random variable.
/// </summary> /// </summary>
/// <param name="rnd">The random number generator.</param> /// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The lambda value.</param> /// <param name="r">The r parameter.</param>
/// <returns>a Poisson distributed random number.</returns> /// <param name="p">The p parameter.</param>
private static int DoPoissonSample(Random rnd, double lambda) public static IEnumerable<int> Samples(Random rnd, double r, double p)
{ {
var c = Math.Exp(-lambda); if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p))
var p = 1.0;
var k = 0;
do
{ {
k = k + 1; throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
p = p * rnd.NextDouble(); }
while (true)
{
yield return SampleUnchecked(rnd, r, p);
} }
while (p >= c);
return k - 1;
} }
} }
} }

186
src/Numerics/Distributions/Discrete/Poisson.cs

@ -43,27 +43,21 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The Poisson distribution parameter λ. /// The Poisson distribution parameter λ.
/// </summary> /// </summary>
private double _lambda; double _lambda;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Gets or sets the Poisson distribution parameter λ. /// Gets or sets the Poisson distribution parameter λ.
/// </summary> /// </summary>
public double Lambda public double Lambda
{ {
get get { return _lambda; }
{
return _lambda;
}
set set { SetParameters(value); }
{
SetParameters(value);
}
} }
/// <summary> /// <summary>
@ -82,7 +76,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="lambda">The mean (λ) of the distribution.</param> /// <param name="lambda">The mean (λ) of the distribution.</param>
/// <exception cref="System.ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception> /// <exception cref="System.ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double lambda) void SetParameters(double lambda)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
{ {
@ -97,7 +91,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="lambda">The mean (λ) of the distribution.</param> /// <param name="lambda">The mean (λ) of the distribution.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double lambda) static bool IsValidParameterSet(double lambda)
{ {
return lambda > 0.00; return lambda > 0.00;
} }
@ -120,10 +114,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -141,10 +132,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return _lambda; }
{
return _lambda;
}
} }
/// <summary> /// <summary>
@ -152,10 +140,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Variance public double Variance
{ {
get get { return _lambda; }
{
return _lambda;
}
} }
/// <summary> /// <summary>
@ -163,10 +148,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(_lambda); }
{
return Math.Sqrt(_lambda);
}
} }
/// <summary> /// <summary>
@ -175,10 +157,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks> /// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
public double Entropy public double Entropy
{ {
get get { return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda)); }
{
return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda));
}
} }
/// <summary> /// <summary>
@ -186,10 +165,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Skewness public double Skewness
{ {
get get { return 1.0 / Math.Sqrt(_lambda); }
{
return 1.0 / Math.Sqrt(_lambda);
}
} }
/// <summary> /// <summary>
@ -197,10 +173,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 0; }
{
return 0;
}
} }
/// <summary> /// <summary>
@ -208,10 +181,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return int.MaxValue; }
{
return int.MaxValue;
}
} }
/// <summary> /// <summary>
@ -233,10 +203,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return (int)Math.Floor(_lambda); }
{
return (int)Math.Floor(_lambda);
}
} }
/// <summary> /// <summary>
@ -245,10 +212,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks> /// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
public int Median public int Median
{ {
get get { return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda)); }
{
return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda));
}
} }
/// <summary> /// <summary>
@ -271,71 +235,15 @@ namespace MathNet.Numerics.Distributions
return -_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k); return -_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k);
} }
/// <summary>
/// Samples a Poisson distributed random variable.
/// </summary>
/// <returns>A sample from the Poisson distribution.</returns>
public int Sample()
{
return DoSample(RandomSource, _lambda);
}
/// <summary>
/// Samples an array of Poisson distributed random variables.
/// </summary>
/// <returns>a sequence of successes in N trials.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return DoSample(RandomSource, _lambda);
}
}
#endregion #endregion
/// <summary>
/// Samples a Poisson distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>A sample from the Poisson distribution.</returns>
public static int Sample(Random rnd, double lambda)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return DoSample(rnd, lambda);
}
/// <summary>
/// Samples a sequence of Poisson distributed random variables.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<int> Samples(Random rnd, double lambda)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
while (true)
{
yield return DoSample(rnd, lambda);
}
}
/// <summary> /// <summary>
/// Generates one sample from the Poisson distribution. /// Generates one sample from the Poisson distribution.
/// </summary> /// </summary>
/// <param name="rnd">The random source to use.</param> /// <param name="rnd">The random source to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param> /// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>A random sample from the Poisson distribution.</returns> /// <returns>A random sample from the Poisson distribution.</returns>
private static int DoSample(Random rnd, double lambda) internal static int SampleUnchecked(Random rnd, double lambda)
{ {
return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda); return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda);
} }
@ -346,7 +254,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="rnd">The random source to use.</param> /// <param name="rnd">The random source to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param> /// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>A random sample from the Poisson distribution.</returns> /// <returns>A random sample from the Poisson distribution.</returns>
private static int DoSampleShort(Random rnd, double lambda) static int DoSampleShort(Random rnd, double lambda)
{ {
var limit = Math.Exp(-lambda); var limit = Math.Exp(-lambda);
var count = 0; var count = 0;
@ -367,7 +275,7 @@ namespace MathNet.Numerics.Distributions
/// <remarks>"Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson, /// <remarks>"Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson,
/// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979) /// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
/// The article is on pages 29-35. The algorithm given here is on page 32. </remarks> /// The article is on pages 29-35. The algorithm given here is on page 32. </remarks>
private static int DoSampleLarge(Random rnd, double lambda) static int DoSampleLarge(Random rnd, double lambda)
{ {
var c = 0.767 - (3.36 / lambda); var c = 0.767 - (3.36 / lambda);
var beta = Math.PI / Math.Sqrt(3.0 * lambda); var beta = Math.PI / Math.Sqrt(3.0 * lambda);
@ -395,5 +303,61 @@ namespace MathNet.Numerics.Distributions
} }
} }
} }
/// <summary>
/// Samples a Poisson distributed random variable.
/// </summary>
/// <returns>A sample from the Poisson distribution.</returns>
public int Sample()
{
return SampleUnchecked(RandomSource, _lambda);
}
/// <summary>
/// Samples an array of Poisson distributed random variables.
/// </summary>
/// <returns>a sequence of successes in N trials.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return SampleUnchecked(RandomSource, _lambda);
}
}
/// <summary>
/// Samples a Poisson distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>A sample from the Poisson distribution.</returns>
public static int Sample(Random rnd, double lambda)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, lambda);
}
/// <summary>
/// Samples a sequence of Poisson distributed random variables.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lambda">The Poisson distribution parameter λ.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<int> Samples(Random rnd, double lambda)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
while (true)
{
yield return SampleUnchecked(rnd, lambda);
}
}
} }
} }

142
src/Numerics/Distributions/Discrete/Zipf.cs

@ -47,17 +47,17 @@ namespace MathNet.Numerics.Distributions
/// <summary> /// <summary>
/// The s parameter of the distribution. /// The s parameter of the distribution.
/// </summary> /// </summary>
private double _s; double _s;
/// <summary> /// <summary>
/// The n parameter of the distribution. /// The n parameter of the distribution.
/// </summary> /// </summary>
private int _n; int _n;
/// <summary> /// <summary>
/// The distribution's random number generator. /// The distribution's random number generator.
/// </summary> /// </summary>
private Random _random; Random _random;
/// <summary> /// <summary>
/// Initializes a new instance of the <see cref="Zipf"/> class. /// Initializes a new instance of the <see cref="Zipf"/> class.
@ -79,7 +79,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
/// <param name="s">The s parameter of the distribution.</param> /// <param name="s">The s parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
private void SetParameters(double s, int n) void SetParameters(double s, int n)
{ {
if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n)) if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
{ {
@ -96,7 +96,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="s">The s parameter of the distribution.</param> /// <param name="s">The s parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns> /// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
private static bool IsValidParameterSet(double s, int n) static bool IsValidParameterSet(double s, int n)
{ {
if (n <= 0 || s <= 0 || Double.IsNaN(s)) if (n <= 0 || s <= 0 || Double.IsNaN(s))
{ {
@ -111,15 +111,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double S public double S
{ {
get get { return _s; }
{
return _s;
}
set set { SetParameters(value, _n); }
{
SetParameters(value, _n);
}
} }
/// <summary> /// <summary>
@ -127,15 +121,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int N public int N
{ {
get get { return _n; }
{
return _n;
}
set set { SetParameters(_s, value); }
{
SetParameters(_s, value);
}
} }
/// <summary> /// <summary>
@ -154,10 +142,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public Random RandomSource public Random RandomSource
{ {
get get { return _random; }
{
return _random;
}
set set
{ {
@ -175,10 +160,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double Mean public double Mean
{ {
get get { return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s); }
{
return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s);
}
} }
/// <summary> /// <summary>
@ -203,10 +185,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public double StdDev public double StdDev
{ {
get get { return Math.Sqrt(Variance); }
{
return Math.Sqrt(Variance);
}
} }
/// <summary> /// <summary>
@ -266,10 +245,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Mode public int Mode
{ {
get get { return 1; }
{
return 1;
}
} }
/// <summary> /// <summary>
@ -277,10 +253,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Median public int Median
{ {
get get { throw new NotSupportedException(); }
{
throw new NotSupportedException();
}
} }
/// <summary> /// <summary>
@ -288,10 +261,7 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Minimum public int Minimum
{ {
get get { return 1; }
{
return 1;
}
} }
/// <summary> /// <summary>
@ -299,12 +269,9 @@ namespace MathNet.Numerics.Distributions
/// </summary> /// </summary>
public int Maximum public int Maximum
{ {
get get { return _n; }
{
return _n;
}
} }
/// <summary> /// <summary>
/// Computes values of the probability mass function. /// Computes values of the probability mass function.
/// </summary> /// </summary>
@ -329,13 +296,45 @@ namespace MathNet.Numerics.Distributions
return Math.Log(Probability(k)); return Math.Log(Probability(k));
} }
#endregion
/// <summary>
/// Generates a sample from the Zipf distribution without doing parameter checking.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="s">The s parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param>
/// <returns>a random number from the Zipf distribution.</returns>
internal static int SampleUnchecked(Random rnd, double s, int n)
{
var r = 0.0;
while (r == 0.0)
{
r = rnd.NextDouble();
}
var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s);
int i;
var sum = 0.0;
for (i = 1; i <= n; i++)
{
sum += p / Math.Pow(i, s);
if (sum >= r)
{
break;
}
}
return i;
}
/// <summary> /// <summary>
/// Draws a random sample from the distribution. /// Draws a random sample from the distribution.
/// </summary> /// </summary>
/// <returns>a sample from the distribution.</returns> /// <returns>a sample from the distribution.</returns>
public int Sample() public int Sample()
{ {
return DoSample(RandomSource, _s, _n); return SampleUnchecked(RandomSource, _s, _n);
} }
/// <summary> /// <summary>
@ -346,40 +345,43 @@ namespace MathNet.Numerics.Distributions
{ {
while (true) while (true)
{ {
yield return DoSample(RandomSource, _s, _n); yield return SampleUnchecked(RandomSource, _s, _n);
} }
} }
#endregion
/// <summary> /// <summary>
/// Generates a sample from the Zipf distribution without doing parameter checking. /// Samples a random variable.
/// </summary> /// </summary>
/// <param name="rnd">The random number generator to use.</param> /// <param name="rnd">The random number generator to use.</param>
/// <param name="s">The s parameter of the distribution.</param> /// <param name="s">The s parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param> /// <param name="n">The n parameter of the distribution.</param>
/// <returns>a random number from the Zipf distribution.</returns> public static int Sample(Random rnd, double s, int n)
private static int DoSample(Random rnd, double s, int n)
{ {
var r = 0.0; if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
while (r == 0.0)
{ {
r = rnd.NextDouble(); throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
} }
var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s); return SampleUnchecked(rnd, s, n);
int i; }
var sum = 0.0;
for (i = 1; i <= n; i++) /// <summary>
/// Samples a sequence of this random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="s">The s parameter of the distribution.</param>
/// <param name="n">The n parameter of the distribution.</param>
public static IEnumerable<int> Samples(Random rnd, double s, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n))
{ {
sum += p / Math.Pow(i, s); throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
if (sum >= r)
{
break;
}
} }
return i; while (true)
{
yield return SampleUnchecked(rnd, s, n);
}
} }
} }
} }

4
src/Numerics/Distributions/Multivariate/Multinomial.cs

@ -368,7 +368,7 @@ namespace MathNet.Numerics.Distributions
for (var i = 0; i < n; i++) for (var i = 0; i < n; i++)
{ {
ret[Categorical.DoSample(rnd, cp)]++; ret[Categorical.SampleUnchecked(rnd, cp)]++;
} }
return ret; return ret;
@ -399,7 +399,7 @@ namespace MathNet.Numerics.Distributions
for (var i = 0; i < n; i++) for (var i = 0; i < n; i++)
{ {
ret[Categorical.DoSample(rnd, cp)]++; ret[Categorical.SampleUnchecked(rnd, cp)]++;
} }
yield return ret; yield return ret;

Loading…
Cancel
Save