From b2be48f200007b7979df5157e62ace91ad3ee152 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 23 Jun 2013 15:11:49 +0200 Subject: [PATCH] Categorical: consistent pmf vs cdf parameters, do not expect them to be normalized --- .../Distributions/Discrete/Categorical.cs | 257 +++++++++++------- .../Distributions/Multivariate/Multinomial.cs | 4 +- .../Discrete/CategoricalTests.cs | 4 +- 3 files changed, 167 insertions(+), 98 deletions(-) diff --git a/src/Numerics/Distributions/Discrete/Categorical.cs b/src/Numerics/Distributions/Discrete/Categorical.cs index f27ff56c..f69f719e 100644 --- a/src/Numerics/Distributions/Discrete/Categorical.cs +++ b/src/Numerics/Distributions/Discrete/Categorical.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// +// Copyright (c) 2009-2013 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 @@ -12,8 +14,10 @@ // 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 @@ -47,25 +51,19 @@ namespace MathNet.Numerics.Distributions /// to false, all parameter checks can be turned off. public class Categorical : IDiscreteDistribution { - /// - /// Stores the unnormalized categorical probabilities. - /// - double[] _p; - - /// - /// The distribution's random number generator. - /// Random _random; + double[] _pmfNormalized; + double[] _cdfUnnormalized; /// /// Initializes a new instance of the Categorical class. /// - /// An array of nonnegative ratios: this array does not need to be normalized + /// 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) + public Categorical(double[] probabilityMass) { - SetParameters(p); + SetParameters(probabilityMass); RandomSource = new Random(); } @@ -101,20 +99,20 @@ namespace MathNet.Numerics.Distributions /// a string representation of the distribution. public override string ToString() { - return "Categorical(Dimension = " + _p.Length + ")"; + return "Categorical(Dimension = " + _pmfNormalized.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. + /// 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 - static bool IsValidParameterSet(IEnumerable p) + static bool IsValidProbabilityMass(double[] p) { var sum = 0.0; - foreach (double t in p) + for (int i = 0; i < p.Length; i++) { + double t = p[i]; if (t < 0.0 || Double.IsNaN(t)) { return false; @@ -123,7 +121,29 @@ namespace MathNet.Numerics.Distributions sum += t; } - return sum != 0.0; + return sum > 0.0; + } + + /// + /// 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 + static bool IsValidCumulativeDistribution(double[] cdf) + { + var last = 0.0; + for (int i = 0; i < cdf.Length; i++) + { + double t = cdf[i]; + if (t < 0.0 || Double.IsNaN(t) || t < last) + { + return false; + } + + last = t; + } + + return last > 0.0; } /// @@ -131,15 +151,30 @@ 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. - /// When the parameters don't pass the function. + /// When the parameters don't pass the function. void SetParameters(double[] p) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + if (Control.CheckDistributionParameters && !IsValidProbabilityMass(p)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - _p = (double[])p.Clone(); + // Extract unnormalized cumulative distribution + _cdfUnnormalized = new double[p.Length]; + _cdfUnnormalized[0] = p[0]; + for (int i = 1; i < p.Length; i++) + { + _cdfUnnormalized[i] = _cdfUnnormalized[i - 1] + p[i]; + } + + // Extract normalized probability mass + var sum = _cdfUnnormalized[_cdfUnnormalized.Length - 1]; + _pmfNormalized = new double[p.Length]; + for (int i = 0; i < p.Length; i++) + { + _pmfNormalized[i] = p[i]/sum; + } + } /// @@ -149,20 +184,7 @@ namespace MathNet.Numerics.Distributions /// exactly in a floating point representation. public double[] P { - get - { - var p = (double[])_p.Clone(); - - var sum = p.Sum(); - - for (var i = 0; i < p.Length; i++) - { - p[i] /= sum; - } - - return p; - } - + get { return (double[])_pmfNormalized.Clone(); } set { SetParameters(value); } } @@ -191,7 +213,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get { return _p.Mean(); } + get { return _pmfNormalized.Mean(); } } /// @@ -199,7 +221,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get { return _p.StandardDeviation(); } + get { return _pmfNormalized.StandardDeviation(); } } /// @@ -207,7 +229,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get { return _p.Variance(); } + get { return _pmfNormalized.Variance(); } } /// @@ -215,7 +237,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get { return _p.Sum(p => p * Math.Log(p)); } + get { return _pmfNormalized.Sum(p => p * Math.Log(p)); } } /// @@ -240,7 +262,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get { return _p.Length - 1; } + get { return _pmfNormalized.Length - 1; } } /// @@ -255,13 +277,12 @@ namespace MathNet.Numerics.Distributions return 0.0; } - if (x >= _p.Length) + if (x >= _cdfUnnormalized.Length) { return 1.0; } - var cdf = UnnormalizedCdf(_p); - return cdf[(int)Math.Floor(x)] / cdf[_p.Length - 1]; + return _cdfUnnormalized[(int) Math.Floor(x)]/_cdfUnnormalized[_cdfUnnormalized.Length - 1]; } #endregion @@ -282,7 +303,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get { return (int)_p.Median(); } + get { return (int)_pmfNormalized.Median(); } } /// @@ -297,12 +318,12 @@ namespace MathNet.Numerics.Distributions return 0.0; } - if (k >= _p.Length) + if (k >= _pmfNormalized.Length) { return 0.0; } - return _p[k]; + return _pmfNormalized[k]; } /// @@ -317,48 +338,48 @@ namespace MathNet.Numerics.Distributions return 0.0; } - if (k >= _p.Length) + if (k >= _pmfNormalized.Length) { return 0.0; } - return Math.Log(_p[k]); + return Math.Log(_pmfNormalized[k]); } #endregion /// - /// Computes the unnormalized cumulative distribution function. This method performs no - /// parameter checking. + /// Computes the cumulative distribution function. This method performs no parameter checking. + /// If the probability mass was normalized, the resulting cumulative distribution is normalized as well (up to numerical errors). /// - /// An array of nonnegative ratios: this array does not need to be normalized + /// 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) + internal static double[] ProbabilityMassToCumulativeDistribution(double[] pmfUnnormalized) { - var cp = (double[])p.Clone(); - - for (var i = 1; i < p.Length; i++) + var cdfUnnormalized = new double[pmfUnnormalized.Length]; + cdfUnnormalized[0] = pmfUnnormalized[0]; + for (int i = 1; i < pmfUnnormalized.Length; i++) { - cp[i] += cp[i - 1]; + cdfUnnormalized[i] = cdfUnnormalized[i - 1] + pmfUnnormalized[i]; } - return cp; + return cdfUnnormalized; } /// /// 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 SampleUnchecked(Random rnd, double[] cdf) + /// The (unnormalized) cumulative distribution of the probability distribution. + /// One sample from the categorical distribution implied by . + internal static int SampleUnchecked(Random rnd, double[] cdfUnnormalized) { // TODO : use binary search to speed up this procedure. - var u = rnd.NextDouble() * cdf[cdf.Length - 1]; + var u = rnd.NextDouble() * cdfUnnormalized[cdfUnnormalized.Length - 1]; var idx = 0; - while (u > cdf[idx]) + while (u > cdfUnnormalized[idx]) { idx++; } @@ -372,7 +393,7 @@ namespace MathNet.Numerics.Distributions /// The number of successful trials. public int Sample() { - return Sample(RandomSource, _p); + return SampleUnchecked(RandomSource, _cdfUnnormalized); } /// @@ -381,27 +402,56 @@ namespace MathNet.Numerics.Distributions /// a sequence of successful trial counts. public IEnumerable Samples() { - return Samples(RandomSource, _p); + while (true) + { + yield return SampleUnchecked(RandomSource, _cdfUnnormalized); + } } /// /// 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 + /// 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(Random rnd, double[] p) + [Obsolete("Use SampleWithProbabilityMass instead (or SampleWithCumulativeDistribution which is faster). Scheduled for removal in v3.0.")] + public static int Sample(Random rnd, double[] pmfUnnormalized) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + return SampleWithProbabilityMass(rnd, pmfUnnormalized); + } + + /// + /// Samples one categorical distributed random variable; also known as the Discrete distribution. + /// + /// The random number generator to use. + /// An array of the cumulative distribution. Not assumed to be normalized. + /// One random integer between 0 and the size of the categorical (exclusive). + public static int SampleWithCumulativeDistribution(Random rnd, double[] cdfUnnormalized) + { + if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - // The cumulative density of p. - var cp = UnnormalizedCdf(p); + return SampleUnchecked(rnd, cdfUnnormalized); + } - return SampleUnchecked(rnd, cp); + /// + /// Samples one categorical distributed random variable; also known as the Discrete distribution. + /// + /// The random number generator to use. + /// An array of nonnegative ratios. Not assumed to be normalized. + /// One random integer between 0 and the size of the categorical (exclusive). + public static int SampleWithProbabilityMass(Random rnd, double[] pmfUnnormalized) + { + if (Control.CheckDistributionParameters && !IsValidProbabilityMass(pmfUnnormalized)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + var cdf = ProbabilityMassToCumulativeDistribution(pmfUnnormalized); + return SampleUnchecked(rnd, cdf); } /// @@ -411,19 +461,48 @@ 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. /// random integers between 0 and the size of the categorical (exclusive). + [Obsolete("Use SamplesWithProbabilityMass instead (or SamplesWithCumulativeDistribution which is faster). Scheduled for removal in v3.0.")] public static IEnumerable Samples(Random rnd, double[] p) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + return SamplesWithProbabilityMass(rnd, p); + } + + /// + /// Samples a categorically distributed random variable. + /// + /// The random number generator to use. + /// An array of the cumulative distribution. Not assumed to be normalized. + /// random integers between 0 and the size of the categorical (exclusive). + public static IEnumerable SamplesWithCumulativeDistribution(Random rnd, double[] cdfUnnormalized) + { + if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - // The cumulative density of p. - var cp = UnnormalizedCdf(p); + while (true) + { + yield return SampleUnchecked(rnd, cdfUnnormalized); + } + } + /// + /// Samples a categorically distributed random variable. + /// + /// The random number generator to use. + /// An array of nonnegative ratios. Not assumed to be normalized. + /// random integers between 0 and the size of the categorical (exclusive). + public static IEnumerable SamplesWithProbabilityMass(Random rnd, double[] pmfUnnormalized) + { + if (Control.CheckDistributionParameters && !IsValidProbabilityMass(pmfUnnormalized)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + var cdf = ProbabilityMassToCumulativeDistribution(pmfUnnormalized); while (true) { - yield return SampleUnchecked(rnd, cp); + yield return SampleUnchecked(rnd, cdf); } } @@ -431,34 +510,24 @@ namespace MathNet.Numerics.Distributions /// Returns the inverse of the distribution function for the categorical distribution /// specified by the given normalized CDF, for the given probability. /// - /// An array corresponding to a normalized CDF for a categorical distribution. + /// An array corresponding to a CDF for a categorical distribution. Not assumed to be normalized. /// A real number between 0 and 1. /// An integer between 0 and the size of the categorical (exclusive), /// that corresponds to the inverse CDF for the given probability. - public static int InverseCumulativeDistribution(double[] normalizedCDF, double probability) + public static int InverseCumulativeDistribution(double[] cdfUnnormalized, double probability) { - if (Control.CheckDistributionParameters) + if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized)) { - if (probability < 0.0 || probability > 1.0 || Double.IsNaN(probability)) - { - throw new ArgumentOutOfRangeException("probability"); - } - - if (normalizedCDF[0] < 0.0 || normalizedCDF[0] > 1.0 || Double.IsNaN(normalizedCDF[0])) - throw new ArgumentOutOfRangeException("normalizedCDF"); - - for (var i = 1; i < normalizedCDF.Length; i++) - { - var cd = normalizedCDF[i]; - if (cd < 0.0 || cd > 1.0 || cd < normalizedCDF[i - 1] || Double.IsNaN(cd)) - { - throw new ArgumentOutOfRangeException("normalizedCDF"); - } - } + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - int idx = Array.BinarySearch(normalizedCDF, probability); + if (probability < 0.0 || probability > 1.0 || Double.IsNaN(probability)) + { + throw new ArgumentOutOfRangeException("probability"); + } + var denormalizedProbability = probability*cdfUnnormalized[cdfUnnormalized.Length - 1]; + int idx = Array.BinarySearch(cdfUnnormalized, denormalizedProbability); if (idx < 0) { idx = ~idx; diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs index 6816902d..43d89f24 100644 --- a/src/Numerics/Distributions/Multivariate/Multinomial.cs +++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs @@ -361,7 +361,7 @@ namespace MathNet.Numerics.Distributions } // The cumulative density of p. - var cp = Categorical.UnnormalizedCdf(p); + var cp = Categorical.ProbabilityMassToCumulativeDistribution(p); // The variable that stores the counts. var ret = new int[p.Length]; @@ -390,7 +390,7 @@ namespace MathNet.Numerics.Distributions } // The cumulative density of p. - var cp = Categorical.UnnormalizedCdf(p); + var cp = Categorical.ProbabilityMassToCumulativeDistribution(p); while (true) { diff --git a/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs b/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs index c0d98c6e..8c2037f0 100644 --- a/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs +++ b/src/UnitTests/DistributionTests/Discrete/CategoricalTests.cs @@ -161,7 +161,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Discrete [Test] public void CanSampleStatic() { - Categorical.Sample(new Random(), _largeP); + Categorical.SampleWithProbabilityMass(new Random(), _largeP); } /// @@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Discrete [Test] public void FailSampleStatic() { - Assert.Throws(() => Categorical.Sample(new Random(), _badP)); + Assert.Throws(() => Categorical.SampleWithProbabilityMass(new Random(), _badP)); } ///