diff --git a/src/Numerics/Distributions/Discrete/Bernoulli.cs b/src/Numerics/Distributions/Discrete/Bernoulli.cs index 15f86424..941c9f80 100644 --- a/src/Numerics/Distributions/Discrete/Bernoulli.cs +++ b/src/Numerics/Distributions/Discrete/Bernoulli.cs @@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions /// /// The probability of generating a one. /// - private double _p; + double _p; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Bernoulli class. @@ -77,7 +77,7 @@ namespace MathNet.Numerics.Distributions /// /// The probability of generating a one. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double p) + static bool IsValidParameterSet(double p) { if (p >= 0.0 && p <= 1.0) { @@ -92,7 +92,7 @@ namespace MathNet.Numerics.Distributions /// /// The probability of generating a one. /// When the parameters don't pass the function. - private void SetParameters(double p) + void SetParameters(double p) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) { @@ -107,15 +107,9 @@ namespace MathNet.Numerics.Distributions /// public double P { - get - { - return _p; - } + get { return _p; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } #region IDistribution Members @@ -125,10 +119,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -146,10 +137,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _p; - } + get { return _p; } } /// @@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(_p * (1.0 - _p)); - } + get { return Math.Sqrt(_p * (1.0 - _p)); } } /// @@ -168,10 +153,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _p * (1.0 - _p); - } + get { return _p * (1.0 - _p); } } /// @@ -179,10 +161,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p)); - } + get { return -(_p * Math.Log(_p)) - ((1.0 - _p) * Math.Log(1.0 - _p)); } } /// @@ -190,10 +169,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p)); - } + get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_p * (1.0 - _p)); } } /// @@ -201,10 +177,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -212,10 +185,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return 1; - } + get { return 1; } } /// @@ -229,7 +199,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x < 1.0) { return 1.0 - _p; @@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return _p > 0.5 ? 1 : 0; - } + get { return _p > 0.5 ? 1 : 0; } } /// @@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - throw new NotSupportedException("The median of the Bernoulli distribution is undefined."); - } + get { throw new NotSupportedException("The median of the Bernoulli distribution is undefined."); } } /// @@ -299,13 +263,31 @@ namespace MathNet.Numerics.Distributions return k == 1 ? Math.Log(_p) : Double.NegativeInfinity; } + #endregion + + /// + /// Generates one sample from the Bernoulli distribution. + /// + /// The random source to use. + /// The probability of generating a one. + /// A random sample from the Bernoulli distribution. + internal static int SampleUnchecked(Random rnd, double p) + { + if (rnd.NextDouble() < p) + { + return 1; + } + + return 0; + } + /// /// Samples a Bernoulli distributed random variable. /// /// A sample from the Bernoulli distribution. public int Sample() { - return DoSample(RandomSource, _p); + return SampleUnchecked(RandomSource, _p); } /// @@ -316,12 +298,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _p); + yield return SampleUnchecked(RandomSource, _p); } } - #endregion - /// /// Samples a Bernoulli distributed random variable. /// @@ -335,7 +315,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return DoSample(rnd, p); + return SampleUnchecked(rnd, p); } /// @@ -353,24 +333,9 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return DoSample(rnd, p); + yield return SampleUnchecked(rnd, p); } } - /// - /// Generates one sample from the Bernoulli distribution. - /// - /// The random source to use. - /// The probability of generating a one. - /// A random sample from the Bernoulli distribution. - private static int DoSample(Random rnd, double p) - { - if (rnd.NextDouble() < p) - { - return 1; - } - - return 0; - } } } diff --git a/src/Numerics/Distributions/Discrete/Binomial.cs b/src/Numerics/Distributions/Discrete/Binomial.cs index 5f836591..f37a5881 100644 --- a/src/Numerics/Distributions/Discrete/Binomial.cs +++ b/src/Numerics/Distributions/Discrete/Binomial.cs @@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions /// /// Stores the normalized binomial probability. /// - private double _p; + double _p; /// /// The number of trials. /// - private int _n; + int _n; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Binomial class. @@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions /// 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) + static bool IsValidParameterSet(double p, int n) { if (p < 0.0 || p > 1.0 || Double.IsNaN(p)) { @@ -107,7 +107,7 @@ namespace MathNet.Numerics.Distributions /// The number of trials. /// If is not in the interval [0.0,1.0]. /// If is negative. - private void SetParameters(double p, int n) + void SetParameters(double p, int n) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n)) { @@ -123,15 +123,9 @@ namespace MathNet.Numerics.Distributions /// public double P { - get - { - return _p; - } + get { return _p; } - set - { - SetParameters(value, _n); - } + set { SetParameters(value, _n); } } /// @@ -139,15 +133,9 @@ namespace MathNet.Numerics.Distributions /// public int N { - get - { - return _n; - } + get { return _n; } - set - { - SetParameters(_p, value); - } + set { SetParameters(_p, value); } } #region IDistribution Members @@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _p * _n; - } + get { return _p * _n; } } /// @@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(_p * (1.0 - _p) * _n); - } + get { return Math.Sqrt(_p * (1.0 - _p) * _n); } } /// @@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _p * (1.0 - _p) * _n; - } + get { return _p * (1.0 - _p) * _n; } } /// @@ -234,10 +210,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p)); - } + get { return (1.0 - (2.0 * _p)) / Math.Sqrt(_n * _p * (1.0 - _p)); } } /// @@ -245,10 +218,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -256,10 +226,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return _n; - } + get { return _n; } } /// @@ -273,7 +240,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x > _n) { return 1.0; @@ -303,7 +270,7 @@ namespace MathNet.Numerics.Distributions { return _n; } - + if (_p == 0.0) { return 0; @@ -318,10 +285,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - return (int)Math.Floor(_p * _n); - } + get { return (int)Math.Floor(_p * _n); } } /// @@ -345,7 +309,7 @@ namespace MathNet.Numerics.Distributions { return 1.0; } - + if (_p == 0.0) { return 0.0; @@ -355,7 +319,7 @@ namespace MathNet.Numerics.Distributions { return 1.0; } - + if (_p == 1.0) { return 0.0; @@ -385,7 +349,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_p == 0.0) { return Double.NegativeInfinity; @@ -395,7 +359,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_p == 1.0) { return Double.NegativeInfinity; @@ -404,13 +368,33 @@ namespace MathNet.Numerics.Distributions return SpecialFunctions.BinomialLn(_n, k) + (k * Math.Log(_p)) + ((_n - k) * Math.Log(1.0 - _p)); } + #endregion + + /// + /// 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. + internal static int SampleUnchecked(Random rnd, double p, int n) + { + var k = 0; + for (var i = 0; i < n; i++) + { + k += rnd.NextDouble() < p ? 1 : 0; + } + + return k; + } + /// /// Samples a Binomially distributed random variable. /// /// The number of successes in N trials. public int Sample() { - return DoSample(RandomSource, _p, _n); + return SampleUnchecked(RandomSource, _p, _n); } /// @@ -421,12 +405,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _p, _n); + yield return SampleUnchecked(RandomSource, _p, _n); } } - #endregion - /// /// Samples a binomially distributed random variable. /// @@ -441,7 +423,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return DoSample(rnd, p, n); + return SampleUnchecked(rnd, p, n); } /// @@ -460,26 +442,8 @@ namespace MathNet.Numerics.Distributions 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(Random rnd, double p, int n) - { - var k = 0; - for (var i = 0; i < n; i++) - { - k += rnd.NextDouble() < p ? 1 : 0; + yield return SampleUnchecked(rnd, p, n); } - - return k; } } } diff --git a/src/Numerics/Distributions/Discrete/Categorical.cs b/src/Numerics/Distributions/Discrete/Categorical.cs index 7f8674e1..97e51e86 100644 --- a/src/Numerics/Distributions/Discrete/Categorical.cs +++ b/src/Numerics/Distributions/Discrete/Categorical.cs @@ -50,12 +50,12 @@ namespace MathNet.Numerics.Distributions /// /// Stores the unnormalized categorical probabilities. /// - private double[] _p; + double[] _p; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Categorical class. @@ -110,7 +110,7 @@ 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(IEnumerable p) + static bool IsValidParameterSet(IEnumerable p) { var sum = 0.0; foreach (double t in p) @@ -119,7 +119,7 @@ namespace MathNet.Numerics.Distributions { return false; } - + sum += t; } @@ -132,7 +132,7 @@ 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. - private void SetParameters(double[] p) + void SetParameters(double[] p) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) { @@ -163,10 +163,7 @@ namespace MathNet.Numerics.Distributions return p; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } #region IDistribution Members @@ -176,10 +173,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -197,10 +191,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _p.Mean(); - } + get { return _p.Mean(); } } /// @@ -208,10 +199,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return _p.StandardDeviation(); - } + get { return _p.StandardDeviation(); } } /// @@ -219,10 +207,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _p.Variance(); - } + get { return _p.Variance(); } } /// @@ -230,10 +215,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return _p.Sum(p => p * Math.Log(p)); - } + get { return _p.Sum(p => p * Math.Log(p)); } } /// @@ -242,10 +224,7 @@ namespace MathNet.Numerics.Distributions /// Throws a . public double Skewness { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -253,10 +232,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -264,10 +240,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return _p.Length - 1; - } + get { return _p.Length - 1; } } /// @@ -281,7 +254,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x >= _p.Length) { return 1.0; @@ -301,10 +274,7 @@ namespace MathNet.Numerics.Distributions /// Throws a . public int Mode { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -312,10 +282,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - return (int)_p.Median(); - } + get { return (int)_p.Median(); } } /// @@ -358,6 +325,47 @@ namespace MathNet.Numerics.Distributions return Math.Log(_p[k]); } + #endregion + + /// + /// 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) + { + var cp = (double[])p.Clone(); + + for (var 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 SampleUnchecked(Random rnd, double[] cdf) + { + // TODO : use binary search to speed up this procedure. + var u = rnd.NextDouble() * cdf[cdf.Length - 1]; + + var idx = 0; + while (u > cdf[idx]) + { + idx++; + } + + return idx; + } + /// /// Samples a Binomially distributed random variable. /// @@ -376,8 +384,6 @@ namespace MathNet.Numerics.Distributions return Samples(RandomSource, _p); } - #endregion - /// /// Samples one categorical distributed random variable; also known as the Discrete distribution. /// @@ -395,7 +401,7 @@ namespace MathNet.Numerics.Distributions // The cumulative density of p. var cp = UnnormalizedCdf(p); - return DoSample(rnd, cp); + return SampleUnchecked(rnd, cp); } /// @@ -417,47 +423,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return DoSample(rnd, cp); + yield return SampleUnchecked(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) - { - var cp = (double[])p.Clone(); - - for (var 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(Random rnd, double[] cdf) - { - // TODO : use binary search to speed up this procedure. - var u = rnd.NextDouble() * cdf[cdf.Length - 1]; - - var idx = 0; - while (u > cdf[idx]) - { - idx++; - } - - return idx; - } } } diff --git a/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs b/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs index 719d8455..9dcdf47f 100644 --- a/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs +++ b/src/Numerics/Distributions/Discrete/ConwayMaxwellPoisson.cs @@ -52,37 +52,37 @@ namespace MathNet.Numerics.Distributions /// Since many properties of the distribution can only be computed approximately, the tolerance /// level specifies how much error we accept. /// - private const double Tolerance = 1e-12; + const double Tolerance = 1e-12; /// /// The mean of the distribution. /// - private double _mean = double.MinValue; + double _mean = double.MinValue; /// /// The variance of the distribution. /// - private double _variance = double.MinValue; + double _variance = double.MinValue; /// /// Caches the value of the normalization constant. /// - private double _z = double.MinValue; + double _z = double.MinValue; /// /// The lambda parameter. /// - private double _lambda; + double _lambda; /// /// The nu parameter. /// - private double _nu; + double _nu; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -105,7 +105,7 @@ namespace MathNet.Numerics.Distributions /// The lambda parameter. /// The nu parameter. /// When the parameters don't pass the function. - private void SetParameters(double lambda, double nu) + void SetParameters(double lambda, double nu) { if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu)) { @@ -122,7 +122,7 @@ namespace MathNet.Numerics.Distributions /// The lambda parameter. /// The nu parameter. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double lambda, double nu) + static bool IsValidParameterSet(double lambda, double nu) { if (lambda <= 0.0) { @@ -143,15 +143,9 @@ namespace MathNet.Numerics.Distributions /// The value of the lambda parameter. public double Lambda { - get - { - return _lambda; - } + get { return _lambda; } - set - { - SetParameters(value, _nu); - } + set { SetParameters(value, _nu); } } /// @@ -160,15 +154,9 @@ namespace MathNet.Numerics.Distributions /// The value of the Nu parameter. public double Nu { - get - { - return _nu; - } + get { return _nu; } - set - { - SetParameters(_lambda, value); - } + set { SetParameters(_lambda, value); } } /// @@ -189,10 +177,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -333,10 +318,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(Variance); - } + get { return Math.Sqrt(Variance); } } /// @@ -344,10 +326,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -355,10 +334,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -386,10 +362,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -397,10 +370,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -408,10 +378,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -419,10 +386,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -449,35 +413,12 @@ namespace MathNet.Numerics.Distributions return Math.Log(Probability(k)); } - /// - /// Samples a Conway-Maxwell-Poisson distributed random variable. - /// - /// a sample from the distribution. - public int Sample() - { - return DoSample(RandomSource, _lambda, _nu, Z); - } - - /// - /// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables. - /// - /// - /// a sequence of samples from a Conway-Maxwell-Poisson distribution. - /// - public IEnumerable Samples() - { - while (true) - { - yield return DoSample(RandomSource, _lambda, _nu, Z); - } - } - #endregion /// /// Gets the normalization constant of the Conway-Maxwell-Poisson distribution. /// - private double Z + double Z { get { @@ -488,7 +429,7 @@ namespace MathNet.Numerics.Distributions _z = Normalization(_lambda, _nu); return _z; - } + } } /// @@ -499,7 +440,7 @@ namespace MathNet.Numerics.Distributions /// /// an approximate normalization constant for the CMP distribution. /// - private static double Normalization(double lambda, double nu) + static double Normalization(double lambda, double nu) { // Initialize Z with the first two terms. var z = 1.0 + lambda; @@ -515,7 +456,7 @@ namespace MathNet.Numerics.Distributions // The new term. t = t * e; - + // The updated normalization constant. z = z + t; @@ -542,7 +483,7 @@ namespace MathNet.Numerics.Distributions /// /// One sample from the distribution implied by , , and . /// - private static int DoSample(Random rnd, double lambda, double nu, double z) + internal static int SampleUnchecked(Random rnd, double lambda, double nu, double z) { var u = rnd.NextDouble(); var p = 1.0 / z; @@ -558,5 +499,65 @@ namespace MathNet.Numerics.Distributions return i; } + + /// + /// Samples a Conway-Maxwell-Poisson distributed random variable. + /// + /// a sample from the distribution. + public int Sample() + { + return SampleUnchecked(RandomSource, _lambda, _nu, Z); + } + + /// + /// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables. + /// + /// + /// a sequence of samples from a Conway-Maxwell-Poisson distribution. + /// + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(RandomSource, _lambda, _nu, Z); + } + } + + /// + /// Samples a random variable. + /// + /// The random number generator to use. + /// The lambda parameter + /// The nu parameter. + public static int Sample(Random rnd, double lambda, double nu) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + var z = Normalization(lambda, nu); + return SampleUnchecked(rnd, lambda, nu, z); + } + + /// + /// Samples a sequence of this random variable. + /// + /// The random number generator to use. + /// The lambda parameter + /// The nu parameter. + public static IEnumerable Samples(Random rnd, double lambda, double nu) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda, nu)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + var z = Normalization(lambda, nu); + while (true) + { + yield return SampleUnchecked(rnd, lambda, nu, z); + } + } } } diff --git a/src/Numerics/Distributions/Discrete/DiscreteUniform.cs b/src/Numerics/Distributions/Discrete/DiscreteUniform.cs index e096f2d9..c4d20347 100644 --- a/src/Numerics/Distributions/Discrete/DiscreteUniform.cs +++ b/src/Numerics/Distributions/Discrete/DiscreteUniform.cs @@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions /// /// The distribution's lower bound. /// - private int _lower; + int _lower; /// /// The distribution's upper bound. /// - private int _upper; + int _upper; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the DiscreteUniform class. @@ -85,7 +85,7 @@ namespace MathNet.Numerics.Distributions /// Lower bound. /// Upper bound; must be at least as large as . /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(int lower, int upper) + static bool IsValidParameterSet(int lower, int upper) { if (lower <= upper) { @@ -101,7 +101,7 @@ namespace MathNet.Numerics.Distributions /// Lower bound. /// Upper bound; must be at least as large as . /// When the parameters don't pass the function. - private void SetParameters(int lower, int upper) + void SetParameters(int lower, int upper) { if (Control.CheckDistributionParameters && !IsValidParameterSet(lower, upper)) { @@ -117,15 +117,9 @@ namespace MathNet.Numerics.Distributions /// public int LowerBound { - get - { - return _lower; - } + get { return _lower; } - set - { - SetParameters(value, _upper); - } + set { SetParameters(value, _upper); } } /// @@ -133,15 +127,9 @@ namespace MathNet.Numerics.Distributions /// public int UpperBound { - get - { - return _upper; - } + get { return _upper; } - set - { - SetParameters(_lower, value); - } + set { SetParameters(_lower, value); } } #region IDistribution Members @@ -151,10 +139,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -172,10 +157,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return (_lower + _upper) / 2.0; - } + get { return (_lower + _upper) / 2.0; } } /// @@ -183,10 +165,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0); - } + get { return Math.Sqrt((((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0); } } /// @@ -194,10 +173,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0; - } + get { return (((_upper - _lower + 1.0) * (_upper - _lower + 1.0)) - 1.0) / 12.0; } } /// @@ -205,10 +181,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return Math.Log(_upper - _lower + 1.0); - } + get { return Math.Log(_upper - _lower + 1.0); } } /// @@ -216,10 +189,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -227,10 +197,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return _lower; - } + get { return _lower; } } /// @@ -238,10 +205,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return _upper; - } + get { return _upper; } } /// @@ -255,7 +219,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x >= _upper) { return 1.0; @@ -273,10 +237,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return (int)Math.Floor((_lower + _upper) / 2.0); - } + get { return (int)Math.Floor((_lower + _upper) / 2.0); } } /// @@ -284,10 +245,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - return (int)Math.Floor((_lower + _upper) / 2.0); - } + get { return (int)Math.Floor((_lower + _upper) / 2.0); } } /// @@ -324,13 +282,27 @@ namespace MathNet.Numerics.Distributions return Double.NegativeInfinity; } + #endregion + + /// + /// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking. + /// + /// The random source to use. + /// The lower bound of the uniform random variable. + /// The upper bound of the uniform random variable. + /// A random sample from the discrete uniform distribution. + internal static int SampleUnchecked(Random rnd, int lower, int upper) + { + return (rnd.Next() % (upper - lower + 1)) + lower; + } + /// /// Draws a random sample from the distribution. /// /// a sample from the distribution. public int Sample() { - return DoSample(RandomSource, _lower, _upper); + return SampleUnchecked(RandomSource, _lower, _upper); } /// @@ -341,12 +313,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _lower, _upper); + yield return SampleUnchecked(RandomSource, _lower, _upper); } } - #endregion - /// /// Samples a uniformly distributed random variable. /// @@ -361,7 +331,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return DoSample(rnd, lower, upper); + return SampleUnchecked(rnd, lower, upper); } /// @@ -380,20 +350,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return DoSample(rnd, lower, upper); + yield return SampleUnchecked(rnd, lower, upper); } } - - /// - /// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking. - /// - /// The random source to use. - /// The lower bound of the uniform random variable. - /// The upper bound of the uniform random variable. - /// A random sample from the discrete uniform distribution. - private static int DoSample(Random rnd, int lower, int upper) - { - return (rnd.Next() % (upper - lower + 1)) + lower; - } } } diff --git a/src/Numerics/Distributions/Discrete/Geometric.cs b/src/Numerics/Distributions/Discrete/Geometric.cs index 93b67cac..a0f01b5a 100644 --- a/src/Numerics/Distributions/Discrete/Geometric.cs +++ b/src/Numerics/Distributions/Discrete/Geometric.cs @@ -45,12 +45,12 @@ namespace MathNet.Numerics.Distributions /// /// The geometric distribution parameter. /// - private double _p; + double _p; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Geometric class. @@ -68,7 +68,7 @@ namespace MathNet.Numerics.Distributions /// /// The probability of generating a one. /// When the parameters don't pass the function. - private void SetParameters(double p) + void SetParameters(double p) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) { @@ -83,7 +83,7 @@ namespace MathNet.Numerics.Distributions /// /// The probability of generating a one. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double p) + static bool IsValidParameterSet(double p) { if (p >= 0.0 && p <= 1.0) { @@ -98,15 +98,9 @@ namespace MathNet.Numerics.Distributions /// public double P { - get - { - return _p; - } + get { return _p; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -127,10 +121,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -148,10 +139,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return 1.0 / _p; - } + get { return 1.0 / _p; } } /// @@ -159,10 +147,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return (1.0 - _p) / (_p * _p); - } + get { return (1.0 - _p) / (_p * _p); } } /// @@ -170,10 +155,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(1.0 - _p) / _p; - } + get { return Math.Sqrt(1.0 - _p) / _p; } } /// @@ -181,10 +163,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p; - } + get { return ((-_p * Math.Log(_p, 2.0)) - ((1.0 - _p) * Math.Log(1.0 - _p, 2.0))) / _p; } } /// @@ -193,10 +172,7 @@ namespace MathNet.Numerics.Distributions /// Throws a not supported exception. public double Skewness { - get - { - return (2.0 - _p) / Math.Sqrt(1.0 - _p); - } + get { return (2.0 - _p) / Math.Sqrt(1.0 - _p); } } /// @@ -218,10 +194,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return 1; - } + get { return 1; } } /// @@ -229,10 +202,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p)); - } + get { return (int)Math.Ceiling(-Constants.Ln2 / Math.Log(1 - _p)); } } /// @@ -240,10 +210,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 1; - } + get { return 1; } } /// @@ -251,10 +218,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return int.MaxValue; - } + get { return int.MaxValue; } } /// @@ -291,13 +255,28 @@ namespace MathNet.Numerics.Distributions return ((k - 1) * Math.Log(1.0 - _p)) + Math.Log(_p); } + #endregion + + /// + /// Returns one sample from the distribution. + /// + /// The random number generator to use. + /// The p parameter + /// + /// One sample from the distribution implied by . + /// + internal static int SampleUnchecked(Random rnd, double p) + { + return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p)); + } + /// /// Samples a Geometric distributed random variable. /// /// A sample from the Geometric distribution. public int Sample() { - return DoSample(RandomSource, _p); + return SampleUnchecked(RandomSource, _p); } /// @@ -308,23 +287,41 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _p); + yield return SampleUnchecked(RandomSource, _p); } } - #endregion + /// + /// Samples a random variable. + /// + /// The random number generator to use. + /// The p parameter + public static int Sample(Random rnd, double p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, p); + } /// - /// Returns one sample from the distribution. + /// Samples a sequence of this random variable. /// /// The random number generator to use. /// The p parameter - /// - /// One sample from the distribution implied by . - /// - private static int DoSample(Random rnd, double p) + public static IEnumerable Samples(Random rnd, double p) { - return p == 1.0 ? 1 : (int)Math.Ceiling(-Math.Log(1.0 - rnd.NextDouble(), 1.0 - p)); + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, p); + } } } } diff --git a/src/Numerics/Distributions/Discrete/Hypergeometric.cs b/src/Numerics/Distributions/Discrete/Hypergeometric.cs index 68ad0679..4901135b 100644 --- a/src/Numerics/Distributions/Discrete/Hypergeometric.cs +++ b/src/Numerics/Distributions/Discrete/Hypergeometric.cs @@ -52,22 +52,22 @@ namespace MathNet.Numerics.Distributions /// /// The size of the population. /// - private int _populationSize; + int _populationSize; /// /// The m parameter of the distribution. /// - private int _m; + int _m; /// /// The n parameter (number to draw) of the distribution. /// - private int _n; + int _n; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Hypergeometric class. @@ -87,7 +87,7 @@ namespace MathNet.Numerics.Distributions /// The Total parameter of the distribution. /// The m parameter of the distribution. /// The n parameter of the distribution. - private void SetParameters(int total, int m, int n) + void SetParameters(int total, int m, int n) { if (Control.CheckDistributionParameters && !IsValidParameterSet(total, m, n)) { @@ -106,7 +106,7 @@ namespace MathNet.Numerics.Distributions /// The m parameter of the distribution. /// The n parameter of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(int total, int m, int n) + static bool IsValidParameterSet(int total, int m, int n) { if (total < 0 || m < 0 || n < 0) { @@ -166,10 +166,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -187,10 +184,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return (double)_m * _n / _populationSize; - } + get { return (double)_m * _n / _populationSize; } } /// @@ -198,10 +192,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0)); - } + get { return _n * _m * (_populationSize - _n) * (_populationSize - _m) / (_populationSize * _populationSize * (_populationSize - 1.0)); } } /// @@ -225,10 +216,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0)); - } + get { return (Math.Sqrt(_populationSize - 1.0) * (_populationSize - (2 * _n)) * (_populationSize - (2 * _m))) / (Math.Sqrt(_n * _m * (_populationSize - _m) * (_populationSize - _n)) * (_populationSize - 2.0)); } } /// @@ -320,27 +308,6 @@ namespace MathNet.Numerics.Distributions return Math.Log(Probability(k)); } - /// - /// Samples a Hypergeometric distributed random variable. - /// - /// The number of successes in n trials. - public int Sample() - { - return DoSample(RandomSource, _populationSize, _m, _n); - } - - /// - /// Samples an array of Hypergeometric distributed random variables. - /// - /// a sequence of successes in n trials. - public IEnumerable Samples() - { - while (true) - { - yield return DoSample(RandomSource, _populationSize, _m, _n); - } - } - #endregion /// @@ -351,7 +318,7 @@ namespace MathNet.Numerics.Distributions /// The m parameter of the distribution. /// The n parameter of the distribution. /// a random number from the Hypergeometric distribution. - private static int DoSample(Random rnd, int size, int m, int n) + internal static int SampleUnchecked(Random rnd, int size, int m, int n) { var x = 0; @@ -367,10 +334,68 @@ namespace MathNet.Numerics.Distributions size--; n--; - } + } while (0 < n); return x; } + + /// + /// Samples a Hypergeometric distributed random variable. + /// + /// The number of successes in n trials. + public int Sample() + { + return SampleUnchecked(RandomSource, _populationSize, _m, _n); + } + + /// + /// Samples an array of Hypergeometric distributed random variables. + /// + /// a sequence of successes in n trials. + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(RandomSource, _populationSize, _m, _n); + } + } + + /// + /// Samples a random variable. + /// + /// The random number generator to use. + /// The population size. + /// The m parameter of the distribution. + /// The n parameter of the distribution. + public static int Sample(Random rnd, int populationSize, int m, int n) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, populationSize, m, n); + } + + /// + /// Samples a sequence of this random variable. + /// + /// The random number generator to use. + /// The population size. + /// The m parameter of the distribution. + /// The n parameter of the distribution. + public static IEnumerable Samples(Random rnd, int populationSize, int m, int n) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(populationSize, m, n)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, populationSize, m, n); + } + } } } diff --git a/src/Numerics/Distributions/Discrete/NegativeBinomial.cs b/src/Numerics/Distributions/Discrete/NegativeBinomial.cs index e9bfdf29..c619a974 100644 --- a/src/Numerics/Distributions/Discrete/NegativeBinomial.cs +++ b/src/Numerics/Distributions/Discrete/NegativeBinomial.cs @@ -46,32 +46,26 @@ namespace MathNet.Numerics.Distributions /// /// The r parameter of the distribution. /// - private double _r; + double _r; /// /// The p parameter of the distribution. /// - private double _p; + double _p; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Gets or sets the number of trials. /// public double R { - get - { - return _r; - } + get { return _r; } - set - { - SetParameters(value, _p); - } + set { SetParameters(value, _p); } } /// @@ -79,15 +73,9 @@ namespace MathNet.Numerics.Distributions /// public double P { - get - { - return _p; - } + get { return _p; } - set - { - SetParameters(_r, value); - } + set { SetParameters(_r, value); } } /// @@ -107,7 +95,7 @@ namespace MathNet.Numerics.Distributions /// The number of trials. /// The probability of a trial resulting in success. /// When the parameters don't pass the function. - private void SetParameters(double r, double p) + void SetParameters(double r, double p) { if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p)) { @@ -124,7 +112,7 @@ namespace MathNet.Numerics.Distributions /// The number of trials. /// The probability of a trial resulting in success. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double r, double p) + static bool IsValidParameterSet(double r, double p) { if (r < 0.0 || Double.IsNaN(r)) { @@ -157,10 +145,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -178,10 +163,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _r * (1.0 - _p) / _p; - } + get { return _r * (1.0 - _p) / _p; } } /// @@ -189,10 +171,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _r * (1.0 - _p) / (_p * _p); - } + get { return _r * (1.0 - _p) / (_p * _p); } } /// @@ -200,10 +179,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(_r * (1.0 - _p)) / _p; - } + get { return Math.Sqrt(_r * (1.0 - _p)) / _p; } } /// @@ -211,10 +187,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -222,10 +195,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p)); - } + get { return (2.0 - _p) / Math.Sqrt(_r * (1.0 - _p)); } } /// @@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0; - } + get { return _r > 1.0 ? (int)Math.Floor((_r - 1.0) * (1.0 - _p) / _p) : 0; } } /// @@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -269,10 +233,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -280,10 +241,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return int.MaxValue; - } + get { return int.MaxValue; } } /// @@ -296,10 +254,10 @@ namespace MathNet.Numerics.Distributions public double Probability(int k) { var ln = SpecialFunctions.GammaLn(_r + k) - - SpecialFunctions.GammaLn(_r) - - SpecialFunctions.GammaLn(k + 1.0) - + (_r * Math.Log(_p)) - + (k * Math.Log(1.0 - _p)); + - SpecialFunctions.GammaLn(_r) + - SpecialFunctions.GammaLn(k + 1.0) + + (_r * Math.Log(_p)) + + (k * Math.Log(1.0 - _p)); return Math.Exp(ln); } @@ -313,20 +271,44 @@ namespace MathNet.Numerics.Distributions public double ProbabilityLn(int k) { var ln = SpecialFunctions.GammaLn(_r + k) - - SpecialFunctions.GammaLn(_r) - - SpecialFunctions.GammaLn(k + 1.0) - + (_r * Math.Log(_p)) - + (k * Math.Log(1.0 - _p)); + - SpecialFunctions.GammaLn(_r) + - SpecialFunctions.GammaLn(k + 1.0) + + (_r * Math.Log(_p)) + + (k * Math.Log(1.0 - _p)); return ln; } + #endregion + + /// + /// Samples a negative binomial distributed random variable. + /// + /// The random number generator to use. + /// The r parameter. + /// The p parameter. + /// a sample from the distribution. + internal static int SampleUnchecked(Random rnd, double r, double p) + { + var lambda = Gamma.SampleUnchecked(rnd, r, p); + var c = Math.Exp(-lambda); + var p1 = 1.0; + var k = 0; + do + { + k = k + 1; + p1 = p1 * rnd.NextDouble(); + } + while (p1 >= c); + return k - 1; + } + /// /// Samples a NegativeBinomial distributed random variable. /// /// a sample from the distribution. public int Sample() { - return DoSample(RandomSource, _r, _p); + return SampleUnchecked(RandomSource, _r, _p); } /// @@ -337,43 +319,43 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _r, _p); + yield return SampleUnchecked(RandomSource, _r, _p); } } - #endregion - /// - /// Samples a negative binomial distributed random variable. + /// Samples a random variable. /// /// The random number generator to use. /// The r parameter. /// The p parameter. - /// a sample from the distribution. - private static int DoSample(Random rnd, double r, double p) + public static int Sample(Random rnd, double r, double p) { - var lambda = Gamma.Sample(rnd, r, p); - return DoPoissonSample(rnd, lambda); + if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, r, p); } /// - /// Use Knuth's method to generate Poisson distributed random variables. + /// Samples a sequence of this random variable. /// - /// The random number generator. - /// The lambda value. - /// a Poisson distributed random number. - private static int DoPoissonSample(Random rnd, double lambda) + /// The random number generator to use. + /// The r parameter. + /// The p parameter. + public static IEnumerable Samples(Random rnd, double r, double p) { - var c = Math.Exp(-lambda); - var p = 1.0; - var k = 0; - do + if (Control.CheckDistributionParameters && !IsValidParameterSet(r, p)) { - k = k + 1; - p = p * rnd.NextDouble(); + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, r, p); } - while (p >= c); - return k - 1; } } } diff --git a/src/Numerics/Distributions/Discrete/Poisson.cs b/src/Numerics/Distributions/Discrete/Poisson.cs index f3fc4b22..39b3d432 100644 --- a/src/Numerics/Distributions/Discrete/Poisson.cs +++ b/src/Numerics/Distributions/Discrete/Poisson.cs @@ -43,27 +43,21 @@ namespace MathNet.Numerics.Distributions /// /// The Poisson distribution parameter λ. /// - private double _lambda; + double _lambda; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Gets or sets the Poisson distribution parameter λ. /// public double Lambda { - get - { - return _lambda; - } + get { return _lambda; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -82,7 +76,7 @@ namespace MathNet.Numerics.Distributions /// /// The mean (λ) of the distribution. /// When the parameters don't pass the function. - private void SetParameters(double lambda) + void SetParameters(double lambda) { if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) { @@ -97,7 +91,7 @@ namespace MathNet.Numerics.Distributions /// /// The mean (λ) of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double lambda) + static bool IsValidParameterSet(double lambda) { return lambda > 0.00; } @@ -120,10 +114,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -141,10 +132,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _lambda; - } + get { return _lambda; } } /// @@ -152,10 +140,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _lambda; - } + get { return _lambda; } } /// @@ -163,10 +148,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(_lambda); - } + get { return Math.Sqrt(_lambda); } } /// @@ -175,10 +157,7 @@ namespace MathNet.Numerics.Distributions /// Approximation, see Wikipedia Poisson distribution public double Entropy { - get - { - return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda)); - } + get { return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda)); } } /// @@ -186,10 +165,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return 1.0 / Math.Sqrt(_lambda); - } + get { return 1.0 / Math.Sqrt(_lambda); } } /// @@ -197,10 +173,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 0; - } + get { return 0; } } /// @@ -208,10 +181,7 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return int.MaxValue; - } + get { return int.MaxValue; } } /// @@ -233,10 +203,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return (int)Math.Floor(_lambda); - } + get { return (int)Math.Floor(_lambda); } } /// @@ -245,10 +212,7 @@ namespace MathNet.Numerics.Distributions /// Approximation, see Wikipedia Poisson distribution public int Median { - get - { - return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda)); - } + get { return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda)); } } /// @@ -271,71 +235,15 @@ namespace MathNet.Numerics.Distributions return -_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k); } - /// - /// Samples a Poisson distributed random variable. - /// - /// A sample from the Poisson distribution. - public int Sample() - { - return DoSample(RandomSource, _lambda); - } - - /// - /// Samples an array of Poisson distributed random variables. - /// - /// a sequence of successes in N trials. - public IEnumerable Samples() - { - while (true) - { - yield return DoSample(RandomSource, _lambda); - } - } - #endregion - /// - /// Samples a Poisson distributed random variable. - /// - /// The random number generator to use. - /// The Poisson distribution parameter λ. - /// A sample from the Poisson distribution. - public static int Sample(Random rnd, double lambda) - { - if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } - - return DoSample(rnd, lambda); - } - - /// - /// Samples a sequence of Poisson distributed random variables. - /// - /// The random number generator to use. - /// The Poisson distribution parameter λ. - /// a sequence of samples from the distribution. - public static IEnumerable Samples(Random rnd, double lambda) - { - if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } - - while (true) - { - yield return DoSample(rnd, lambda); - } - } - /// /// Generates one sample from the Poisson distribution. /// /// The random source to use. /// The Poisson distribution parameter λ. /// A random sample from the Poisson distribution. - private static int DoSample(Random rnd, double lambda) + internal static int SampleUnchecked(Random rnd, double lambda) { return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda); } @@ -346,7 +254,7 @@ namespace MathNet.Numerics.Distributions /// The random source to use. /// The Poisson distribution parameter λ. /// A random sample from the Poisson distribution. - private static int DoSampleShort(Random rnd, double lambda) + static int DoSampleShort(Random rnd, double lambda) { var limit = Math.Exp(-lambda); var count = 0; @@ -367,7 +275,7 @@ namespace MathNet.Numerics.Distributions /// "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson, /// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979) /// The article is on pages 29-35. The algorithm given here is on page 32. - private static int DoSampleLarge(Random rnd, double lambda) + static int DoSampleLarge(Random rnd, double lambda) { var c = 0.767 - (3.36 / lambda); var beta = Math.PI / Math.Sqrt(3.0 * lambda); @@ -395,5 +303,61 @@ namespace MathNet.Numerics.Distributions } } } + + /// + /// Samples a Poisson distributed random variable. + /// + /// A sample from the Poisson distribution. + public int Sample() + { + return SampleUnchecked(RandomSource, _lambda); + } + + /// + /// Samples an array of Poisson distributed random variables. + /// + /// a sequence of successes in N trials. + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(RandomSource, _lambda); + } + } + + /// + /// Samples a Poisson distributed random variable. + /// + /// The random number generator to use. + /// The Poisson distribution parameter λ. + /// A sample from the Poisson distribution. + public static int Sample(Random rnd, double lambda) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, lambda); + } + + /// + /// Samples a sequence of Poisson distributed random variables. + /// + /// The random number generator to use. + /// The Poisson distribution parameter λ. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double lambda) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, lambda); + } + } } } diff --git a/src/Numerics/Distributions/Discrete/Zipf.cs b/src/Numerics/Distributions/Discrete/Zipf.cs index 34849ea7..8bd0bc88 100644 --- a/src/Numerics/Distributions/Discrete/Zipf.cs +++ b/src/Numerics/Distributions/Discrete/Zipf.cs @@ -47,17 +47,17 @@ namespace MathNet.Numerics.Distributions /// /// The s parameter of the distribution. /// - private double _s; + double _s; /// /// The n parameter of the distribution. /// - private int _n; + int _n; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -79,7 +79,7 @@ namespace MathNet.Numerics.Distributions /// /// The s parameter of the distribution. /// The n parameter of the distribution. - private void SetParameters(double s, int n) + void SetParameters(double s, int n) { if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n)) { @@ -96,7 +96,7 @@ namespace MathNet.Numerics.Distributions /// The s parameter of the distribution. /// The n parameter of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double s, int n) + static bool IsValidParameterSet(double s, int n) { if (n <= 0 || s <= 0 || Double.IsNaN(s)) { @@ -111,15 +111,9 @@ namespace MathNet.Numerics.Distributions /// public double S { - get - { - return _s; - } + get { return _s; } - set - { - SetParameters(value, _n); - } + set { SetParameters(value, _n); } } /// @@ -127,15 +121,9 @@ namespace MathNet.Numerics.Distributions /// public int N { - get - { - return _n; - } + get { return _n; } - set - { - SetParameters(_s, value); - } + set { SetParameters(_s, value); } } /// @@ -154,10 +142,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -175,10 +160,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s); - } + get { return SpecialFunctions.GeneralHarmonic(_n, _s - 1.0) / SpecialFunctions.GeneralHarmonic(_n, _s); } } /// @@ -203,10 +185,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(Variance); - } + get { return Math.Sqrt(Variance); } } /// @@ -266,10 +245,7 @@ namespace MathNet.Numerics.Distributions /// public int Mode { - get - { - return 1; - } + get { return 1; } } /// @@ -277,10 +253,7 @@ namespace MathNet.Numerics.Distributions /// public int Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -288,10 +261,7 @@ namespace MathNet.Numerics.Distributions /// public int Minimum { - get - { - return 1; - } + get { return 1; } } /// @@ -299,12 +269,9 @@ namespace MathNet.Numerics.Distributions /// public int Maximum { - get - { - return _n; - } + get { return _n; } } - + /// /// Computes values of the probability mass function. /// @@ -329,13 +296,45 @@ namespace MathNet.Numerics.Distributions return Math.Log(Probability(k)); } + #endregion + + /// + /// Generates a sample from the Zipf distribution without doing parameter checking. + /// + /// The random number generator to use. + /// The s parameter of the distribution. + /// The n parameter of the distribution. + /// a random number from the Zipf distribution. + internal static int SampleUnchecked(Random rnd, double s, int n) + { + var r = 0.0; + while (r == 0.0) + { + r = rnd.NextDouble(); + } + + var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s); + int i; + var sum = 0.0; + for (i = 1; i <= n; i++) + { + sum += p / Math.Pow(i, s); + if (sum >= r) + { + break; + } + } + + return i; + } + /// /// Draws a random sample from the distribution. /// /// a sample from the distribution. public int Sample() { - return DoSample(RandomSource, _s, _n); + return SampleUnchecked(RandomSource, _s, _n); } /// @@ -346,40 +345,43 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _s, _n); + yield return SampleUnchecked(RandomSource, _s, _n); } } - #endregion - /// - /// Generates a sample from the Zipf distribution without doing parameter checking. + /// Samples a random variable. /// /// The random number generator to use. /// The s parameter of the distribution. /// The n parameter of the distribution. - /// a random number from the Zipf distribution. - private static int DoSample(Random rnd, double s, int n) + public static int Sample(Random rnd, double s, int n) { - var r = 0.0; - while (r == 0.0) + if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n)) { - r = rnd.NextDouble(); + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s); - int i; - var sum = 0.0; - for (i = 1; i <= n; i++) + return SampleUnchecked(rnd, s, n); + } + + /// + /// Samples a sequence of this random variable. + /// + /// The random number generator to use. + /// The s parameter of the distribution. + /// The n parameter of the distribution. + public static IEnumerable Samples(Random rnd, double s, int n) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(s, n)) { - sum += p / Math.Pow(i, s); - if (sum >= r) - { - break; - } + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return i; + while (true) + { + yield return SampleUnchecked(rnd, s, n); + } } } } diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs index ba424745..6816902d 100644 --- a/src/Numerics/Distributions/Multivariate/Multinomial.cs +++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs @@ -368,7 +368,7 @@ namespace MathNet.Numerics.Distributions for (var i = 0; i < n; i++) { - ret[Categorical.DoSample(rnd, cp)]++; + ret[Categorical.SampleUnchecked(rnd, cp)]++; } return ret; @@ -399,7 +399,7 @@ namespace MathNet.Numerics.Distributions for (var i = 0; i < n; i++) { - ret[Categorical.DoSample(rnd, cp)]++; + ret[Categorical.SampleUnchecked(rnd, cp)]++; } yield return ret;