diff --git a/src/Numerics/Distributions/Discrete/Binomial.cs b/src/Numerics/Distributions/Discrete/Binomial.cs new file mode 100644 index 00000000..5c197dbc --- /dev/null +++ b/src/Numerics/Distributions/Discrete/Binomial.cs @@ -0,0 +1,391 @@ +// +// 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. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + + /// + /// Implements the binomial distribution. For details about this distribution, see + /// Wikipedia - Binomial distribution. + /// + /// The distribution is parameterized by a probability (between 0.0 and 1.0). + /// The distribution will use the by default. + /// Users can set the random number generator by using the property. + /// 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. + public class Binomial : IDiscreteDistribution + { + /// + /// Stores the normalized binomial probability. + /// + private double _p; + + /// + /// The number of trials. + /// + private int _n; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Initializes a new instance of the Binomial class. + /// + /// The success probability of a trial. + /// The number of trials. + /// If is not in the interval [0.0,1.0]. + /// If is negative. + public Binomial(double p, int n) + { + SetParameters(p, n); + RandomSource = new System.Random(); + } + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Binomial(Success Probability = " + _p + ", Number of Trials = " + _n + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// The success probability of a trial. + /// The number of trials. + /// false is not in the interval [0.0,1.0] or is negative, true otherwise. + private static bool IsValidParameterSet(double p, int n) + { + if(p < 0.0 || p > 1.0) + { + return false; + } + + if(n < 0) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// The success probability of a trial. + /// The number of trials. + /// If is not in the interval [0.0,1.0]. + /// If is negative. + private void SetParameters(double p, int n) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _p = p; + _n = n; + } + + /// + /// Gets or sets the success probability. + /// + public double P + { + get + { + return _p; + } + + set + { + SetParameters(value, _n); + } + } + + /// + /// Gets or sets the number of trials. + /// + public int N + { + get + { + return _n; + } + + set + { + SetParameters(_p, value); + } + } + + #region IDistribution Members + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } + + /// + /// Gets the mean of the distribution. + /// + public double Mean + { + get { return _p * _n; } + } + + /// + /// Gets the standard deviation of the distribution. + /// + public double StdDev + { + get { return Math.Sqrt(_p * (1.0 - _p) * _n); } + } + + /// + /// Gets the variance of the distribution. + /// + public double Variance + { + get { return _p * (1.0 - _p) * _n; } + } + + /// + /// Gets the entropy of the distribution. + /// + 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; + } + } + + /// + /// Gets the skewness of the distribution. + /// + public double Skewness + { + get { return (1.0 - 2.0 * _p) / Math.Sqrt(_n * _p * (1.0 - _p)); } + } + + /// + /// Gets the smallest element in the domain of the distributions which can be represented by an integer. + /// + public int Minimum { get { return 0; } } + + /// + /// Gets the largest element in the domain of the distributions which can be represented by an integer. + /// + public int Maximum { get { return _n; } } + + /// + /// Computes the cumulative distribution function of the Binomial distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + 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 + + /// + /// The mode of the distribution. + /// + public int Mode + { + get { return (int) Math.Floor((_n + 1) * _p); } + } + + /// + /// The median of the distribution. + /// + public int Median + { + get { throw new NotImplementedException(); } + } + + /// + /// Computes the probability of a specific value. + /// + 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); + } + + /// + /// Computes the probability of a specific value. + /// + 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); + } + + /// + /// Samples a Binomially distributed random variable. + /// + /// The number of successful trials. + public int Sample() + { + return DoSample(RandomSource, _p, _n); + } + + /// + /// Samples an array of Bernoulli distributed random variables. + /// + /// a sequence of successful trial counts. + public IEnumerable Samples() + { + while (true) + { + yield return DoSample(RandomSource, _p, _n); + } + } + + #endregion + + /// + /// Samples a binomially distributed random variable. + /// + /// The random number generator to use. + /// The success probability of a trial; must be in the interval [0.0, 1.0]. + /// The number of trials; must be positive. + /// The number of successes in trials. + 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); + } + + /// + /// Samples a sequence of binomially distributed random variable. + /// + /// The random number generator to use. + /// The success probability of a trial; must be in the interval [0.0, 1.0]. + /// The number of trials; must be positive. + /// a sequence of successful trial counts. + public static IEnumerable 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); + } + } + + /// + /// Generates a sample from the Binomial distribution without doing parameter checking. + /// + /// The random number generator to use. + /// The success probability of a trial; must be in the interval [0.0, 1.0]. + /// The number of trials; must be positive. + /// The number of successful trials. + 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; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Distributions/Discrete/Categorical.cs b/src/Numerics/Distributions/Discrete/Categorical.cs new file mode 100644 index 00000000..d532e6cc --- /dev/null +++ b/src/Numerics/Distributions/Discrete/Categorical.cs @@ -0,0 +1,422 @@ +// +// 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. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + using MathNet.Numerics.Statistics; + + + /// + /// Implements the categorical distribution. For details about this distribution, see + /// Wikipedia - Categorical distribution. This + /// distribution is sometimes called the Discrete distribution. + /// + /// 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. + /// The distribution will use the by default. + /// Users can set the random number generator by using the property. + /// 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. + public class Categorical : IDiscreteDistribution + { + /// + /// Stores the normalized categorical probabilities. + /// + private double[] _p; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Initializes a new instance of the Categorical class. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// If any of the probabilities are negative or do not sum to one. + public Categorical(double[] p) + { + SetParameters(p); + RandomSource = new System.Random(); + } + + /* TODO + /// + /// Generate a categorical distribution from histogram . The distribution will + /// not be automatically updated when the histogram changes. + /// + 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(); + }*/ + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Categorical(Dimension = " + _p.Length + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true + 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; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// When the parameters don't pass the function. + private void SetParameters(double[] p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _p = (double[])p.Clone(); + } + + /// + /// Gets or sets the probability of generating a one. + /// + public double[] P + { + get + { + return (double[]) _p.Clone(); + } + + set + { + SetParameters(value); + } + } + + #region IDistribution Members + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } + + /// + /// Gets the mean of the distribution. + /// + public double Mean + { + get { return _p.Mean(); } + } + + /// + /// Gets the standard deviation of the distribution. + /// + public double StdDev + { + get { return _p.StandardDeviation(); } + } + + /// + /// Gets the variance of the distribution. + /// + public double Variance + { + get { return _p.Variance(); } + } + + /// + /// Gets the entropy of the distribution. + /// + 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; + } + } + + /// + /// Gets the skewness of the distribution. + /// + public double Skewness + { + get { throw new NotImplementedException(); } + } + + /// + /// Gets the smallest element in the domain of the distributions which can be represented by an integer. + /// + public int Minimum { get { return 0; } } + + /// + /// Gets the largest element in the domain of the distributions which can be represented by an integer. + /// + public int Maximum { get { return _p.Length-1; } } + + /// + /// Computes the cumulative distribution function of the Binomial distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + 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 + + /// + /// The mode of the distribution. + /// + public int Mode + { + get { throw new NotImplementedException(); } + } + + /// + /// The median of the distribution. + /// + public int Median + { + get { return (int) _p.Median(); } + } + + /// + /// Computes the probability of a specific value. + /// + public double Probability(int val) + { + if (val < 0) + { + return 0.0; + } + + if (val >= _p.Length) + { + return 0.0; + } + + return _p[val]; + } + + /// + /// Computes the probability of a specific value. + /// + public double ProbabilityLn(int val) + { + if (val < 0) + { + return 0.0; + } + + if (val >= _p.Length) + { + return 0.0; + } + + return Math.Log(_p[val]); + } + + /// + /// Samples a Binomially distributed random variable. + /// + /// The number of successful trials. + public int Sample() + { + return DoSample(RandomSource, _p); + } + + /// + /// Samples an array of Bernoulli distributed random variables. + /// + /// a sequence of successful trial counts. + public IEnumerable Samples() + { + while (true) + { + yield return DoSample(RandomSource, _p); + } + } + + #endregion + + /// + /// Samples one categorical distributed random variable; also known as the Discrete distribution. + /// + /// The random number generator to use. + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// One random integer between 0 and the size of the categorical (exclusive). + 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); + } + + /// + /// Samples a categorically distributed random variable. + /// + /// The random number generator to use. + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// random integers between 0 and the size of the categorical (exclusive). + public static IEnumerable 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); + } + } + + /// + /// Computes the unnormalized cumulative distribution function. This method performs no + /// parameter checking. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// An array representing the unnormalized cumulative distribution function. + 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; + } + + /// + /// Returns one trials from the categorical distribution. + /// + /// The random number generator to use. + /// The cumulative distribution of the probability distribution. + /// One sample from the categorical distribution implied by . + 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; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs index 0b3e00b3..cad0ceaf 100644 --- a/src/Numerics/Distributions/Multivariate/Multinomial.cs +++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -52,6 +52,11 @@ namespace MathNet.Numerics.Distributions /// private double[] _p; + /// + /// The number of trials. + /// + private int _n; + /// /// The distribution's random number generator. /// @@ -62,10 +67,12 @@ namespace MathNet.Numerics.Distributions /// /// An array of nonnegative ratios: this array does not need to be normalized /// as this is often impossible using floating point arithmetic. - /// If any of the probabilities are negative or do not sum to one. - public Multinomial(double[] p) + /// The number of trials. + /// If any of the probabilities are negative or do not sum to one. + /// If is negative. + public Multinomial(double[] p, int n) { - SetParameters(p); + SetParameters(p, n); RandomSource = new System.Random(); } @@ -93,7 +100,7 @@ namespace MathNet.Numerics.Distributions /// public override string ToString() { - return "Multinomial(Dimension = " + _p.Length + ")"; + return "Multinomial(Dimension = " + _p.Length + ", Number of Trails = " + _n + ")"; } /// @@ -101,8 +108,10 @@ namespace MathNet.Numerics.Distributions /// /// An array of nonnegative ratios: this array does not need to be normalized /// as this is often impossible using floating point arithmetic. - /// If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true - private static bool IsValidParameterSet(double[] p) + /// The number of trials. + /// 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 + 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 /// /// An array of nonnegative ratios: this array does not need to be normalized /// as this is often impossible using floating point arithmetic. + /// The number of trials. /// When the parameters don't pass the function. - 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; } /// - /// Gets or sets the probability of generating a one. + /// Gets or sets the proportion of ratios. /// public double[] P { @@ -153,7 +169,23 @@ namespace MathNet.Numerics.Distributions set { - SetParameters(value); + SetParameters(value, _n); + } + } + + /// + /// Gets or sets the number of trials. + /// + public int N + { + get + { + return _n; + } + + set + { + SetParameters(_p, value); } } @@ -179,100 +211,84 @@ namespace MathNet.Numerics.Distributions } /// - /// Samples one multinomial distributed random variable; also known as the Discrete distribution. + /// Samples one multinomial distributed random variable. /// - /// One random integer between 0 and the size of the multinomial (exclusive). - public int Sample() + /// the counts for each of the different possible values. + public int[] Sample() { - return Sample(RandomSource, _p); + return Sample(RandomSource, _p, _n); } /// - /// Samples a multinomially distributed random variable. + /// Samples a sequence multinomially distributed random variables. /// - /// The number of variables needed. - /// random integers between 0 and the size of the multinomial (exclusive). - public int[] Sample(int n) + /// a sequence of counts for each of the different possible values. + public IEnumerable Samples() { - return Sample(RandomSource, n, _p); + while (true) + { + yield return Sample(RandomSource, _p, _n); + } } /// - /// Samples one multinomial distributed random variable; also known as the Discrete distribution. + /// Samples one multinomial distributed random variable. /// /// The random number generator to use. /// An array of nonnegative ratios: this array does not need to be normalized /// as this is often impossible using floating point arithmetic. - /// One random integer between 0 and the size of the multinomial (exclusive). - public static int Sample(System.Random rnd, double[] p) + /// The number of trials. + /// the counts for each of the different possible values. + 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; } /// /// Samples a multinomially distributed random variable. /// /// The random number generator to use. - /// The number of variables needed. /// An array of nonnegative ratios: this array does not need to be normalized /// as this is often impossible using floating point arithmetic. - /// random integers between 0 and the size of the multinomial (exclusive). - public static int[] Sample(System.Random rnd, int n, double[] p) + /// The number of variables needed. + /// a sequence of counts for each of the different possible values. + public static IEnumerable 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; - } - - /// - /// Computes the unnormalized cumulative distribution function. This method performs no - /// parameter checking. - /// - /// An array of nonnegative ratios: this array does not need to be normalized - /// as this is often impossible using floating point arithmetic. - /// An array representing the unnormalized cumulative distribution function. - 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; } } } \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 70730d90..25f1d94f 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -79,6 +79,8 @@ + + diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index b6018479..b7456df1 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/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] diff --git a/src/UnitTests/DistributionTests/Discrete/BinomialTests.cs b/src/UnitTests/DistributionTests/Discrete/BinomialTests.cs new file mode 100644 index 00000000..9ffc4960 --- /dev/null +++ b/src/UnitTests/DistributionTests/Discrete/BinomialTests.cs @@ -0,0 +1,248 @@ +// +// 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. +// + +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(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("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((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(m, b.Mode); + } + + [Test] + public void ValidateMinimum() + { + var b = new Binomial(0.3, 10); + AssertEx.AreEqual(0, b.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var b = new Binomial(0.3, 10); + AssertEx.AreEqual(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(); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs b/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs new file mode 100644 index 00000000..d748f792 --- /dev/null +++ b/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs @@ -0,0 +1,117 @@ +// +// 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. +// + +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(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("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(); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs index 7bfd0fdb..a0f46115 100644 --- a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs +++ b/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(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("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(); } } diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 45b72e39..d669620c 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -72,6 +72,8 @@ + +