Browse Source

Some bug fixes on Binomial, Categorical and Multinomial distributions.

la-knuth
Jurgen Van Gael 17 years ago
parent
commit
1ce3b42491
  1. 391
      src/Numerics/Distributions/Discrete/Binomial.cs
  2. 422
      src/Numerics/Distributions/Discrete/Categorical.cs
  3. 140
      src/Numerics/Distributions/Multivariate/Multinomial.cs
  4. 2
      src/Numerics/Numerics.csproj
  5. 4
      src/UnitTests/DistributionTests/CommonDistributionTests.cs
  6. 248
      src/UnitTests/DistributionTests/Discrete/BinomialTests.cs
  7. 117
      src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs
  8. 18
      src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs
  9. 2
      src/UnitTests/UnitTests.csproj

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

@ -0,0 +1,391 @@
// <copyright file="Binomial.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.Distributions
{
using System;
using System.Collections.Generic;
using Properties;
/// <summary>
/// Implements the binomial distribution. For details about this distribution, see
/// <a href="http://en.wikipedia.org/wiki/Binomial_distribution">Wikipedia - Binomial distribution</a>.
/// </summary>
/// <remarks><para>The distribution is parameterized by a probability (between 0.0 and 1.0).</para>
/// <para>The distribution will use the <see cref="System.Random"/> by default.
/// Users can set the random number generator by using the <see cref="RandomSource"/> property.</para>
/// <para>The statistics classes will check all the incoming parameters whether they are in the allowed
/// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
/// to false, all parameter checks can be turned off.</para></remarks>
public class Binomial : IDiscreteDistribution
{
/// <summary>
/// Stores the normalized binomial probability.
/// </summary>
private double _p;
/// <summary>
/// The number of trials.
/// </summary>
private int _n;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// Initializes a new instance of the Binomial class.
/// </summary>
/// <param name="p">The success probability of a trial.</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="n"/> is negative.</exception>
public Binomial(double p, int n)
{
SetParameters(p, n);
RandomSource = new System.Random();
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
public override string ToString()
{
return "Binomial(Success Probability = " + _p + ", Number of Trials = " + _n + ")";
}
/// <summary>
/// Checks whether the parameters of the distribution are valid.
/// </summary>
/// <param name="p">The success probability of a trial.</param>
/// <param name="n">The number of trials.</param>
/// <returns>false <paramref name="p"/> is not in the interval [0.0,1.0] or <paramref name="n"/> is negative, true otherwise.</exception>
private static bool IsValidParameterSet(double p, int n)
{
if(p < 0.0 || p > 1.0)
{
return false;
}
if(n < 0)
{
return false;
}
return true;
}
/// <summary>
/// Sets the parameters of the distribution after checking their validity.
/// </summary>
/// <param name="p">The success probability of a trial.</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="n"/> is negative.</exception>
private void SetParameters(double p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
_p = p;
_n = n;
}
/// <summary>
/// Gets or sets the success probability.
/// </summary>
public double P
{
get
{
return _p;
}
set
{
SetParameters(value, _n);
}
}
/// <summary>
/// Gets or sets the number of trials.
/// </summary>
public int N
{
get
{
return _n;
}
set
{
SetParameters(_p, value);
}
}
#region IDistribution Members
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// Gets the mean of the distribution.
/// </summary>
public double Mean
{
get { return _p * _n; }
}
/// <summary>
/// Gets the standard deviation of the distribution.
/// </summary>
public double StdDev
{
get { return Math.Sqrt(_p * (1.0 - _p) * _n); }
}
/// <summary>
/// Gets the variance of the distribution.
/// </summary>
public double Variance
{
get { return _p * (1.0 - _p) * _n; }
}
/// <summary>
/// Gets the entropy of the distribution.
/// </summary>
public double Entropy
{
get
{
double E = 0.0;
for(int i = 0; i < _n; i++)
{
double p = Probability(i);
E += p * Math.Log(p);
}
return E;
}
}
/// <summary>
/// Gets the skewness of the distribution.
/// </summary>
public double Skewness
{
get { return (1.0 - 2.0 * _p) / Math.Sqrt(_n * _p * (1.0 - _p)); }
}
/// <summary>
/// Gets the smallest element in the domain of the distributions which can be represented by an integer.
/// </summary>
public int Minimum { get { return 0; } }
/// <summary>
/// Gets the largest element in the domain of the distributions which can be represented by an integer.
/// </summary>
public int Maximum { get { return _n; } }
/// <summary>
/// Computes the cumulative distribution function of the Binomial distribution.
/// </summary>
/// <param name="x">The location at which to compute the cumulative density.</param>
/// <returns>the cumulative density at <paramref name="x"/>.</returns>
public double CumulativeDistribution(double x)
{
if (x < 0.0)
{
return 0.0;
}
else if (x > _n)
{
return 1.0;
}
int k = (int) Math.Floor(x);
return (_n - k) * Combinatorics.Combinations(_n,k) * SpecialFunctions.BetaRegularized(_n - k, 1 + k, 1-_p);
}
#endregion
#region IDiscreteDistribution Members
/// <summary>
/// The mode of the distribution.
/// </summary>
public int Mode
{
get { return (int) Math.Floor((_n + 1) * _p); }
}
/// <summary>
/// The median of the distribution.
/// </summary>
public int Median
{
get { throw new NotImplementedException(); }
}
/// <summary>
/// Computes the probability of a specific value.
/// </summary>
public double Probability(int val)
{
if (val < 0)
{
return 0.0;
}
if (val > _n)
{
return 0.0;
}
return SpecialFunctions.Binomial(_n, val) * Math.Pow(_p, val) * Math.Pow(1.0 - _p, _n - val);
}
/// <summary>
/// Computes the probability of a specific value.
/// </summary>
public double ProbabilityLn(int val)
{
if (val < 0)
{
return 0.0;
}
if (val > _n)
{
return 0.0;
}
return SpecialFunctions.BinomialLn(_n, val) + val * Math.Log(_p) + (_n - val) * Math.Log(1.0 - _p);
}
/// <summary>
/// Samples a Binomially distributed random variable.
/// </summary>
/// <returns>The number of successful trials.</returns>
public int Sample()
{
return DoSample(RandomSource, _p, _n);
}
/// <summary>
/// Samples an array of Bernoulli distributed random variables.
/// </summary>
/// <returns>a sequence of successful trial counts.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return DoSample(RandomSource, _p, _n);
}
}
#endregion
/// <summary>
/// Samples a binomially distributed random variable.
/// </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 successes in <see cref="N"/> trials.</returns>
public static int Sample(System.Random rnd, double p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
return DoSample(rnd, p, n);
}
/// <summary>
/// Samples a sequence of binomially distributed random variable.
/// </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>a sequence of successful trial counts.</returns>
public static IEnumerable<int> Samples(System.Random rnd, double p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
while (true)
{
yield return DoSample(rnd, p, n);
}
}
/// <summary>
/// Generates a sample from the Binomial distribution without doing parameter checking.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="p">The success probability of a trial; must be in the interval [0.0, 1.0].</param>
/// <param name="n">The number of trials; must be positive.</param>
/// <returns>The number of successful trials.</returns>
private static int DoSample(System.Random rnd, double p, int n)
{
int k = 0;
for (int i = 0; i < n; i++)
{
k += (rnd.NextDouble() < p ? 1 : 0);
}
return k;
}
}
}

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

@ -0,0 +1,422 @@
// <copyright file="Categorical.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.Distributions
{
using System;
using System.Collections.Generic;
using Properties;
using MathNet.Numerics.Statistics;
/// <summary>
/// Implements the categorical distribution. For details about this distribution, see
/// <a href="http://en.wikipedia.org/wiki/Categorical_distribution">Wikipedia - Categorical distribution</a>. This
/// distribution is sometimes called the Discrete distribution.
/// </summary>
/// <remarks><para>The distribution is parameterized by a vector of ratios: in other words, the parameter
/// does not have to be normalized and sum to 1. The reason is that some vectors can't be exactly normalized
/// to sum to 1 in floating point representation.</para>
/// <para>The distribution will use the <see cref="System.Random"/> by default.
/// Users can set the random number generator by using the <see cref="RandomSource"/> property.</para>
/// <para>The statistics classes will check all the incoming parameters whether they are in the allowed
/// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
/// to false, all parameter checks can be turned off.</para></remarks>
public class Categorical : IDiscreteDistribution
{
/// <summary>
/// Stores the normalized categorical probabilities.
/// </summary>
private double[] _p;
/// <summary>
/// The distribution's random number generator.
/// </summary>
private Random _random;
/// <summary>
/// Initializes a new instance of the Categorical class.
/// </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>
/// <exception cref="ArgumentException">If any of the probabilities are negative or do not sum to one.</exception>
public Categorical(double[] p)
{
SetParameters(p);
RandomSource = new System.Random();
}
/* TODO
/// <summary>
/// Generate a categorical distribution from histogram <paramref name="h"/>. The distribution will
/// not be automatically updated when the histogram changes.
/// </summary>
public Categorical(Histogram h)
{
// The probability distribution vector.
_p = new double[h.BinCount];
// Fill in the distribution vector.
for (int i = 0; i < h.BinCount; i++)
{
_p[i] = h[i];
}
RandomNumberGenerator = new System.Random();
}*/
/// <summary>
/// A string representation of the distribution.
/// </summary>
public override string ToString()
{
return "Categorical(Dimension = " + _p.Length + ")";
}
/// <summary>
/// Checks whether the parameters of the distribution are valid.
/// </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>If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true</returns>
private static bool IsValidParameterSet(double[] p)
{
double sum = 0.0;
for (int i = 0; i < p.Length; i++)
{
if (p[i] < 0.0 || Double.IsNaN(p[i]))
{
return false;
}
else
{
sum += p[i];
}
}
if (sum == 0.0)
{
return false;
}
return true;
}
/// <summary>
/// Sets the parameters of the distribution after checking their validity.
/// </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>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double[] p)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
_p = (double[])p.Clone();
}
/// <summary>
/// Gets or sets the probability of generating a one.
/// </summary>
public double[] P
{
get
{
return (double[]) _p.Clone();
}
set
{
SetParameters(value);
}
}
#region IDistribution Members
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public Random RandomSource
{
get
{
return _random;
}
set
{
if (value == null)
{
throw new ArgumentNullException();
}
_random = value;
}
}
/// <summary>
/// Gets the mean of the distribution.
/// </summary>
public double Mean
{
get { return _p.Mean(); }
}
/// <summary>
/// Gets the standard deviation of the distribution.
/// </summary>
public double StdDev
{
get { return _p.StandardDeviation(); }
}
/// <summary>
/// Gets the variance of the distribution.
/// </summary>
public double Variance
{
get { return _p.Variance(); }
}
/// <summary>
/// Gets the entropy of the distribution.
/// </summary>
public double Entropy
{
get {
double E = 0.0;
for (int i = 0; i < _p.Length; i++)
{
double p = _p[i];
E += p * Math.Log(p);
}
return E;
}
}
/// <summary>
/// Gets the skewness of the distribution.
/// </summary>
public double Skewness
{
get { throw new NotImplementedException(); }
}
/// <summary>
/// Gets the smallest element in the domain of the distributions which can be represented by an integer.
/// </summary>
public int Minimum { get { return 0; } }
/// <summary>
/// Gets the largest element in the domain of the distributions which can be represented by an integer.
/// </summary>
public int Maximum { get { return _p.Length-1; } }
/// <summary>
/// Computes the cumulative distribution function of the Binomial distribution.
/// </summary>
/// <param name="x">The location at which to compute the cumulative density.</param>
/// <returns>the cumulative density at <paramref name="x"/>.</returns>
public double CumulativeDistribution(double x)
{
if (x < 0.0)
{
return 0.0;
}
else if (x >= _p.Length)
{
return 1.0;
}
var cdf = UnnormalizedCDF(_p);
return cdf[(int) Math.Floor(x)] / cdf[_p.Length - 1];
}
#endregion
#region IDiscreteDistribution Members
/// <summary>
/// The mode of the distribution.
/// </summary>
public int Mode
{
get { throw new NotImplementedException(); }
}
/// <summary>
/// The median of the distribution.
/// </summary>
public int Median
{
get { return (int) _p.Median(); }
}
/// <summary>
/// Computes the probability of a specific value.
/// </summary>
public double Probability(int val)
{
if (val < 0)
{
return 0.0;
}
if (val >= _p.Length)
{
return 0.0;
}
return _p[val];
}
/// <summary>
/// Computes the probability of a specific value.
/// </summary>
public double ProbabilityLn(int val)
{
if (val < 0)
{
return 0.0;
}
if (val >= _p.Length)
{
return 0.0;
}
return Math.Log(_p[val]);
}
/// <summary>
/// Samples a Binomially distributed random variable.
/// </summary>
/// <returns>The number of successful trials.</returns>
public int Sample()
{
return DoSample(RandomSource, _p);
}
/// <summary>
/// Samples an array of Bernoulli distributed random variables.
/// </summary>
/// <returns>a sequence of successful trial counts.</returns>
public IEnumerable<int> Samples()
{
while (true)
{
yield return DoSample(RandomSource, _p);
}
}
#endregion
/// <summary>
/// Samples one categorical distributed random variable; also known as the Discrete distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <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>One random integer between 0 and the size of the categorical (exclusive).</returns>
public static int Sample(System.Random rnd, double[] p)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
// The cumulative density of p.
double[] cp = UnnormalizedCDF(p);
return DoSample(rnd, cp);
}
/// <summary>
/// Samples a categorically distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <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><paramref name="n"/> random integers between 0 and the size of the categorical (exclusive).</returns>
public static IEnumerable<int> Samples(System.Random rnd, double[] p)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
// The cumulative density of p.
double[] cp = UnnormalizedCDF(p);
while (true)
{
yield return DoSample(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)
{
double[] cp = (double[]) p.Clone();
for (int 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 <see cref="cdf"/>.</returns>
internal static int DoSample(System.Random rnd, double[] cdf)
{
// TODO : use binary search to speed up this procedure.
double u = rnd.NextDouble() * cdf[cdf.Length - 1];
int idx = 0;
while (u > cdf[idx])
{
idx++;
}
return idx;
}
}
}

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

@ -1,4 +1,4 @@
// <copyright file="Bernoulli.cs" company="Math.NET">
// <copyright file="Multinomial.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -52,6 +52,11 @@ namespace MathNet.Numerics.Distributions
/// </summary>
private double[] _p;
/// <summary>
/// The number of trials.
/// </summary>
private int _n;
/// <summary>
/// The distribution's random number generator.
/// </summary>
@ -62,10 +67,12 @@ namespace MathNet.Numerics.Distributions
/// </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>
/// <exception cref="ArgumentException">If any of the probabilities are negative or do not sum to one.</exception>
public Multinomial(double[] p)
/// <param name="n">The number of trials.</param>
/// <exception cref="ArgumentOutOfRangeException">If any of the probabilities are negative or do not sum to one.</exception>
/// <exception cref="ArgumentOutOfRangeException">If <paramref name="n"/> is negative.</exception>
public Multinomial(double[] p, int n)
{
SetParameters(p);
SetParameters(p, n);
RandomSource = new System.Random();
}
@ -93,7 +100,7 @@ namespace MathNet.Numerics.Distributions
/// </summary>
public override string ToString()
{
return "Multinomial(Dimension = " + _p.Length + ")";
return "Multinomial(Dimension = " + _p.Length + ", Number of Trails = " + _n + ")";
}
/// <summary>
@ -101,8 +108,10 @@ namespace MathNet.Numerics.Distributions
/// </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>If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true</returns>
private static bool IsValidParameterSet(double[] p)
/// <param name="n">The number of trials.</param>
/// <returns>If any of the probabilities are negative returns false,
/// if the sum of parameters is 0.0, or if the number of trials is negative; otherwise true</returns>
private static bool IsValidParameterSet(double[] p, int n)
{
double sum = 0.0;
for (int i = 0; i < p.Length; i++)
@ -122,6 +131,11 @@ namespace MathNet.Numerics.Distributions
return false;
}
if (n < 0)
{
return false;
}
return true;
}
@ -130,19 +144,21 @@ namespace MathNet.Numerics.Distributions
/// </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>
/// <param name="n">The number of trials.</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
private void SetParameters(double[] p)
private void SetParameters(double[] p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
_p = (double[])p.Clone();
_n = n;
}
/// <summary>
/// Gets or sets the probability of generating a one.
/// Gets or sets the proportion of ratios.
/// </summary>
public double[] P
{
@ -153,7 +169,23 @@ namespace MathNet.Numerics.Distributions
set
{
SetParameters(value);
SetParameters(value, _n);
}
}
/// <summary>
/// Gets or sets the number of trials.
/// </summary>
public int N
{
get
{
return _n;
}
set
{
SetParameters(_p, value);
}
}
@ -179,100 +211,84 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// Samples one multinomial distributed random variable; also known as the Discrete distribution.
/// Samples one multinomial distributed random variable.
/// </summary>
/// <returns>One random integer between 0 and the size of the multinomial (exclusive).</returns>
public int Sample()
/// <returns>the counts for each of the different possible values.</returns>
public int[] Sample()
{
return Sample(RandomSource, _p);
return Sample(RandomSource, _p, _n);
}
/// <summary>
/// Samples a multinomially distributed random variable.
/// Samples a sequence multinomially distributed random variables.
/// </summary>
/// <param name="n">The number of variables needed.</param>
/// <returns><paramref name="n"/> random integers between 0 and the size of the multinomial (exclusive).</returns>
public int[] Sample(int n)
/// <returns>a sequence of counts for each of the different possible values.</returns>
public IEnumerable<int[]> Samples()
{
return Sample(RandomSource, n, _p);
while (true)
{
yield return Sample(RandomSource, _p, _n);
}
}
/// <summary>
/// Samples one multinomial distributed random variable; also known as the Discrete distribution.
/// Samples one multinomial distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <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>One random integer between 0 and the size of the multinomial (exclusive).</returns>
public static int Sample(System.Random rnd, double[] p)
/// <param name="n">The number of trials.</param>
/// <returns>the counts for each of the different possible values.</returns>
public static int[] Sample(System.Random rnd, double[] p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
// The cumulative density of p.
double[] cp = UnnormalizedCDF(p);
double[] cp = Categorical.UnnormalizedCDF(p);
// The variable that stores the counts.
int[] ret = new int[p.Length];
double u = rnd.NextDouble()*cp[cp.Length - 1];
int idx = 0;
while (u > cp[idx])
for (int i = 0; i < n; i++)
{
idx++;
ret[Categorical.DoSample(rnd, cp)]++;
}
return idx;
return ret;
}
/// <summary>
/// Samples a multinomially distributed random variable.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="n">The number of variables needed.</param>
/// <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><paramref name="n"/> random integers between 0 and the size of the multinomial (exclusive).</returns>
public static int[] Sample(System.Random rnd, int n, double[] p)
/// <param name="n">The number of variables needed.</param>
/// <returns>a sequence of counts for each of the different possible values.</returns>
public static IEnumerable<int[]> Samples(System.Random rnd, double[] p, int n)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
// The cumulative density of p.
double[] cp = UnnormalizedCDF(p);
double[] cp = Categorical.UnnormalizedCDF(p);
int[] arr = new int[n];
for (int i = 0; i < n; i++)
while (true)
{
double u = rnd.NextDouble()*cp[cp.Length - 1];
int idx = 0;
while (u > cp[idx])
// The variable that stores the counts.
int[] ret = new int[p.Length];
for (int i = 0; i < n; i++)
{
idx++;
ret[Categorical.DoSample(rnd, cp)]++;
}
arr[i] = idx;
}
return arr;
}
/// <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>
private static double[] UnnormalizedCDF(double[] p)
{
double[] cp = (double[]) p.Clone();
for (int i = 1; i < p.Length; i++)
{
cp[i] += cp[i - 1];
yield return ret;
}
return cp;
}
}
}

2
src/Numerics/Numerics.csproj

@ -79,6 +79,8 @@
<Compile Include="Distributions\Continuous\Gamma.cs" />
<Compile Include="Distributions\Continuous\Normal.cs" />
<Compile Include="Distributions\Discrete\Bernoulli.cs" />
<Compile Include="Distributions\Discrete\Binomial.cs" />
<Compile Include="Distributions\Discrete\Categorical.cs" />
<Compile Include="Distributions\Discrete\DiscreteUniform.cs" />
<Compile Include="Distributions\IContinuousDistribution.cs" />
<Compile Include="Distributions\IDiscreteDistribution.cs" />

4
src/UnitTests/DistributionTests/CommonDistributionTests.cs

@ -41,7 +41,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
[SetUp]
public void SetupDistributions()
{
dists = new IDistribution[8];
dists = new IDistribution[10];
dists[0] = new Beta(1.0, 1.0);
dists[1] = new ContinuousUniform(0.0, 1.0);
@ -51,6 +51,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
dists[5] = new Weibull(1.0, 1.0);
dists[6] = new DiscreteUniform(1, 10);
dists[7] = new LogNormal(1.0, 1.0);
dists[8] = new Binomial(0.7, 10);
dists[9] = new Categorical(0.7);
}
[Test]

248
src/UnitTests/DistributionTests/Discrete/BinomialTests.cs

@ -0,0 +1,248 @@
// <copyright file="BinomialTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.UnitTests.DistributionTests
{
using System;
using System.Linq;
using MbUnit.Framework;
using MathNet.Numerics.Distributions;
[TestFixture]
public class BinomialTests
{
[SetUp]
public void SetUp()
{
Control.CheckDistributionParameters = true;
}
[Test]
[Row(0.0, 4)]
[Row(0.3, 3)]
[Row(1.0, 2)]
public void CanCreateBinomial(double p, int n)
{
var bernoulli = new Binomial(p,n);
AssertEx.AreEqual<double>(p, bernoulli.P);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
[Row(Double.NaN, 1)]
[Row(-1.0, 1)]
[Row(2.0, 1)]
[Row(0.3, -2)]
public void BinomialCreateFailsWithBadParameters(double p, int n)
{
var bernoulli = new Binomial(p,n);
}
[Test]
public void ValidateToString()
{
var b = new Binomial(0.3, 2);
AssertEx.AreEqual<string>("Binomial(Success Probability = 0.3, Number of Trials = 2)", b.ToString());
}
[Test]
[Row(0.0, 4)]
[Row(0.3, 3)]
[Row(1.0, 2)]
public void CanSetSuccessProbability(double p, int n)
{
var b = new Binomial(0.3, n);
b.P = p;
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
[Row(Double.NaN, 1)]
[Row(-1.0, 1)]
[Row(2.0, 1)]
public void SetProbabilityOfOneFails(double p, int n)
{
var b = new Binomial(0.3, n);
b.P = p;
}
[Test]
[Row(0.0, 4)]
[Row(0.3, 3)]
[Row(1.0, 2)]
public void ValidateEntropy(double p, int n)
{
var b = new Binomial(p,n);
AssertHelpers.AlmostEqual(n * (-(1.0 - p) * Math.Log(1.0 - p) - p * Math.Log(p)), b.Entropy, 14);
}
[Test]
[Row(0.0, 4)]
[Row(0.3, 3)]
[Row(1.0, 2)]
public void ValidateSkewness(double p, int n)
{
var b = new Binomial(p,n);
AssertEx.AreEqual<double>((1.0 - 2.0 * p) / Math.Sqrt(n * p * (1.0 - p)), b.Skewness);
}
[Test]
[Row(0.0, 4, 0.0)]
[Row(0.3, 3, 1.0)]
[Row(1.0, 2, 1.0)]
public void ValidateMode(double p, int n, double m)
{
var b = new Binomial(p,n);
AssertEx.AreEqual<double>(m, b.Mode);
}
[Test]
public void ValidateMinimum()
{
var b = new Binomial(0.3, 10);
AssertEx.AreEqual<int>(0, b.Minimum);
}
[Test]
public void ValidateMaximum()
{
var b = new Binomial(0.3, 10);
AssertEx.AreEqual<int>(10, b.Maximum);
}
[Test]
[Row(0.000000, 1, 0, 1.0)]
[Row(0.000000, 1, 1, 0.0)]
[Row(0.000000, 1, 1, 0.0)]
[Row(0.000000, 3, 0, 1.0)]
[Row(0.000000, 3, 1, 0.0)]
[Row(0.000000, 3, 3, 0.0)]
[Row(0.000000, 10, 0, 1.0)]
[Row(0.000000, 10, 1, 0.0)]
[Row(0.000000, 10, 10, 0.0)]
[Row(0.300000, 1, 0, 0.69999999999999995559107901499373838305473327636719)]
[Row(0.300000, 1, 1, 0.2999999999999999888977697537484345957636833190918)]
[Row(0.300000, 1, 1, 0.2999999999999999888977697537484345957636833190918)]
[Row(0.300000, 3, 0, 0.34299999999999993471888615204079956461021032657166)]
[Row(0.300000, 3, 1, 0.44099999999999992772448109690231306411849135972008)]
[Row(0.300000, 3, 3, 0.026999999999999997002397833512077451789759292859569)]
[Row(0.300000, 10, 0, 0.02824752489999998207939855277004937778546385011091)]
[Row(0.300000, 10, 1, 0.12106082099999992639752977030555903089040470780077)]
[Row(0.300000, 10, 10, 0.0000059048999999999978147480206303047454017251032868501)]
[Row(1.000000, 1, 0, 0.0)]
[Row(1.000000, 1, 1, 1.0)]
[Row(1.000000, 1, 1, 1.0)]
[Row(1.000000, 3, 0, 0.0)]
[Row(1.000000, 3, 1, 0.0)]
[Row(1.000000, 3, 3, 1.0)]
[Row(1.000000, 10, 0, 0.0)]
[Row(1.000000, 10, 1, 0.0)]
[Row(1.000000, 10, 10, 1.0)]
public void ValidateProbability(double p, int n, int x, double d)
{
var b = new Binomial(p,n);
AssertEx.AreEqual(d, b.Probability(x));
}
[Test]
[Row(0.000000, 1, 0, 0.0)]
[Row(0.000000, 1, 1, -inf)]
[Row(0.000000, 1, 1, -inf)]
[Row(0.000000, 3, 0, 0.0)]
[Row(0.000000, 3, 1, -inf)]
[Row(0.000000, 3, 3, -inf)]
[Row(0.000000, 10, 0, 0.0)]
[Row(0.000000, 10, 1, -inf)]
[Row(0.000000, 10, 10, -inf)]
[Row(0.300000, 1, 0, -0.3566749439387324423539544041072745145718090708995)]
[Row(0.300000, 1, 1, -1.2039728043259360296301803719337238685164245381839)]
[Row(0.300000, 1, 1, -1.2039728043259360296301803719337238685164245381839)]
[Row(0.300000, 3, 0, -1.0700248318161973270618632123218235437154272126985)]
[Row(0.300000, 3, 1, -0.81871040353529122294284394322574719301255212216016)]
[Row(0.300000, 3, 3, -3.6119184129778080888905411158011716055492736145517)]
[Row(0.300000, 10, 0, -3.566749439387324423539544041072745145718090708995)]
[Row(0.300000, 10, 1, -2.1114622067804823267977785542148302920616046876506)]
[Row(0.300000, 10, 10, -12.039728043259360296301803719337238685164245381839)]
[Row(1.000000, 1, 0, -inf)]
[Row(1.000000, 1, 1, 0.0)]
[Row(1.000000, 1, 1, 0.0)]
[Row(1.000000, 3, 0, -inf)]
[Row(1.000000, 3, 1, -inf)]
[Row(1.000000, 3, 3, 0.0)]
[Row(1.000000, 10, 0, -inf)]
[Row(1.000000, 10, 1, -inf)]
[Row(1.000000, 10, 10, 0.0)]
public void ValidateProbabilityLn(double p, int n, int x, double dln)
{
var b = new Binomial(p,n);
AssertEx.AreEqual(dln, b.ProbabilityLn(x));
}
[Test]
public void CanSampleStatic()
{
var d = Binomial.Sample(new Random(), 0.3, 5);
}
[Test]
public void CanSampleSequenceStatic()
{
var ied = Binomial.Samples(new Random(), 0.3, 5);
var arr = ied.Take(5).ToArray();
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleStatic()
{
var d = Binomial.Sample(new Random(), -1.0, 5);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleSequenceStatic()
{
var ied = Binomial.Samples(new Random(), -1.0, 5).First();
}
[Test]
public void CanSample()
{
var n = new Binomial(0.3, 5);
var d = n.Sample();
}
[Test]
public void CanSampleSequence()
{
var n = new Binomial(0.3, 5);
var ied = n.Samples();
var e = ied.Take(5).ToArray();
}
}
}

117
src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs

@ -0,0 +1,117 @@
// <copyright file="CategorialTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.UnitTests.DistributionTests
{
using System;
using System.Linq;
using MbUnit.Framework;
using MathNet.Numerics.Distributions;
[TestFixture]
public class CategoricalTests
{
double[] badP;
double[] badP2;
double[] smallP;
double[] largeP;
[SetUp]
public void SetUp()
{
Control.CheckDistributionParameters = true;
badP = new double[] { -1.0, 1.0 };
badP2 = new double[] { 0.0, 0.0 };
smallP = new double[] { 1.0, 1.0, 1.0 };
largeP = new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };
}
[Test]
public void CanCreateCategorical()
{
var m = new Categorical(largeP, 4);
AssertEx.AreEqual<double[]>(largeP, m.P);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void CategoricalCreateFailsWithNegativeRatios()
{
var m = new Categorical(badP, 4);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void CategoricalCreateFailsWithAllZeroRatios()
{
var m = new Categorical(badP2, 4);
}
[Test]
public void ValidateToString()
{
var b = new Categorical(smallP);
AssertEx.AreEqual<string>("Categorical(Dimension = 3)", b.ToString());
}
[Test]
public void CanSetProbability()
{
var b = new Categorical(largeP, 4);
b.P = smallP;
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void SetProbabilityFails()
{
var b = new Categorical(largeP, 4);
b.P = badP;
}
[Test]
public void CanSampleStatic()
{
var d = Categorical.Sample(new Random(), largeP, 4);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleStatic()
{
var d = Categorical.Sample(new Random(), badP, 4);
}
[Test]
public void CanSample()
{
var n = new Categorical(largeP, 4);
var d = n.Sample();
}
}
}

18
src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs

@ -54,7 +54,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
[Test]
public void CanCreateMultinomial()
{
var m = new Multinomial(largeP);
var m = new Multinomial(largeP, 4);
AssertEx.AreEqual<double[]>(largeP, m.P);
}
@ -62,27 +62,27 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void MultinomialCreateFailsWithNegativeRatios()
{
var m = new Multinomial(badP);
var m = new Multinomial(badP, 4);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void MultinomialCreateFailsWithAllZeroRatios()
{
var m = new Multinomial(badP2);
var m = new Multinomial(badP2, 4);
}
[Test]
public void ValidateToString()
{
var b = new Multinomial(smallP);
var b = new Multinomial(smallP, 4);
AssertEx.AreEqual<string>("Multinomial(Dimension = 3)", b.ToString());
}
[Test]
public void CanSetProbability()
{
var b = new Multinomial(largeP);
var b = new Multinomial(largeP, 4);
b.P = smallP;
}
@ -90,27 +90,27 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void SetProbabilityFails()
{
var b = new Multinomial(largeP);
var b = new Multinomial(largeP, 4);
b.P = badP;
}
[Test]
public void CanSampleStatic()
{
var d = Multinomial.Sample(new Random(), largeP);
var d = Multinomial.Sample(new Random(), largeP, 4);
}
[Test]
[ExpectedException(typeof(ArgumentOutOfRangeException))]
public void FailSampleStatic()
{
var d = Multinomial.Sample(new Random(), badP);
var d = Multinomial.Sample(new Random(), badP, 4);
}
[Test]
public void CanSample()
{
var n = new Multinomial(largeP);
var n = new Multinomial(largeP, 4);
var d = n.Sample();
}
}

2
src/UnitTests/UnitTests.csproj

@ -72,6 +72,8 @@
<Compile Include="DistributionTests\Continuous\GammaTests.cs" />
<Compile Include="DistributionTests\Continuous\NormalTests.cs" />
<Compile Include="DistributionTests\Discrete\BernoulliTests.cs" />
<Compile Include="DistributionTests\Discrete\BinomialTests.cs" />
<Compile Include="DistributionTests\Discrete\CategoricalTests.cs" />
<Compile Include="DistributionTests\Discrete\DiscreteUniformTests.cs" />
<Compile Include="DistributionTests\Multivariate\DirichletTests.cs" />
<Compile Include="DistributionTests\Multivariate\MultinomialTests.cs" />

Loading…
Cancel
Save