From 6c32ba172b7490221ddb53e381154bd6f5dc5798 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 25 Nov 2012 14:23:52 +0100 Subject: [PATCH] Distributions: consistent sample methods for continuous distributions #32 --- src/Numerics/Distributions/Continuous/Beta.cs | 140 +++++------- .../Distributions/Continuous/Cauchy.cs | 136 ++++++----- src/Numerics/Distributions/Continuous/Chi.cs | 113 +++++----- .../Distributions/Continuous/ChiSquare.cs | 91 ++++---- .../Continuous/ContinuousUniform.cs | 79 +++---- .../Distributions/Continuous/Erlang.cs | 174 +++++++------- .../Distributions/Continuous/Exponential.cs | 115 ++++------ .../Continuous/FisherSnedecor.cs | 111 ++++----- .../Distributions/Continuous/Gamma.cs | 213 ++++++++---------- .../Distributions/Continuous/InverseGamma.cs | 126 +++++------ .../Distributions/Continuous/Laplace.cs | 135 ++++++----- .../Distributions/Continuous/LogNormal.cs | 77 ++----- .../Distributions/Continuous/Normal.cs | 179 +++++++-------- .../Distributions/Continuous/Pareto.cs | 123 +++++----- .../Distributions/Continuous/Rayleigh.cs | 118 +++++----- .../Distributions/Continuous/Stable.cs | 164 +++++++------- .../Distributions/Continuous/StudentT.cs | 131 ++++------- .../Distributions/Continuous/Weibull.cs | 103 +++------ .../Multivariate/MatrixNormal.cs | 2 +- 19 files changed, 1065 insertions(+), 1265 deletions(-) diff --git a/src/Numerics/Distributions/Continuous/Beta.cs b/src/Numerics/Distributions/Continuous/Beta.cs index 374450a1..bffd2db1 100644 --- a/src/Numerics/Distributions/Continuous/Beta.cs +++ b/src/Numerics/Distributions/Continuous/Beta.cs @@ -51,17 +51,17 @@ namespace MathNet.Numerics.Distributions /// /// Beta shape parameter a. /// - private double _shapeA; + double _shapeA; /// /// Beta shape parameter b. /// - private double _shapeB; + double _shapeB; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Beta class. @@ -90,7 +90,7 @@ namespace MathNet.Numerics.Distributions /// The a shape parameter of the Beta distribution. /// The b shape parameter of the Beta distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double a, double b) + static bool IsValidParameterSet(double a, double b) { if (a < 0.0 || b < 0.0 || Double.IsNaN(a) || Double.IsNaN(b)) { @@ -106,7 +106,7 @@ namespace MathNet.Numerics.Distributions /// The a shape parameter of the Beta distribution. /// The b shape parameter of the Beta distribution. /// When the parameters don't pass the function. - private void SetParameters(double a, double b) + void SetParameters(double a, double b) { if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) { @@ -122,15 +122,8 @@ namespace MathNet.Numerics.Distributions /// public double A { - get - { - return _shapeA; - } - - set - { - SetParameters(value, _shapeB); - } + get { return _shapeA; } + set { SetParameters(value, _shapeB); } } /// @@ -138,15 +131,8 @@ namespace MathNet.Numerics.Distributions /// public double B { - get - { - return _shapeB; - } - - set - { - SetParameters(_shapeA, value); - } + get { return _shapeB; } + set { SetParameters(_shapeA, value); } } #region IDistribution implementation @@ -156,10 +142,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -218,10 +201,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return (_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0)); - } + get { return (_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0)); } } /// @@ -229,10 +209,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt((_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0))); - } + get { return Math.Sqrt((_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0))); } } /// @@ -258,9 +235,9 @@ namespace MathNet.Numerics.Distributions } return SpecialFunctions.BetaLn(_shapeA, _shapeB) - - ((_shapeA - 1.0) * SpecialFunctions.DiGamma(_shapeA)) - - ((_shapeB - 1.0) * SpecialFunctions.DiGamma(_shapeB)) - + ((_shapeA + _shapeB - 2.0) * SpecialFunctions.DiGamma(_shapeA + _shapeB)); + - ((_shapeA - 1.0) * SpecialFunctions.DiGamma(_shapeA)) + - ((_shapeB - 1.0) * SpecialFunctions.DiGamma(_shapeB)) + + ((_shapeA + _shapeB - 2.0) * SpecialFunctions.DiGamma(_shapeA + _shapeB)); } } @@ -302,7 +279,7 @@ namespace MathNet.Numerics.Distributions } return 2.0 * (_shapeB - _shapeA) * Math.Sqrt(_shapeA + _shapeB + 1.0) - / ((_shapeA + _shapeB + 2.0) * Math.Sqrt(_shapeA * _shapeB)); + / ((_shapeA + _shapeB + 2.0) * Math.Sqrt(_shapeA * _shapeB)); } } @@ -361,10 +338,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -372,10 +346,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -383,10 +354,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return 1.0; - } + get { return 1.0; } } /// @@ -515,27 +483,27 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x >= 1.0) { return 1.0; } - + if (Double.IsPositiveInfinity(_shapeA) && Double.IsPositiveInfinity(_shapeB)) { return x < 0.5 ? 0.0 : 1.0; } - + if (Double.IsPositiveInfinity(_shapeA)) { return x < 1.0 ? 0.0 : 1.0; } - + if (Double.IsPositiveInfinity(_shapeB)) { return x >= 0.0 ? 1.0 : 0.0; } - + if (_shapeA == 0.0 && _shapeB == 0.0) { if (x >= 0.0 && x < 1.0) @@ -545,32 +513,48 @@ namespace MathNet.Numerics.Distributions return 1.0; } - + if (_shapeA == 0.0) { return 1.0; } - + if (_shapeB == 0.0) { return x >= 1.0 ? 1.0 : 0.0; } - + if (_shapeA == 1.0 && _shapeB == 1.0) { return x; } - + return SpecialFunctions.BetaRegularized(_shapeA, _shapeB, x); } + #endregion + + /// + /// Samples Beta distributed random variables by sampling two Gamma variables and normalizing. + /// + /// The random number generator to use. + /// The A shape parameter. + /// The B shape parameter. + /// a random number from the Beta distribution. + internal static double SampleUnchecked(Random rnd, double a, double b) + { + var x = Gamma.SampleUnchecked(rnd, a, 1.0); + var y = Gamma.SampleUnchecked(rnd, b, 1.0); + return x / (x + y); + } + /// /// Generates a sample from the Beta distribution. /// /// a sample from the distribution. public double Sample() { - return SampleBeta(RandomSource, _shapeA, _shapeB); + return SampleUnchecked(RandomSource, _shapeA, _shapeB); } /// @@ -581,37 +565,35 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return SampleBeta(RandomSource, _shapeA, _shapeB); + yield return SampleUnchecked(RandomSource, _shapeA, _shapeB); } } - #endregion - /// - /// Generates a sample from the normal distribution using the Box-Muller algorithm. + /// Generates a sample from the distribution. /// - /// The random number generator to use. + /// The random number generator to use. /// The a shape parameter of the Beta distribution. /// The b shape parameter of the Beta distribution. /// a sample from the distribution. - public static double Sample(Random rng, double a, double b) + public static double Sample(Random rnd, double a, double b) { if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return SampleBeta(rng, a, b); + return SampleUnchecked(rnd, a, b); } /// - /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. + /// Generates a sequence of samples from the distribution. /// - /// The random number generator to use. + /// The random number generator to use. /// The a shape parameter of the Beta distribution. /// The b shape parameter of the Beta distribution. /// a sequence of samples from the distribution. - public static IEnumerable Samples(Random rng, double a, double b) + public static IEnumerable Samples(Random rnd, double a, double b) { if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) { @@ -620,22 +602,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return SampleBeta(rng, a, b); + yield return SampleUnchecked(rnd, a, b); } } - - /// - /// Samples Beta distributed random variables by sampling two Gamma variables and normalizing. - /// - /// The random number generator to use. - /// The A shape parameter. - /// The B shape parameter. - /// a random number from the Beta distribution. - internal static double SampleBeta(Random rnd, double a, double b) - { - var x = Gamma.SampleGamma(rnd, a, 1.0); - var y = Gamma.SampleGamma(rnd, b, 1.0); - return x / (x + y); - } } } diff --git a/src/Numerics/Distributions/Continuous/Cauchy.cs b/src/Numerics/Distributions/Continuous/Cauchy.cs index 531b6163..a8e206b4 100644 --- a/src/Numerics/Distributions/Continuous/Cauchy.cs +++ b/src/Numerics/Distributions/Continuous/Cauchy.cs @@ -44,17 +44,18 @@ namespace MathNet.Numerics.Distributions /// /// The scale of the Cauchy distribution. /// - private double _scale; + double _scale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class with the location parameter set to 0 and the scale parameter set to 1 /// - public Cauchy() : this(0, 1) + public Cauchy() + : this(0, 1) { } @@ -82,7 +83,7 @@ namespace MathNet.Numerics.Distributions /// Location parameter. /// Scale parameter. Must be greater than 0. /// When the parameters don't pass the function. - private void SetParameters(double location, double scale) + void SetParameters(double location, double scale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) { @@ -99,7 +100,7 @@ namespace MathNet.Numerics.Distributions /// Location parameter. /// Scale parameter. Must be greater than 0. /// True when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double location, double scale) + static bool IsValidParameterSet(double location, double scale) { if (scale <= 0) { @@ -119,15 +120,9 @@ namespace MathNet.Numerics.Distributions /// public double Location { - get - { - return Median; - } + get { return Median; } - set - { - SetParameters(value, _scale); - } + set { SetParameters(value, _scale); } } /// @@ -135,15 +130,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(Median, value); - } + set { SetParameters(Median, value); } } /// @@ -162,11 +151,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } - + get { return _random; } set { if (value == null) @@ -183,10 +168,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -194,10 +176,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -205,10 +184,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -216,10 +192,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return Math.Log(4.0 * Constants.Pi * _scale); - } + get { return Math.Log(4.0 * Constants.Pi * _scale); } } /// @@ -227,10 +200,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -252,30 +222,20 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return Median; - } + get { return Median; } } /// /// Gets the median of the distribution. /// - public double Median - { - get; - private set; - } + public double Median { get; private set; } /// /// Gets the minimum of the distribution. /// public double Minimum { - get - { - return Double.NegativeInfinity; - } + get { return Double.NegativeInfinity; } } /// @@ -283,10 +243,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -309,13 +266,28 @@ namespace MathNet.Numerics.Distributions return -Math.Log(Constants.Pi * _scale * (1.0 + (((x - Median) / _scale) * ((x - Median) / _scale)))); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The location shape parameter. + /// The scale parameter. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double location, double scale) + { + var u = rnd.NextDouble(); + return location + (scale * Math.Tan(Constants.Pi * (u - 0.5))); + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource, Median, _scale); + return SampleUnchecked(RandomSource, Median, _scale); } /// @@ -326,23 +298,45 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, Median, _scale); + yield return SampleUnchecked(RandomSource, Median, _scale); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The location shape parameter. + /// The scale parameter. + /// a sample from the distribution. + public static double Sample(Random rnd, double location, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, location, scale); + } /// - /// Samples the distribution. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. /// The location shape parameter. /// The scale parameter. - /// a random number from the distribution. - private static double DoSample(Random rnd, double location, double scale) + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double location, double scale) { - var u = rnd.NextDouble(); - return location + (scale * Math.Tan(Constants.Pi * (u - 0.5))); + if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, location, scale); + } } } } diff --git a/src/Numerics/Distributions/Continuous/Chi.cs b/src/Numerics/Distributions/Continuous/Chi.cs index 62c7f441..c7ac9073 100644 --- a/src/Numerics/Distributions/Continuous/Chi.cs +++ b/src/Numerics/Distributions/Continuous/Chi.cs @@ -47,12 +47,12 @@ namespace MathNet.Numerics.Distributions /// /// Keeps track of the degrees of freedom for the Chi distribution. /// - private double _dof; + double _dof; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -71,7 +71,7 @@ namespace MathNet.Numerics.Distributions /// /// The degrees of freedom for the Chi distribution. /// When the parameters don't pass the function. - private void SetParameters(double dof) + void SetParameters(double dof) { if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) { @@ -86,7 +86,7 @@ namespace MathNet.Numerics.Distributions /// /// The degrees of freedom for the Chi distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double dof) + static bool IsValidParameterSet(double dof) { if (dof <= 0 || Double.IsNaN(dof)) { @@ -101,15 +101,9 @@ namespace MathNet.Numerics.Distributions /// public double DegreesOfFreedom { - get - { - return _dof; - } + get { return _dof; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -128,10 +122,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -149,10 +140,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return Math.Sqrt(2) * (SpecialFunctions.Gamma((_dof + 1.0) / 2.0) / SpecialFunctions.Gamma(_dof / 2.0)); - } + get { return Math.Sqrt(2) * (SpecialFunctions.Gamma((_dof + 1.0) / 2.0) / SpecialFunctions.Gamma(_dof / 2.0)); } } /// @@ -160,10 +148,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _dof - (Mean * Mean); - } + get { return _dof - (Mean * Mean); } } /// @@ -171,10 +156,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(Variance); - } + get { return Math.Sqrt(Variance); } } /// @@ -182,10 +164,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return SpecialFunctions.GammaLn(_dof / 2.0) + ((_dof - Math.Log(2) - ((_dof - 1.0) * SpecialFunctions.DiGamma(_dof / 2.0))) / 2.0); - } + get { return SpecialFunctions.GammaLn(_dof / 2.0) + ((_dof - Math.Log(2) - ((_dof - 1.0) * SpecialFunctions.DiGamma(_dof / 2.0))) / 2.0); } } /// @@ -235,10 +214,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -246,10 +222,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -257,10 +230,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -283,13 +253,31 @@ namespace MathNet.Numerics.Distributions return ((1.0 - (_dof / 2.0)) * Math.Log(2.0)) + ((_dof - 1.0) * Math.Log(x)) - (x * x / 2.0) - SpecialFunctions.GammaLn(_dof / 2.0); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, int dof) + { + double sum = 0; + for (var i = 0; i < dof; i++) + { + sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2); + } + + return Math.Sqrt(sum); + } + /// /// Generates a sample from the Chi distribution. /// /// a sample from the distribution. public double Sample() { - return DoSample(RandomSource); + return SampleUnchecked(RandomSource, (int)_dof); } /// @@ -298,29 +286,44 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public IEnumerable Samples() { + var dof = (int)_dof; while (true) { - yield return DoSample(RandomSource); + yield return SampleUnchecked(RandomSource, dof); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// a sample from the distribution. + public static double Sample(Random rnd, int dof) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, dof); + } /// - /// Samples the distribution. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. - /// a random number from the distribution. - private double DoSample(Random rnd) + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, int dof) { - double sum = 0; - var n = (int)_dof; - for (var i = 0; i < n; i++) + if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) { - sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2); + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return Math.Sqrt(sum); + while (true) + { + yield return SampleUnchecked(rnd, dof); + } } } } diff --git a/src/Numerics/Distributions/Continuous/ChiSquare.cs b/src/Numerics/Distributions/Continuous/ChiSquare.cs index cd05120f..1c9a2d0a 100644 --- a/src/Numerics/Distributions/Continuous/ChiSquare.cs +++ b/src/Numerics/Distributions/Continuous/ChiSquare.cs @@ -49,7 +49,7 @@ namespace MathNet.Numerics.Distributions /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -68,7 +68,7 @@ namespace MathNet.Numerics.Distributions /// /// The degrees of freedom for the ChiSquare distribution. /// When the parameters don't pass the function. - private void SetParameters(double dof) + void SetParameters(double dof) { if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) { @@ -83,7 +83,7 @@ namespace MathNet.Numerics.Distributions /// /// The degrees of freedom for the ChiSquare distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double dof) + static bool IsValidParameterSet(double dof) { return dof > 0 && !Double.IsNaN(dof); } @@ -93,15 +93,9 @@ namespace MathNet.Numerics.Distributions /// public double DegreesOfFreedom { - get - { - return Mean; - } + get { return Mean; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -120,11 +114,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } - + get { return _random; } set { if (value == null) @@ -139,11 +129,7 @@ namespace MathNet.Numerics.Distributions /// /// Gets the mean of the distribution. /// - public double Mean - { - get; - private set; - } + public double Mean { get; private set; } /// /// Gets the variance of the distribution. @@ -243,13 +229,39 @@ namespace MathNet.Numerics.Distributions return (-x / 2.0) + (((Mean / 2.0) - 1.0) * Math.Log(x)) - ((Mean / 2.0) * Math.Log(2)) - SpecialFunctions.GammaLn(Mean / 2.0); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The degrees of freedom. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double dof) + { + //Use the simple method if the dof is an integer anyway + if (Math.Floor(dof) == dof && dof < Int32.MaxValue) + { + double sum = 0; + var n = (int)dof; + for (var i = 0; i < n; i++) + { + sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2); + } + return sum; + } + //Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations + //for a justification) + return Gamma.SampleUnchecked(rnd, dof / 2.0, .5); + } + /// /// Generates a sample from the ChiSquare distribution. /// /// a sample from the distribution. public double Sample() { - return DoSample(RandomSource, Mean); + return SampleUnchecked(RandomSource, Mean); } /// @@ -260,50 +272,43 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, Mean); + yield return SampleUnchecked(RandomSource, Mean); } } - #endregion - /// - /// Samples the distribution. + /// Generates a sample from the ChiSquare distribution. /// /// The random number generator to use. /// The degrees of freedom. - /// a random number from the distribution. - private static double DoSample(Random rnd, double dof) + /// a sample from the distribution. + public static double Sample(Random rnd, double dof) { - //Use the simple method if the dof is an integer anyway - if (Math.Floor(dof) == dof && dof < Int32.MaxValue) + if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) { - double sum = 0; - var n = (int)dof; - for (var i = 0; i < n; i++) - { - sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2); - } - return sum; + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - //Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations - //for a justification) - return Gamma.Sample(rnd, dof / 2.0, .5); + + return SampleUnchecked(rnd, dof); } /// - /// Generates a sample from the ChiSquare distribution. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. /// The degrees of freedom. /// a sample from the distribution. - public static double Sample(Random rnd, double dof) + public static IEnumerable Samples(Random rnd, double dof) { if (Control.CheckDistributionParameters && !IsValidParameterSet(dof)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return DoSample(rnd, dof); + while (true) + { + yield return SampleUnchecked(rnd, dof); + } } } } diff --git a/src/Numerics/Distributions/Continuous/ContinuousUniform.cs b/src/Numerics/Distributions/Continuous/ContinuousUniform.cs index 25701c89..b8314a5d 100644 --- a/src/Numerics/Distributions/Continuous/ContinuousUniform.cs +++ b/src/Numerics/Distributions/Continuous/ContinuousUniform.cs @@ -48,22 +48,23 @@ namespace MathNet.Numerics.Distributions /// /// The distribution's lower bound. /// - private double _lower; + double _lower; /// /// The distribution's upper bound. /// - private double _upper; + double _upper; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the ContinuousUniform class with lower bound 0 and upper bound 1. /// - public ContinuousUniform() : this(0.0, 1.0) + public ContinuousUniform() + : this(0.0, 1.0) { } @@ -94,13 +95,13 @@ 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(double lower, double upper) + static bool IsValidParameterSet(double lower, double upper) { if (upper < lower) { return false; } - + if (Double.IsNaN(upper) || Double.IsNaN(lower)) { return false; @@ -115,7 +116,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(double lower, double upper) + void SetParameters(double lower, double upper) { if (Control.CheckDistributionParameters && !IsValidParameterSet(lower, upper)) { @@ -131,15 +132,9 @@ namespace MathNet.Numerics.Distributions /// public double Lower { - get - { - return _lower; - } + get { return _lower; } - set - { - SetParameters(value, _upper); - } + set { SetParameters(value, _upper); } } /// @@ -147,15 +142,9 @@ namespace MathNet.Numerics.Distributions /// public double Upper { - get - { - return _upper; - } + get { return _upper; } - set - { - SetParameters(_lower, value); - } + set { SetParameters(_lower, value); } } #region IDistribution Members @@ -165,10 +154,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -221,6 +207,7 @@ namespace MathNet.Numerics.Distributions { get { return 0.0; } } + #endregion #region IContinuousDistribution Members @@ -300,7 +287,7 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (x >= _upper) { return 1.0; @@ -309,13 +296,27 @@ namespace MathNet.Numerics.Distributions return (x - _lower) / (_upper - _lower); } + #endregion + + /// + /// Generates one sample from the ContinuousUniform distribution without parameter checking. + /// + /// The random number generator to use. + /// The lower bound of the uniform random variable. + /// The upper bound of the uniform random variable. + /// a uniformly distributed random number. + internal static double SampleUnchecked(Random rnd, double lower, double upper) + { + return lower + (rnd.NextDouble() * (upper - lower)); + } + /// /// Generates a sample from the ContinuousUniform distribution. /// /// a sample from the distribution. public double Sample() { - return DoSample(RandomSource, _lower, _upper); + return SampleUnchecked(RandomSource, _lower, _upper); } /// @@ -326,12 +327,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _lower, _upper); + yield return SampleUnchecked(RandomSource, _lower, _upper); } } - #endregion - /// /// Generates a sample from the ContinuousUniform distribution. /// @@ -346,7 +345,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return DoSample(rnd, lower, upper); + return SampleUnchecked(rnd, lower, upper); } /// @@ -365,20 +364,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return DoSample(rnd, lower, upper); + yield return SampleUnchecked(rnd, lower, upper); } } - - /// - /// Generates one sample from the ContinuousUniform distribution without parameter checking. - /// - /// The random number generator to use. - /// The lower bound of the uniform random variable. - /// The upper bound of the uniform random variable. - /// a uniformly distributed random number. - private static double DoSample(Random rnd, double lower, double upper) - { - return lower + (rnd.NextDouble() * (upper - lower)); - } } } \ No newline at end of file diff --git a/src/Numerics/Distributions/Continuous/Erlang.cs b/src/Numerics/Distributions/Continuous/Erlang.cs index 8e8057a4..bfb220ae 100644 --- a/src/Numerics/Distributions/Continuous/Erlang.cs +++ b/src/Numerics/Distributions/Continuous/Erlang.cs @@ -46,17 +46,17 @@ namespace MathNet.Numerics.Distributions /// /// Erlang shape parameter. /// - private double _shape; + double _shape; /// /// Erlang inverse scale parameter. /// - private double _invScale; + double _invScale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -102,7 +102,7 @@ namespace MathNet.Numerics.Distributions /// /// The shape of the Erlang distribution. /// The inverse scale of the Erlang distribution. - private void SetParameters(double shape, double invScale) + void SetParameters(double shape, double invScale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) { @@ -119,7 +119,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Erlang distribution. /// The inverse scale of the Erlang distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double shape, double invScale) + static bool IsValidParameterSet(double shape, double invScale) { if (shape < 0.0 || invScale < 0.0 || Double.IsNaN(shape) || Double.IsNaN(invScale)) { @@ -134,15 +134,9 @@ namespace MathNet.Numerics.Distributions /// public int Shape { - get - { - return (int)_shape; - } + get { return (int)_shape; } - set - { - SetParameters(value, _invScale); - } + set { SetParameters(value, _invScale); } } /// @@ -150,10 +144,7 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return 1.0 / _invScale; - } + get { return 1.0 / _invScale; } set { @@ -173,15 +164,9 @@ namespace MathNet.Numerics.Distributions /// public double InvScale { - get - { - return _invScale; - } + get { return _invScale; } - set - { - SetParameters(_shape, value); - } + set { SetParameters(_shape, value); } } /// @@ -200,10 +185,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -227,12 +209,12 @@ namespace MathNet.Numerics.Distributions { return _shape; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return _shape / _invScale; } } @@ -248,12 +230,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return _shape / (_invScale * _invScale); } } @@ -269,12 +251,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return Math.Sqrt(_shape) / _invScale; } } @@ -295,7 +277,7 @@ namespace MathNet.Numerics.Distributions { return Double.NaN; } - + return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + ((1.0 - _shape) * SpecialFunctions.DiGamma(_shape)); } } @@ -311,12 +293,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return 2.0 / Math.Sqrt(_shape); } } @@ -332,12 +314,12 @@ namespace MathNet.Numerics.Distributions { return x >= _shape ? 1.0 : 0.0; } - + if (_shape == 0.0 && _invScale == 0.0) { return 0.0; } - + return SpecialFunctions.GammaLowerRegularized(_shape, x * _invScale); } @@ -361,12 +343,12 @@ namespace MathNet.Numerics.Distributions { return _shape; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return (_shape - 1.0) / _invScale; } } @@ -376,10 +358,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -387,10 +366,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -398,10 +374,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return double.PositiveInfinity; - } + get { return double.PositiveInfinity; } } /// @@ -420,12 +393,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_shape == 1.0) { return _invScale * Math.Exp(-_invScale * x); } - + return Math.Pow(_invScale, _shape) * Math.Pow(x, _shape - 1.0) * Math.Exp(-_invScale * x) / SpecialFunctions.Gamma(_shape); } @@ -440,39 +413,18 @@ namespace MathNet.Numerics.Distributions { return x == _shape ? Double.PositiveInfinity : Double.NegativeInfinity; } - + if (_shape == 0.0 && _invScale == 0.0) { return Double.NegativeInfinity; } - + if (_shape == 1.0) { return Math.Log(_invScale) - (_invScale * x); } - - return (_shape * Math.Log(_invScale)) + ((_shape - 1.0) * Math.Log(x)) - (_invScale * x) - SpecialFunctions.GammaLn(_shape); - } - - /// - /// Generates a sample from the Erlang distribution. - /// - /// a sample from the distribution. - public double Sample() - { - return DoSample(RandomSource, _shape, _invScale); - } - /// - /// Generates a sequence of samples from the Erlang distribution. - /// - /// a sequence of samples from the distribution. - public IEnumerable Samples() - { - while (true) - { - yield return DoSample(RandomSource, _shape, _invScale); - } + return (_shape * Math.Log(_invScale)) + ((_shape - 1.0) * Math.Log(x)) - (_invScale * x) - SpecialFunctions.GammaLn(_shape); } #endregion @@ -487,13 +439,13 @@ namespace MathNet.Numerics.Distributions /// The shape of the Gamma distribution. /// The inverse scale of the Gamma distribution. /// A sample from a Erlang distributed random variable. - private static double DoSample(Random rnd, double shape, double invScale) + internal static double SampleUnchecked(Random rnd, double shape, double invScale) { if (Double.IsPositiveInfinity(invScale)) { return shape; } - + var a = shape; var alphafix = 1.0; @@ -530,5 +482,63 @@ namespace MathNet.Numerics.Distributions } } } + + /// + /// Generates a sample from the Erlang distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(RandomSource, _shape, _invScale); + } + + /// + /// Generates a sequence of samples from the Erlang distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(RandomSource, _shape, _invScale); + } + } + + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The shape of the Gamma distribution. + /// The inverse scale of the Gamma distribution. + /// a sample from the distribution. + public static double Sample(Random rnd, double shape, double invScale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, shape, invScale); + } + + /// + /// Generates a sequence of samples from the distribution. + /// + /// The random number generator to use. + /// The shape of the Gamma distribution. + /// The inverse scale of the Gamma distribution. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double shape, double invScale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, shape, invScale); + } + } } } diff --git a/src/Numerics/Distributions/Continuous/Exponential.cs b/src/Numerics/Distributions/Continuous/Exponential.cs index 5855854c..0e639c65 100644 --- a/src/Numerics/Distributions/Continuous/Exponential.cs +++ b/src/Numerics/Distributions/Continuous/Exponential.cs @@ -44,12 +44,12 @@ namespace MathNet.Numerics.Distributions /// /// The lambda parameter of the Exponential distribution. /// - private double _lambda; + double _lambda; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -68,7 +68,7 @@ namespace MathNet.Numerics.Distributions /// /// Lambda parameter. /// When the parameters don't pass the function. - private void SetParameters(double lambda) + void SetParameters(double lambda) { if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) { @@ -83,7 +83,7 @@ namespace MathNet.Numerics.Distributions /// /// Lambda parameter. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double lambda) + static bool IsValidParameterSet(double lambda) { if (lambda < 0) { @@ -103,15 +103,9 @@ namespace MathNet.Numerics.Distributions /// public double Lambda { - get - { - return _lambda; - } + get { return _lambda; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -130,10 +124,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -151,10 +142,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return 1.0 / _lambda; - } + get { return 1.0 / _lambda; } } /// @@ -162,10 +150,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return 1.0 / (_lambda * _lambda); - } + get { return 1.0 / (_lambda * _lambda); } } /// @@ -173,10 +158,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return 1.0 / _lambda; - } + get { return 1.0 / _lambda; } } /// @@ -184,10 +166,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return 1.0 - Math.Log(_lambda); - } + get { return 1.0 - Math.Log(_lambda); } } /// @@ -195,10 +174,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return 2.0; - } + get { return 2.0; } } /// @@ -225,10 +201,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -236,10 +209,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return Math.Log(2.0) / _lambda; - } + get { return Math.Log(2.0) / _lambda; } } /// @@ -247,10 +217,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -258,10 +225,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -289,13 +253,32 @@ namespace MathNet.Numerics.Distributions return Math.Log(_lambda) - (_lambda * x); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The lambda parameter of the Exponential distribution. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double lambda) + { + var r = rnd.NextDouble(); + while (r == 0.0) + { + r = rnd.NextDouble(); + } + + return -Math.Log(r) / lambda; + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource, _lambda); + return SampleUnchecked(RandomSource, _lambda); } /// @@ -306,7 +289,7 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _lambda); + yield return SampleUnchecked(RandomSource, _lambda); } } @@ -322,8 +305,8 @@ namespace MathNet.Numerics.Distributions { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - - return DoSample(rnd, lambda); + + return SampleUnchecked(rnd, lambda); } /// @@ -341,26 +324,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return DoSample(rnd, lambda); - } - } - #endregion - - /// - /// Samples the distribution. - /// - /// The random number generator to use. - /// The lambda parameter of the Exponential distribution. - /// a random number from the distribution. - private static double DoSample(Random rnd, double lambda) - { - var r = rnd.NextDouble(); - while (r == 0.0) - { - r = rnd.NextDouble(); + yield return SampleUnchecked(rnd, lambda); } - - return -Math.Log(r) / lambda; } } } diff --git a/src/Numerics/Distributions/Continuous/FisherSnedecor.cs b/src/Numerics/Distributions/Continuous/FisherSnedecor.cs index 9e382890..638889ba 100644 --- a/src/Numerics/Distributions/Continuous/FisherSnedecor.cs +++ b/src/Numerics/Distributions/Continuous/FisherSnedecor.cs @@ -44,17 +44,17 @@ namespace MathNet.Numerics.Distributions /// /// The first parameter - degree of freedom. /// - private double _d1; + double _d1; /// /// The second parameter - degree of freedom. /// - private double _d2; + double _d2; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -76,7 +76,7 @@ namespace MathNet.Numerics.Distributions /// /// The first parameter - degree of freedom. /// The second parameter - degree of freedom. - private void SetParameters(double d1, double d2) + void SetParameters(double d1, double d2) { if (Control.CheckDistributionParameters && !IsValidParameterSet(d1, d2)) { @@ -93,7 +93,7 @@ namespace MathNet.Numerics.Distributions /// The first parameter - degree of freedom. /// The second parameter - degree of freedom. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double d1, double d2) + static bool IsValidParameterSet(double d1, double d2) { if (d1 <= 0 || d2 <= 0) { @@ -113,15 +113,9 @@ namespace MathNet.Numerics.Distributions /// public double DegreeOfFreedom1 { - get - { - return _d1; - } + get { return _d1; } - set - { - SetParameters(value, _d2); - } + set { SetParameters(value, _d2); } } /// @@ -129,15 +123,9 @@ namespace MathNet.Numerics.Distributions /// public double DegreeOfFreedom2 { - get - { - return _d2; - } + get { return _d2; } - set - { - SetParameters(_d1, value); - } + set { SetParameters(_d1, value); } } /// @@ -156,10 +144,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -209,10 +194,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(Variance); - } + get { return Math.Sqrt(Variance); } } /// @@ -220,10 +202,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -277,10 +256,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -288,10 +264,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -299,10 +272,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return double.PositiveInfinity; - } + get { return double.PositiveInfinity; } } /// @@ -325,13 +295,27 @@ namespace MathNet.Numerics.Distributions return Math.Log(Density(x)); } + #endregion + + /// + /// Generates one sample from the FisherSnedecor distribution without parameter checking. + /// + /// The random number generator to use. + /// The first parameter - degree of freedom. + /// The second parameter - degree of freedom. + /// a FisherSnedecor distributed random number. + internal static double SampleUnchecked(Random rnd, double d1, double d2) + { + return (ChiSquare.Sample(rnd, d1) / d1) / (ChiSquare.Sample(rnd, d2) / d2); + } + /// /// Generates a sample from the FisherSnedecor distribution. /// /// a sample from the distribution. public double Sample() { - return DoSample(RandomSource, _d1, _d2); + return SampleUnchecked(RandomSource, _d1, _d2); } /// @@ -342,22 +326,45 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _d1, _d2); + yield return SampleUnchecked(RandomSource, _d1, _d2); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The first parameter - degree of freedom. + /// The second parameter - degree of freedom. + /// a sample from the distribution. + public static double Sample(Random rnd, double d1, double d2) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(d1, d2)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, d1, d2); + } /// - /// Generates one sample from the FisherSnedecor distribution without parameter checking. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. /// The first parameter - degree of freedom. /// The second parameter - degree of freedom. - /// a FisherSnedecor distributed random number. - private static double DoSample(Random rnd, double d1, double d2) + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double d1, double d2) { - return (ChiSquare.Sample(rnd, d1) / d1) / (ChiSquare.Sample(rnd, d2) / d2); + if (Control.CheckDistributionParameters && !IsValidParameterSet(d1, d2)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, d1, d2); + } } } } diff --git a/src/Numerics/Distributions/Continuous/Gamma.cs b/src/Numerics/Distributions/Continuous/Gamma.cs index 39a11a86..755d4ffb 100644 --- a/src/Numerics/Distributions/Continuous/Gamma.cs +++ b/src/Numerics/Distributions/Continuous/Gamma.cs @@ -52,17 +52,17 @@ namespace MathNet.Numerics.Distributions /// /// Gamma shape parameter. /// - private double _shape; + double _shape; /// /// Gamma inverse scale parameter. /// - private double _invScale; + double _invScale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Gamma class. @@ -114,7 +114,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Gamma distribution. /// The inverse scale of the Gamma distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double shape, double invScale) + static bool IsValidParameterSet(double shape, double invScale) { if (shape < 0.0 || invScale < 0.0 || Double.IsNaN(shape) || Double.IsNaN(invScale)) { @@ -130,7 +130,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Gamma distribution. /// The inverse scale of the Gamma distribution. /// When the parameters don't pass the function. - private void SetParameters(double shape, double invScale) + void SetParameters(double shape, double invScale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) { @@ -146,15 +146,9 @@ namespace MathNet.Numerics.Distributions /// public double Shape { - get - { - return _shape; - } + get { return _shape; } - set - { - SetParameters(value, _invScale); - } + set { SetParameters(value, _invScale); } } /// @@ -162,10 +156,7 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return 1.0 / _invScale; - } + get { return 1.0 / _invScale; } set { @@ -185,15 +176,9 @@ namespace MathNet.Numerics.Distributions /// public double InvScale { - get - { - return _invScale; - } + get { return _invScale; } - set - { - SetParameters(_shape, value); - } + set { SetParameters(_shape, value); } } #region IDistribution implementation @@ -203,10 +188,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -230,12 +212,12 @@ namespace MathNet.Numerics.Distributions { return _shape; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return _shape / _invScale; } } @@ -251,12 +233,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return _shape / (_invScale * _invScale); } } @@ -272,12 +254,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return Math.Sqrt(_shape / (_invScale * _invScale)); } } @@ -293,12 +275,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + ((1.0 - _shape) * SpecialFunctions.DiGamma(_shape)); } } @@ -314,12 +296,12 @@ namespace MathNet.Numerics.Distributions { return 0.0; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return 2.0 / Math.Sqrt(_shape); } } @@ -339,12 +321,12 @@ namespace MathNet.Numerics.Distributions { return _shape; } - + if (_invScale == 0.0 && _shape == 0.0) { return Double.NaN; } - + return (_shape - 1.0) / _invScale; } } @@ -354,10 +336,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -365,10 +344,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -376,10 +352,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -393,17 +366,17 @@ namespace MathNet.Numerics.Distributions { return x == _shape ? Double.PositiveInfinity : 0.0; } - + if (_shape == 0.0 && _invScale == 0.0) { return 0.0; } - + if (_shape == 1.0) { return _invScale * Math.Exp(-_invScale * x); } - + return Math.Pow(_invScale, _shape) * Math.Pow(x, _shape - 1.0) * Math.Exp(-_invScale * x) / SpecialFunctions.Gamma(_shape); } @@ -418,17 +391,17 @@ namespace MathNet.Numerics.Distributions { return x == _shape ? Double.PositiveInfinity : Double.NegativeInfinity; } - + if (_shape == 0.0 && _invScale == 0.0) { return Double.NegativeInfinity; } - + if (_shape == 1.0) { return Math.Log(_invScale) - (_invScale * x); } - + return (_shape * Math.Log(_invScale)) + ((_shape - 1.0) * Math.Log(x)) - (_invScale * x) - SpecialFunctions.GammaLn(_shape); } @@ -443,75 +416,17 @@ namespace MathNet.Numerics.Distributions { return x >= _shape ? 1.0 : 0.0; } - + if (_shape == 0.0 && _invScale == 0.0) { return 0.0; } - - return SpecialFunctions.GammaLowerRegularized(_shape, x * _invScale); - } - /// - /// Generates a sample from the Gamma distribution. - /// - /// a sample from the distribution. - public double Sample() - { - return SampleGamma(RandomSource, _shape, _invScale); - } - - /// - /// Generates a sequence of samples from the Gamma distribution. - /// - /// a sequence of samples from the distribution. - public IEnumerable Samples() - { - while (true) - { - yield return SampleGamma(RandomSource, _shape, _invScale); - } + return SpecialFunctions.GammaLowerRegularized(_shape, x * _invScale); } #endregion - /// - /// Generates a sample from the Gamma distribution. - /// - /// The random number generator to use. - /// The shape of the Gamma distribution from which to generate samples. - /// The inverse scale of the Gamma distribution from which to generate samples. - /// a sample from the distribution. - public static double Sample(Random rng, double shape, double invScale) - { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } - - return SampleGamma(rng, shape, invScale); - } - - /// - /// Generates a sequence of samples from the Gamma distribution. - /// - /// The random number generator to use. - /// The shape of the Gamma distribution from which to generate samples. - /// The inverse scale of the Gamma distribution from which to generate samples. - /// a sequence of samples from the distribution. - public static IEnumerable Samples(Random rng, double shape, double invScale) - { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } - - while (true) - { - yield return SampleGamma(rng, shape, invScale); - } - } - /// /// Sampling implementation based on: /// "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang @@ -522,7 +437,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Gamma distribution. /// The inverse scale of the Gamma distribution. /// A sample from a Gamma distributed random variable. - internal static double SampleGamma(Random rnd, double shape, double invScale) + internal static double SampleUnchecked(Random rnd, double shape, double invScale) { if (Double.IsPositiveInfinity(invScale)) { @@ -565,5 +480,63 @@ namespace MathNet.Numerics.Distributions } } } + + /// + /// Generates a sample from the Gamma distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleUnchecked(RandomSource, _shape, _invScale); + } + + /// + /// Generates a sequence of samples from the Gamma distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return SampleUnchecked(RandomSource, _shape, _invScale); + } + } + + /// + /// Generates a sample from the Gamma distribution. + /// + /// The random number generator to use. + /// The shape of the Gamma distribution from which to generate samples. + /// The inverse scale of the Gamma distribution from which to generate samples. + /// a sample from the distribution. + public static double Sample(Random rng, double shape, double invScale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rng, shape, invScale); + } + + /// + /// Generates a sequence of samples from the Gamma distribution. + /// + /// The random number generator to use. + /// The shape of the Gamma distribution from which to generate samples. + /// The inverse scale of the Gamma distribution from which to generate samples. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rng, double shape, double invScale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rng, shape, invScale); + } + } } } diff --git a/src/Numerics/Distributions/Continuous/InverseGamma.cs b/src/Numerics/Distributions/Continuous/InverseGamma.cs index 67a56f26..834b715f 100644 --- a/src/Numerics/Distributions/Continuous/InverseGamma.cs +++ b/src/Numerics/Distributions/Continuous/InverseGamma.cs @@ -45,17 +45,17 @@ namespace MathNet.Numerics.Distributions /// /// Inverse Gamma shape parameter. /// - private double _shape; + double _shape; /// /// Inverse Gamma scale parameter scale. /// - private double _scale; + double _scale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -82,7 +82,7 @@ namespace MathNet.Numerics.Distributions /// The scale (beta) parameter of the inverse Gamma distribution. /// /// When the parameters don't pass the function. - private void SetParameters(double shape, double scale) + void SetParameters(double shape, double scale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) { @@ -103,7 +103,7 @@ namespace MathNet.Numerics.Distributions /// The scale (beta) parameter of the inverse Gamma distribution. /// /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double shape, double scale) + static bool IsValidParameterSet(double shape, double scale) { if (shape <= 0 || scale <= 0) { @@ -123,15 +123,9 @@ namespace MathNet.Numerics.Distributions /// public double Shape { - get - { - return _shape; - } + get { return _shape; } - set - { - SetParameters(value, _scale); - } + set { SetParameters(value, _scale); } } /// @@ -139,15 +133,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(_shape, value); - } + set { SetParameters(_shape, value); } } /// @@ -166,10 +154,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -219,10 +204,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return _scale / (Math.Abs(_shape - 1.0) * Math.Sqrt(_shape - 2.0)); - } + get { return _scale / (Math.Abs(_shape - 1.0) * Math.Sqrt(_shape - 2.0)); } } /// @@ -230,10 +212,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return _shape + Math.Log(_scale) + SpecialFunctions.GammaLn(_shape) - ((1 + _shape) * SpecialFunctions.DiGamma(_shape)); - } + get { return _shape + Math.Log(_scale) + SpecialFunctions.GammaLn(_shape) - ((1 + _shape) * SpecialFunctions.DiGamma(_shape)); } } /// @@ -271,10 +250,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return _scale / (_shape + 1.0); - } + get { return _scale / (_shape + 1.0); } } /// @@ -283,10 +259,7 @@ namespace MathNet.Numerics.Distributions /// Throws . public double Median { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -294,10 +267,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -305,10 +275,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -336,13 +303,27 @@ namespace MathNet.Numerics.Distributions return Math.Log(Density(x)); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The shape (alpha) parameter of the inverse Gamma distribution. + /// The scale (beta) parameter of the inverse Gamma distribution. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double shape, double scale) + { + return 1.0 / Gamma.Sample(rnd, shape, scale); + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource, _shape, _scale); + return SampleUnchecked(RandomSource, _shape, _scale); } /// @@ -353,26 +334,45 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _shape, _scale); + yield return SampleUnchecked(RandomSource, _shape, _scale); } } - #endregion - /// - /// Samples the distribution. + /// Generates a sample from the distribution. /// /// The random number generator to use. - /// - /// The shape (alpha) parameter of the inverse Gamma distribution. - /// - /// - /// The scale (beta) parameter of the inverse Gamma distribution. - /// - /// a random number from the distribution. - private static double DoSample(Random rnd, double shape, double scale) + /// The shape (alpha) parameter of the inverse Gamma distribution. + /// The scale (beta) parameter of the inverse Gamma distribution. + /// a sample from the distribution. + public static double Sample(Random rnd, double shape, double scale) { - return 1.0 / Gamma.Sample(rnd, shape, scale); + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, shape, scale); + } + + /// + /// Generates a sequence of samples from the distribution. + /// + /// The random number generator to use. + /// The shape (alpha) parameter of the inverse Gamma distribution. + /// The scale (beta) parameter of the inverse Gamma distribution. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double shape, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, shape, scale); + } } } } diff --git a/src/Numerics/Distributions/Continuous/Laplace.cs b/src/Numerics/Distributions/Continuous/Laplace.cs index f822fbfe..a9a41e13 100644 --- a/src/Numerics/Distributions/Continuous/Laplace.cs +++ b/src/Numerics/Distributions/Continuous/Laplace.cs @@ -46,27 +46,21 @@ namespace MathNet.Numerics.Distributions /// /// The scale of the Laplace distribution. /// - private double _scale; + double _scale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Gets or sets the location of the Laplace distribution. /// public double Location { - get - { - return Mean; - } + get { return Mean; } - set - { - SetParameters(value, _scale); - } + set { SetParameters(value, _scale); } } /// @@ -74,21 +68,16 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(Mean, value); - } + set { SetParameters(Mean, value); } } /// /// Initializes a new instance of the class (location = 0, scale = 1). /// - public Laplace() : this(0.0, 1.0) + public Laplace() + : this(0.0, 1.0) { } @@ -116,7 +105,7 @@ namespace MathNet.Numerics.Distributions /// The location for the Laplace distribution. /// The scale for the Laplace distribution. /// When the parameters don't pass the function. - private void SetParameters(double location, double scale) + void SetParameters(double location, double scale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) { @@ -133,7 +122,7 @@ namespace MathNet.Numerics.Distributions /// The location for the Laplace distribution. /// The scale for the Laplace distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double location, double scale) + static bool IsValidParameterSet(double location, double scale) { if (scale <= 0) { @@ -164,10 +153,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -183,21 +169,14 @@ namespace MathNet.Numerics.Distributions /// /// Gets the mean of the distribution. /// - public double Mean - { - get; - private set; - } + public double Mean { get; private set; } /// /// Gets the variance of the distribution. /// public double Variance { - get - { - return 2.0 * _scale * _scale; - } + get { return 2.0 * _scale * _scale; } } /// @@ -205,10 +184,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(2.0) * _scale; - } + get { return Math.Sqrt(2.0) * _scale; } } /// @@ -216,10 +192,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return Math.Log(2.0 * Constants.E * _scale); - } + get { return Math.Log(2.0 * Constants.E * _scale); } } /// @@ -227,10 +200,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -252,10 +222,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return Mean; - } + get { return Mean; } } /// @@ -263,10 +230,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return Mean; - } + get { return Mean; } } /// @@ -274,10 +238,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return Double.NegativeInfinity; - } + get { return Double.NegativeInfinity; } } /// @@ -285,10 +246,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -311,13 +269,28 @@ namespace MathNet.Numerics.Distributions return Math.Log(Density(x)); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The location shape parameter. + /// The scale parameter. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double location, double scale) + { + var u = rnd.NextDouble() - 0.5; + return location - (scale * Math.Sign(u) * Math.Log(1.0 - (2.0 * Math.Abs(u)))); + } + /// /// Samples a Laplace distributed random variable. /// /// a sample from the distribution. public double Sample() { - return DoSample(RandomSource, Mean, _scale); + return SampleUnchecked(RandomSource, Mean, _scale); } /// @@ -328,23 +301,45 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, Mean, _scale); + yield return SampleUnchecked(RandomSource, Mean, _scale); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The location shape parameter. + /// The scale parameter. + /// a sample from the distribution. + public static double Sample(Random rnd, double location, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, location, scale); + } /// - /// Samples the distribution. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. /// The location shape parameter. /// The scale parameter. - /// a random number from the distribution. - private static double DoSample(Random rnd, double location, double scale) + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double location, double scale) { - var u = rnd.NextDouble() - 0.5; - return location - (scale * Math.Sign(u) * Math.Log(1.0 - (2.0 * Math.Abs(u)))); + if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, location, scale); + } } } } diff --git a/src/Numerics/Distributions/Continuous/LogNormal.cs b/src/Numerics/Distributions/Continuous/LogNormal.cs index 1c723281..2e4f2b7e 100644 --- a/src/Numerics/Distributions/Continuous/LogNormal.cs +++ b/src/Numerics/Distributions/Continuous/LogNormal.cs @@ -44,17 +44,17 @@ namespace MathNet.Numerics.Distributions /// /// Keeps track of the mu of the logarithm of the log-log-normal distribution. /// - private double _mu; + double _mu; /// /// Keeps track of the standard deviation of the logarithm of the log-log-normal distribution. /// - private double _sigma; + double _sigma; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -88,7 +88,7 @@ namespace MathNet.Numerics.Distributions /// The mu of the logarithm of the distribution. /// The standard deviation of the logarithm of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double mu, double sigma) + static bool IsValidParameterSet(double mu, double sigma) { if (sigma < 0.0 || Double.IsNaN(mu) || Double.IsNaN(mu) || Double.IsNaN(sigma)) { @@ -104,7 +104,7 @@ namespace MathNet.Numerics.Distributions /// The mu of the logarithm of the distribution. /// The standard deviation of the logarithm of the distribution. /// When the parameters don't pass the function. - private void SetParameters(double mu, double sigma) + void SetParameters(double mu, double sigma) { if (Control.CheckDistributionParameters && !IsValidParameterSet(mu, sigma)) { @@ -120,15 +120,9 @@ namespace MathNet.Numerics.Distributions /// public double Mu { - get - { - return _mu; - } + get { return _mu; } - set - { - SetParameters(value, _sigma); - } + set { SetParameters(value, _sigma); } } /// @@ -136,15 +130,9 @@ namespace MathNet.Numerics.Distributions /// public double Sigma { - get - { - return _sigma; - } + get { return _sigma; } - set - { - SetParameters(_mu, value); - } + set { SetParameters(_mu, value); } } #region IDistribution implementation @@ -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 Math.Exp(_mu + (_sigma * _sigma / 2.0)); - } + get { return Math.Exp(_mu + (_sigma * _sigma / 2.0)); } } /// @@ -210,10 +192,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return 0.5 + Math.Log(_sigma) + _mu + Constants.LogSqrt2Pi; - } + get { return 0.5 + Math.Log(_sigma) + _mu + Constants.LogSqrt2Pi; } } /// @@ -237,10 +216,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return Math.Exp(_mu - (_sigma * _sigma)); - } + get { return Math.Exp(_mu - (_sigma * _sigma)); } } /// @@ -248,10 +224,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return Math.Exp(_mu); - } + get { return Math.Exp(_mu); } } /// @@ -259,10 +232,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -270,10 +240,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -323,13 +290,15 @@ namespace MathNet.Numerics.Distributions return 0.5 * (1.0 + SpecialFunctions.Erf((Math.Log(x) - _mu) / (_sigma * Constants.Sqrt2))); } + #endregion + /// /// Generates a sample from the log-normal distribution using the Box-Muller algorithm. /// /// a sample from the distribution. public double Sample() { - return Math.Exp(_mu + (_sigma * Normal.SampleBoxMuller(RandomSource).Item1)); + return Math.Exp(Normal.SampleUnchecked(RandomSource, _mu, _sigma)); } /// @@ -340,14 +309,12 @@ namespace MathNet.Numerics.Distributions { while (true) { - var sample = Normal.SampleBoxMuller(RandomSource); + var sample = Normal.SampleUncheckedBoxMuller(RandomSource); yield return Math.Exp(_mu + (_sigma * sample.Item1)); yield return Math.Exp(_mu + (_sigma * sample.Item2)); } } - #endregion - /// /// Generates a sample from the log-normal distribution using the Box-Muller algorithm. /// @@ -362,7 +329,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return Math.Exp(mu + (sigma * Normal.SampleBoxMuller(rng).Item1)); + return Math.Exp(Normal.SampleUnchecked(rng, mu, sigma)); } /// @@ -381,7 +348,7 @@ namespace MathNet.Numerics.Distributions while (true) { - var sample = Normal.SampleBoxMuller(rng); + var sample = Normal.SampleUncheckedBoxMuller(rng); yield return Math.Exp(mu + (sigma * sample.Item1)); yield return Math.Exp(mu + (sigma * sample.Item2)); } diff --git a/src/Numerics/Distributions/Continuous/Normal.cs b/src/Numerics/Distributions/Continuous/Normal.cs index f497e3f5..2658d4e4 100644 --- a/src/Numerics/Distributions/Continuous/Normal.cs +++ b/src/Numerics/Distributions/Continuous/Normal.cs @@ -44,24 +44,25 @@ namespace MathNet.Numerics.Distributions /// /// Keeps track of the mean of the normal distribution. /// - private double _mean; + double _mean; /// /// Keeps track of the standard deviation of the normal distribution. /// - private double _stdDev; + double _stdDev; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Normal class. This is a normal distribution with mean 0.0 /// and standard deviation 1.0. The distribution will /// be initialized with the default random number generator. /// - public Normal() : this(0.0, 1.0) + public Normal() + : this(0.0, 1.0) { } @@ -128,7 +129,7 @@ namespace MathNet.Numerics.Distributions /// The mean of the normal distribution. /// The standard deviation of the normal distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double mean, double stddev) + static bool IsValidParameterSet(double mean, double stddev) { if (stddev < 0.0 || Double.IsNaN(mean) || Double.IsNaN(stddev)) { @@ -144,7 +145,7 @@ namespace MathNet.Numerics.Distributions /// The mean of the normal distribution. /// The standard deviation of the normal distribution. /// When the parameters don't pass the function. - private void SetParameters(double mean, double stddev) + void SetParameters(double mean, double stddev) { if (Control.CheckDistributionParameters && !IsValidParameterSet(mean, stddev)) { @@ -160,10 +161,7 @@ namespace MathNet.Numerics.Distributions /// public double Precision { - get - { - return 1.0 / (_stdDev * _stdDev); - } + get { return 1.0 / (_stdDev * _stdDev); } set { @@ -186,10 +184,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -207,15 +202,9 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _mean; - } + get { return _mean; } - set - { - SetParameters(value, _stdDev); - } + set { SetParameters(value, _stdDev); } } /// @@ -223,15 +212,9 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return _stdDev * _stdDev; - } + get { return _stdDev * _stdDev; } - set - { - SetParameters(_mean, Math.Sqrt(value)); - } + set { SetParameters(_mean, Math.Sqrt(value)); } } /// @@ -239,15 +222,9 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return _stdDev; - } + get { return _stdDev; } - set - { - SetParameters(_mean, value); - } + set { SetParameters(_mean, value); } } /// @@ -255,10 +232,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return Math.Log(_stdDev) + Constants.LogSqrt2PiE; - } + get { return Math.Log(_stdDev) + Constants.LogSqrt2PiE; } } /// @@ -266,10 +240,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return 0.0; - } + get { return 0.0; } } #endregion @@ -281,10 +252,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return _mean; - } + get { return _mean; } } /// @@ -292,10 +260,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return _mean; - } + get { return _mean; } } /// @@ -303,10 +268,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return Double.NegativeInfinity; - } + get { return Double.NegativeInfinity; } } /// @@ -314,10 +276,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -388,13 +347,60 @@ namespace MathNet.Numerics.Distributions return CumulativeDistribution(_mean, _stdDev, x); } + #endregion + + /// + /// Computes the inverse cumulative distribution function of the normal distribution. + /// + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + public double InverseCumulativeDistribution(double p) + { + return _mean - (_stdDev * Math.Sqrt(2.0) * SpecialFunctions.ErfcInv(2.0 * p)); + } + + + + /// + /// Samples a pair of standard normal distributed random variables using the Box-Muller algorithm. + /// + /// The random number generator to use. + /// a pair of random numbers from the standard normal distribution. + internal static Tuple SampleUncheckedBoxMuller(Random rnd) + { + var v1 = (2.0 * rnd.NextDouble()) - 1.0; + var v2 = (2.0 * rnd.NextDouble()) - 1.0; + var r = (v1 * v1) + (v2 * v2); + while (r >= 1.0 || r == 0.0) + { + v1 = (2.0 * rnd.NextDouble()) - 1.0; + v2 = (2.0 * rnd.NextDouble()) - 1.0; + r = (v1 * v1) + (v2 * v2); + } + + var fac = Math.Sqrt(-2.0 * Math.Log(r) / r); + return new Tuple(v1 * fac, v2 * fac); + } + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The mean of the normal distribution from which to generate samples. + /// The standard deviation of the normal distribution from which to generate samples. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double mean, double stddev) + { + return mean + (stddev * SampleUncheckedBoxMuller(rnd).Item1); + } + /// /// Generates a sample from the normal distribution using the Box-Muller algorithm. /// /// a sample from the distribution. public double Sample() { - return _mean + (_stdDev * SampleBoxMuller(RandomSource).Item1); + return SampleUnchecked(RandomSource, _mean, _stdDev); } /// @@ -405,49 +411,37 @@ namespace MathNet.Numerics.Distributions { while (true) { - var sample = SampleBoxMuller(RandomSource); + var sample = SampleUncheckedBoxMuller(RandomSource); yield return _mean + (_stdDev * sample.Item1); yield return _mean + (_stdDev * sample.Item2); } } - #endregion - - /// - /// Computes the inverse cumulative distribution function of the normal distribution. - /// - /// The location at which to compute the inverse cumulative density. - /// the inverse cumulative density at . - public double InverseCumulativeDistribution(double p) - { - return _mean - (_stdDev * Math.Sqrt(2.0) * SpecialFunctions.ErfcInv(2.0 * p)); - } - /// /// Generates a sample from the normal distribution using the Box-Muller algorithm. /// - /// The random number generator to use. + /// The random number generator to use. /// The mean of the normal distribution from which to generate samples. /// The standard deviation of the normal distribution from which to generate samples. /// a sample from the distribution. - public static double Sample(Random rng, double mean, double stddev) + public static double Sample(Random rnd, double mean, double stddev) { if (Control.CheckDistributionParameters && !IsValidParameterSet(mean, stddev)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return mean + (stddev * SampleBoxMuller(rng).Item1); + return SampleUnchecked(rnd, mean, stddev); } /// /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. /// - /// The random number generator to use. + /// The random number generator to use. /// The mean of the normal distribution from which to generate samples. /// The standard deviation of the normal distribution from which to generate samples. /// a sequence of samples from the distribution. - public static IEnumerable Samples(Random rng, double mean, double stddev) + public static IEnumerable Samples(Random rnd, double mean, double stddev) { if (Control.CheckDistributionParameters && !IsValidParameterSet(mean, stddev)) { @@ -456,31 +450,10 @@ namespace MathNet.Numerics.Distributions while (true) { - var sample = SampleBoxMuller(rng); + var sample = SampleUncheckedBoxMuller(rnd); yield return mean + (stddev * sample.Item1); yield return mean + (stddev * sample.Item2); } } - - /// - /// Samples a pair of standard normal distributed random variables using the Box-Muller algorithm. - /// - /// The random number generator to use. - /// a pair of random numbers from the standard normal distribution. - internal static Tuple SampleBoxMuller(Random rnd) - { - var v1 = (2.0 * rnd.NextDouble()) - 1.0; - var v2 = (2.0 * rnd.NextDouble()) - 1.0; - var r = (v1 * v1) + (v2 * v2); - while (r >= 1.0 || r == 0.0) - { - v1 = (2.0 * rnd.NextDouble()) - 1.0; - v2 = (2.0 * rnd.NextDouble()) - 1.0; - r = (v1 * v1) + (v2 * v2); - } - - var fac = Math.Sqrt(-2.0 * Math.Log(r) / r); - return new Tuple(v1 * fac, v2 * fac); - } } } diff --git a/src/Numerics/Distributions/Continuous/Pareto.cs b/src/Numerics/Distributions/Continuous/Pareto.cs index 97c3bc34..276a20a5 100644 --- a/src/Numerics/Distributions/Continuous/Pareto.cs +++ b/src/Numerics/Distributions/Continuous/Pareto.cs @@ -46,17 +46,17 @@ namespace MathNet.Numerics.Distributions /// /// The scale parameter of the distribution. /// - private double _scale; + double _scale; /// /// The shape parameter of the distribution. /// - private double _shape; + double _shape; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -82,7 +82,7 @@ namespace MathNet.Numerics.Distributions /// The scale parameter of the distribution. /// The shape parameter of the distribution. /// When the parameters don't pass the function. - private void SetParameters(double scale, double shape) + void SetParameters(double scale, double shape) { if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) { @@ -99,7 +99,7 @@ namespace MathNet.Numerics.Distributions /// The scale parameter of the distribution. /// The shape parameter of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double scale, double shape) + static bool IsValidParameterSet(double scale, double shape) { if (scale <= 0 || shape <= 0) { @@ -119,15 +119,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(value, _shape); - } + set { SetParameters(value, _shape); } } /// @@ -135,15 +129,9 @@ namespace MathNet.Numerics.Distributions /// public double Shape { - get - { - return _shape; - } + get { return _shape; } - set - { - SetParameters(_scale, value); - } + set { SetParameters(_scale, value); } } /// @@ -162,10 +150,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -215,10 +200,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return (_scale * Math.Sqrt(_shape)) / (Math.Abs(_shape - 1.0) * Math.Sqrt(_shape - 2.0)); - } + get { return (_scale * Math.Sqrt(_shape)) / (Math.Abs(_shape - 1.0) * Math.Sqrt(_shape - 2.0)); } } /// @@ -226,10 +208,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return Math.Log(_shape / _scale) - (1.0 / _shape) - 1.0; - } + get { return Math.Log(_shape / _scale) - (1.0 / _shape) - 1.0; } } /// @@ -237,10 +216,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (2.0 * (_shape + 1.0) / (_shape - 3.0)) * Math.Sqrt((_shape - 2.0) / _shape); - } + get { return (2.0 * (_shape + 1.0) / (_shape - 3.0)) * Math.Sqrt((_shape - 2.0) / _shape); } } /// @@ -262,10 +238,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return _scale; - } + get { return _scale; } } /// @@ -273,10 +246,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return _scale * Math.Pow(2.0, 1.0 / _shape); - } + get { return _scale * Math.Pow(2.0, 1.0 / _shape); } } /// @@ -284,10 +254,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return _scale; - } + get { return _scale; } } /// @@ -295,10 +262,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -321,13 +285,27 @@ namespace MathNet.Numerics.Distributions return Math.Log(Density(x)); } + #endregion + + /// + /// Generates a sample from the Pareto distribution without doing parameter checking. + /// + /// The random number generator to use. + /// The scale parameter. + /// The shape parameter. + /// a random number from the Pareto distribution. + internal static double SampleUnchecked(Random rnd, double scale, double shape) + { + return scale * Math.Pow(rnd.NextDouble(), -1.0 / shape); + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource); + return SampleUnchecked(RandomSource, _scale, _shape); } /// @@ -338,20 +316,45 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource); + yield return SampleUnchecked(RandomSource, _scale, _shape); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The scale parameter. + /// The shape parameter. + /// a sample from the distribution. + public static double Sample(Random rnd, double scale, double shape) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, scale, shape); + } /// - /// Generates a sample from the Pareto distribution without doing parameter checking. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. - /// a random number from the Pareto distribution. - private double DoSample(Random rnd) + /// The scale parameter. + /// The shape parameter. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double scale, double shape) { - return _scale * Math.Pow(rnd.NextDouble(), -1.0 / _shape); + if (Control.CheckDistributionParameters && !IsValidParameterSet(scale, shape)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, scale, shape); + } } } } diff --git a/src/Numerics/Distributions/Continuous/Rayleigh.cs b/src/Numerics/Distributions/Continuous/Rayleigh.cs index 6713305a..b5fcfad6 100644 --- a/src/Numerics/Distributions/Continuous/Rayleigh.cs +++ b/src/Numerics/Distributions/Continuous/Rayleigh.cs @@ -47,12 +47,12 @@ namespace MathNet.Numerics.Distributions /// /// The scale parameter of the distribution. /// - private double _scale; + double _scale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -74,7 +74,7 @@ namespace MathNet.Numerics.Distributions /// /// The scale parameter of the distribution. /// When the parameters don't pass the function. - private void SetParameters(double scale) + void SetParameters(double scale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(scale)) { @@ -89,7 +89,7 @@ namespace MathNet.Numerics.Distributions /// /// The scale parameter of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double scale) + static bool IsValidParameterSet(double scale) { if (scale <= 0) { @@ -109,15 +109,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(value); - } + set { SetParameters(value); } } /// @@ -136,10 +130,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -157,10 +148,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _scale * Math.Sqrt(Constants.PiOver2); - } + get { return _scale * Math.Sqrt(Constants.PiOver2); } } /// @@ -168,10 +156,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return (2.0 - Constants.PiOver2) * _scale * _scale; - } + get { return (2.0 - Constants.PiOver2) * _scale * _scale; } } /// @@ -179,10 +164,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(2.0 - Constants.PiOver2) * _scale; - } + get { return Math.Sqrt(2.0 - Constants.PiOver2) * _scale; } } /// @@ -190,10 +172,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return 1.0 + Math.Log(_scale / Math.Sqrt(2)) + (Constants.EulerMascheroni / 2.0); - } + get { return 1.0 + Math.Log(_scale / Math.Sqrt(2)) + (Constants.EulerMascheroni / 2.0); } } /// @@ -201,10 +180,7 @@ namespace MathNet.Numerics.Distributions /// public double Skewness { - get - { - return (2.0 * Math.Sqrt(Constants.Pi) * (Constants.Pi - 3.0)) / Math.Pow(4.0 - Constants.Pi, 1.5); - } + get { return (2.0 * Math.Sqrt(Constants.Pi) * (Constants.Pi - 3.0)) / Math.Pow(4.0 - Constants.Pi, 1.5); } } /// @@ -226,10 +202,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return _scale; - } + get { return _scale; } } /// @@ -237,10 +210,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return _scale * Math.Sqrt(Math.Log(4.0)); - } + get { return _scale * Math.Sqrt(Math.Log(4.0)); } } /// @@ -248,10 +218,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return 0.0; - } + get { return 0.0; } } /// @@ -259,10 +226,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -285,13 +249,26 @@ namespace MathNet.Numerics.Distributions return Math.Log(x / (_scale * _scale)) - (x * x / (2.0 * _scale * _scale)); } + #endregion + + /// + /// Generates a sample from the Rayleigh distribution without doing parameter checking. + /// + /// The random number generator to use. + /// The scale parameter. + /// a random number from the Rayleigh distribution. + internal static double SampleUnchecked(Random rnd, double scale) + { + return scale * Math.Sqrt(-2.0 * Math.Log(rnd.NextDouble())); + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource); + return SampleUnchecked(RandomSource, _scale); } /// @@ -302,20 +279,43 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource); + yield return SampleUnchecked(RandomSource, _scale); } } - #endregion + /// + /// Generates a sample from the distribution. + /// + /// The random number generator to use. + /// The scale parameter. + /// a sample from the distribution. + public static double Sample(Random rnd, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleUnchecked(rnd, scale); + } /// - /// Generates a sample from the Rayleigh distribution without doing parameter checking. + /// Generates a sequence of samples from the distribution. /// /// The random number generator to use. - /// a random number from the Rayleigh distribution. - private double DoSample(Random rnd) + /// The scale parameter. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double scale) { - return _scale * Math.Sqrt(-2.0 * Math.Log(rnd.NextDouble())); + if (Control.CheckDistributionParameters && !IsValidParameterSet(scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleUnchecked(rnd, scale); + } } } } diff --git a/src/Numerics/Distributions/Continuous/Stable.cs b/src/Numerics/Distributions/Continuous/Stable.cs index 7d8cfeb2..ecdf2f0d 100644 --- a/src/Numerics/Distributions/Continuous/Stable.cs +++ b/src/Numerics/Distributions/Continuous/Stable.cs @@ -47,27 +47,27 @@ namespace MathNet.Numerics.Distributions /// /// The stability parameter of the distribution. /// - private double _alpha; + double _alpha; /// /// The skewness parameter of the distribution. /// - private double _beta; + double _beta; /// /// The scale parameter of the distribution. /// - private double _scale; + double _scale; /// /// The location parameter of the distribution. /// - private double _location; + double _location; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the class. @@ -97,7 +97,7 @@ namespace MathNet.Numerics.Distributions /// The skewness parameter of the distribution. /// The scale parameter of the distribution. /// The location parameter of the distribution. - private void SetParameters(double alpha, double beta, double scale, double location) + void SetParameters(double alpha, double beta, double scale, double location) { if (Control.CheckDistributionParameters && !IsValidParameterSet(alpha, beta, scale, location)) { @@ -118,7 +118,7 @@ namespace MathNet.Numerics.Distributions /// The scale parameter of the distribution. /// The location parameter of the distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double alpha, double beta, double scale, double location) + static bool IsValidParameterSet(double alpha, double beta, double scale, double location) { if (alpha <= 0 || alpha > 2) { @@ -148,15 +148,9 @@ namespace MathNet.Numerics.Distributions /// public double Alpha { - get - { - return _alpha; - } + get { return _alpha; } - set - { - SetParameters(value, _beta, _scale, _location); - } + set { SetParameters(value, _beta, _scale, _location); } } /// @@ -164,15 +158,9 @@ namespace MathNet.Numerics.Distributions /// public double Beta { - get - { - return _beta; - } + get { return _beta; } - set - { - SetParameters(_alpha, value, _scale, _location); - } + set { SetParameters(_alpha, value, _scale, _location); } } /// @@ -180,15 +168,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(_alpha, _beta, value, _location); - } + set { SetParameters(_alpha, _beta, value, _location); } } /// @@ -196,15 +178,9 @@ namespace MathNet.Numerics.Distributions /// public double Location { - get - { - return _location; - } + get { return _location; } - set - { - SetParameters(_alpha, _beta, _scale, value); - } + set { SetParameters(_alpha, _beta, _scale, value); } } /// @@ -223,10 +199,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -293,10 +266,7 @@ namespace MathNet.Numerics.Distributions /// Always throws a not supported exception. public double Entropy { - get - { - throw new NotSupportedException(); - } + get { throw new NotSupportedException(); } } /// @@ -351,7 +321,7 @@ namespace MathNet.Numerics.Distributions /// /// the cumulative density at . /// - private static double LevyCumulativeDistribution(double scale, double location, double x) + static double LevyCumulativeDistribution(double scale, double location, double x) { // The parameters scale and location must be correct return SpecialFunctions.Erfc(Math.Sqrt(scale / (2 * (x - location)))); @@ -416,10 +386,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -454,7 +421,7 @@ namespace MathNet.Numerics.Distributions /// The location parameter of the distribution. /// The location at which to compute the density. /// the density at . - private static double LevyDensity(double scale, double location, double x) + static double LevyDensity(double scale, double location, double x) { // The parameters scale and location must be correct if (x < location) @@ -475,13 +442,51 @@ namespace MathNet.Numerics.Distributions return Math.Log(Density(x)); } + #endregion + + /// + /// Samples the distribution. + /// + /// The random number generator to use. + /// The stability parameter of the distribution. + /// The skewness parameter of the distribution. + /// The scale parameter of the distribution. + /// The location parameter of the distribution. + /// a random number from the distribution. + internal static double SampleUnchecked(Random rnd, double alpha, double beta, double scale, double location) + { + var randTheta = ContinuousUniform.Sample(rnd, -Constants.PiOver2, Constants.PiOver2); + var randW = Exponential.Sample(rnd, 1.0); + + if (!1.0.AlmostEqual(alpha)) + { + var theta = (1.0 / alpha) * Math.Atan(beta * Math.Tan(Constants.PiOver2 * alpha)); + var angle = alpha * (randTheta + theta); + var part1 = beta * Math.Tan(Constants.PiOver2 * alpha); + + var factor = Math.Pow(1.0 + (part1 * part1), 1.0 / (2.0 * alpha)); + var factor1 = Math.Sin(angle) / Math.Pow(Math.Cos(randTheta), (1.0 / alpha)); + var factor2 = Math.Pow(Math.Cos(randTheta - angle) / randW, (1 - alpha) / alpha); + + return location + scale * (factor * factor1 * factor2); + } + else + { + var part1 = Constants.PiOver2 + (beta * randTheta); + var summand = part1 * Math.Tan(randTheta); + var subtrahend = beta * Math.Log(Constants.PiOver2 * randW * Math.Cos(randTheta) / part1); + + return location + scale * ((2.0 / Math.PI) * (summand - subtrahend)); + } + } + /// /// Draws a random sample from the distribution. /// /// A random number from this distribution. public double Sample() { - return DoSample(RandomSource, _alpha, _beta); + return SampleUnchecked(RandomSource, _alpha, _beta, _scale, _location); } /// @@ -492,43 +497,50 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return DoSample(RandomSource, _alpha, _beta); + yield return SampleUnchecked(RandomSource, _alpha, _beta, _scale, _location); } } + /// - /// Samples the distribution. + /// Generates a sample from the distribution. /// /// The random number generator to use. /// The stability parameter of the distribution. /// The skewness parameter of the distribution. - /// a random number from the distribution. - private static double DoSample(Random rnd, double alpha, double beta) + /// The scale parameter of the distribution. + /// The location parameter of the distribution. + /// a sample from the distribution. + public static double Sample(Random rnd, double alpha, double beta, double scale, double location) { - var randTheta = ContinuousUniform.Sample(rnd, -Constants.PiOver2, Constants.PiOver2); - var randW = Exponential.Sample(rnd, 1.0); - - if (!1.0.AlmostEqual(alpha)) + if (Control.CheckDistributionParameters && !IsValidParameterSet(alpha, beta, scale, location)) { - var theta = (1.0 / alpha) * Math.Atan(beta * Math.Tan(Constants.PiOver2 * alpha)); - var angle = alpha * (randTheta + theta); - var part1 = beta * Math.Tan(Constants.PiOver2 * alpha); + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } - var factor = Math.Pow(1.0 + (part1 * part1), 1.0 / (2.0 * alpha)); - var factor1 = Math.Sin(angle) / Math.Pow(Math.Cos(randTheta), (1.0 / alpha)); - var factor2 = Math.Pow(Math.Cos(randTheta - angle) / randW, (1 - alpha) / alpha); + return SampleUnchecked(rnd, location, scale, scale, location); + } - return factor * factor1 * factor2; - } - else + /// + /// Generates a sequence of samples from the distribution. + /// + /// The random number generator to use. + /// The stability parameter of the distribution. + /// The skewness parameter of the distribution. + /// The scale parameter of the distribution. + /// The location parameter of the distribution. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rnd, double alpha, double beta, double scale, double location) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, scale, location)) { - var part1 = Constants.PiOver2 + (beta * randTheta); - var summand = part1 * Math.Tan(randTheta); - var subtrahend = beta * Math.Log(Constants.PiOver2 * randW * Math.Cos(randTheta) / part1); + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } - return (2.0 / Math.PI) * (summand - subtrahend); + while (true) + { + yield return SampleUnchecked(rnd, location, scale, scale, location); } } - #endregion } } diff --git a/src/Numerics/Distributions/Continuous/StudentT.cs b/src/Numerics/Distributions/Continuous/StudentT.cs index 53113176..4512e1f3 100644 --- a/src/Numerics/Distributions/Continuous/StudentT.cs +++ b/src/Numerics/Distributions/Continuous/StudentT.cs @@ -55,29 +55,30 @@ namespace MathNet.Numerics.Distributions /// /// Keeps track of the location of the Student t-distribution. /// - private double _location; + double _location; /// /// Keeps track of the degrees of freedom for the Student t-distribution. /// - private double _dof; + double _dof; /// /// Keeps track of the scale for the Student t-distribution. /// - private double _scale; + double _scale; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the StudentT class. This is a Student t-distribution with location 0.0 /// scale 1.0 and degrees of freedom 1. The distribution will /// be initialized with the default random number generator. /// - public StudentT() : this(0.0, 1.0, 1.0) + public StudentT() + : this(0.0, 1.0, 1.0) { } @@ -111,7 +112,7 @@ namespace MathNet.Numerics.Distributions /// The scale of the Student t-distribution. /// The degrees of freedom for the Student t-distribution. /// true when the parameters are valid, false otherwise. - private static bool IsValidParameterSet(double location, double scale, double dof) + static bool IsValidParameterSet(double location, double scale, double dof) { if (scale <= 0.0 || dof <= 0.0 || Double.IsNaN(scale) || Double.IsNaN(location) || Double.IsNaN(dof)) { @@ -128,7 +129,7 @@ namespace MathNet.Numerics.Distributions /// The scale of the Student t-distribution. /// The degrees of freedom for the Student t-distribution. /// When the parameters don't pass the function. - private void SetParameters(double location, double scale, double dof) + void SetParameters(double location, double scale, double dof) { if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, dof)) { @@ -145,15 +146,9 @@ namespace MathNet.Numerics.Distributions /// public double Location { - get - { - return _location; - } + get { return _location; } - set - { - SetParameters(value, _scale, _dof); - } + set { SetParameters(value, _scale, _dof); } } /// @@ -161,15 +156,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(_location, value, _dof); - } + set { SetParameters(_location, value, _dof); } } /// @@ -177,15 +166,9 @@ namespace MathNet.Numerics.Distributions /// public double DegreesOfFreedom { - get - { - return _dof; - } + get { return _dof; } - set - { - SetParameters(_location, _scale, value); - } + set { SetParameters(_location, _scale, value); } } #region IDistribution implementation @@ -195,10 +178,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -216,10 +196,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _dof > 1.0 ? _location : Double.NaN; - } + get { return _dof > 1.0 ? _location : Double.NaN; } } /// @@ -305,10 +282,7 @@ namespace MathNet.Numerics.Distributions /// public double Mode { - get - { - return _location; - } + get { return _location; } } /// @@ -316,10 +290,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return _location; - } + get { return _location; } } /// @@ -327,10 +298,7 @@ namespace MathNet.Numerics.Distributions /// public double Minimum { - get - { - return Double.NegativeInfinity; - } + get { return Double.NegativeInfinity; } } /// @@ -338,10 +306,7 @@ namespace MathNet.Numerics.Distributions /// public double Maximum { - get - { - return Double.PositiveInfinity; - } + get { return Double.PositiveInfinity; } } /// @@ -359,9 +324,9 @@ namespace MathNet.Numerics.Distributions var d = (x - _location) / _scale; return Math.Exp(SpecialFunctions.GammaLn((_dof + 1.0) / 2.0) - SpecialFunctions.GammaLn(_dof / 2.0)) - * Math.Pow(1.0 + (d * d / _dof), -0.5 * (_dof + 1.0)) - / Math.Sqrt(_dof * Math.PI) - / _scale; + * Math.Pow(1.0 + (d * d / _dof), -0.5 * (_dof + 1.0)) + / Math.Sqrt(_dof * Math.PI) + / _scale; } /// @@ -379,9 +344,9 @@ namespace MathNet.Numerics.Distributions var d = (x - _location) / _scale; return SpecialFunctions.GammaLn((_dof + 1.0) / 2.0) - - (0.5 * ((_dof + 1.0) * Math.Log(1.0 + (d * d / _dof)))) - - SpecialFunctions.GammaLn(_dof / 2.0) - - (0.5 * Math.Log(_dof * Math.PI)) - Math.Log(_scale); + - (0.5 * ((_dof + 1.0) * Math.Log(1.0 + (d * d / _dof)))) + - SpecialFunctions.GammaLn(_dof / 2.0) + - (0.5 * Math.Log(_dof * Math.PI)) - Math.Log(_scale); } /// @@ -403,13 +368,32 @@ namespace MathNet.Numerics.Distributions return x <= _location ? ib : 1.0 - ib; } + #endregion + + /// + /// Samples student-t distributed random variables. + /// + /// The algorithm is method 2 in section 5, chapter 9 + /// in L. Devroye's "Non-Uniform Random Variate Generation" + /// The random number generator to use. + /// The location of the Student t-distribution. + /// The scale of the Student t-distribution. + /// The degrees of freedom for the standard student-t distribution. + /// a random number from the standard student-t distribution. + internal static double SampleUnchecked(Random rnd, double location, double scale, double dof) + { + var n = Normal.SampleUncheckedBoxMuller(rnd).Item1; + var g = Gamma.SampleUnchecked(rnd, 0.5 * dof, 0.5); + return location + (scale * n * Math.Sqrt(dof / g)); + } + /// /// Generates a sample from the Student t-distribution. /// /// a sample from the distribution. public double Sample() { - return _location + (_scale * Sample(RandomSource, _dof)); + return SampleUnchecked(RandomSource, _location, _scale, _dof); } /// @@ -420,12 +404,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return _location + (_scale * Sample(RandomSource, _dof)); + yield return SampleUnchecked(RandomSource, _location, _scale, _dof); } } - #endregion - /// /// Generates a sample from the Student t-distribution. /// @@ -441,7 +423,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return location + (scale * Sample(rng, dof)); + return SampleUnchecked(rng, location, scale, dof); } /// @@ -461,23 +443,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return location + (scale * Sample(rng, dof)); + yield return SampleUnchecked(rng, location, scale, dof); } } - - /// - /// Samples standard student-t distributed random variables. - /// - /// The algorithm is method 2 in section 5, chapter 9 - /// in L. Devroye's "Non-Uniform Random Variate Generation" - /// The random number generator to use. - /// The degrees of freedom for the standard student-t distribution. - /// a random number from the standard student-t distribution. - internal static double Sample(Random rnd, double dof) - { - var n = Normal.SampleBoxMuller(rnd).Item1; - var g = Gamma.Sample(rnd, 0.5 * dof, 0.5); - return Math.Sqrt(dof / g) * n; - } } } diff --git a/src/Numerics/Distributions/Continuous/Weibull.cs b/src/Numerics/Distributions/Continuous/Weibull.cs index c970f98b..c3ccd589 100644 --- a/src/Numerics/Distributions/Continuous/Weibull.cs +++ b/src/Numerics/Distributions/Continuous/Weibull.cs @@ -50,12 +50,12 @@ namespace MathNet.Numerics.Distributions /// /// Weibull shape parameter. /// - private double _shape; + double _shape; /// /// Weibull inverse scale parameter. /// - private double _scale; + double _scale; /// /// Reusable intermediate result 1 / ( ^ ) @@ -64,12 +64,12 @@ namespace MathNet.Numerics.Distributions /// By caching this parameter we can get slightly better numerics precision /// in certain constellations without any additional computations. /// - private double _scalePowShapeInv; + double _scalePowShapeInv; /// /// The distribution's random number generator. /// - private Random _random; + Random _random; /// /// Initializes a new instance of the Weibull class. @@ -97,7 +97,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Weibull distribution. /// The scale of the Weibull distribution. /// true when the parameters positive valid floating point numbers, false otherwise. - private static bool IsValidParameterSet(double shape, double scale) + static bool IsValidParameterSet(double shape, double scale) { if (shape <= 0.0 || scale <= 0.0 || Double.IsNaN(shape) || Double.IsNaN(scale)) { @@ -113,7 +113,7 @@ namespace MathNet.Numerics.Distributions /// The shape of the Weibull distribution. /// The inverse scale of the Weibull distribution. /// When the parameters don't pass the function. - private void SetParameters(double shape, double scale) + void SetParameters(double shape, double scale) { if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) { @@ -130,15 +130,9 @@ namespace MathNet.Numerics.Distributions /// public double Shape { - get - { - return _shape; - } + get { return _shape; } - set - { - SetParameters(value, _scale); - } + set { SetParameters(value, _scale); } } /// @@ -146,15 +140,9 @@ namespace MathNet.Numerics.Distributions /// public double Scale { - get - { - return _scale; - } + get { return _scale; } - set - { - SetParameters(_shape, value); - } + set { SetParameters(_shape, value); } } #region IDistribution implementation @@ -164,10 +152,7 @@ namespace MathNet.Numerics.Distributions /// public Random RandomSource { - get - { - return _random; - } + get { return _random; } set { @@ -185,10 +170,7 @@ namespace MathNet.Numerics.Distributions /// public double Mean { - get - { - return _scale * SpecialFunctions.Gamma(1.0 + (1.0 / _shape)); - } + get { return _scale * SpecialFunctions.Gamma(1.0 + (1.0 / _shape)); } } /// @@ -196,10 +178,7 @@ namespace MathNet.Numerics.Distributions /// public double Variance { - get - { - return (_scale * _scale * SpecialFunctions.Gamma(1.0 + (2.0 / _shape))) - (Mean * Mean); - } + get { return (_scale * _scale * SpecialFunctions.Gamma(1.0 + (2.0 / _shape))) - (Mean * Mean); } } /// @@ -207,10 +186,7 @@ namespace MathNet.Numerics.Distributions /// public double StdDev { - get - { - return Math.Sqrt(Variance); - } + get { return Math.Sqrt(Variance); } } /// @@ -218,10 +194,7 @@ namespace MathNet.Numerics.Distributions /// public double Entropy { - get - { - return (Constants.EulerMascheroni * (1.0 - (1.0 / _shape))) + Math.Log(_scale / _shape) + 1.0; - } + get { return (Constants.EulerMascheroni * (1.0 - (1.0 / _shape))) + Math.Log(_scale / _shape) + 1.0; } } /// @@ -238,6 +211,7 @@ namespace MathNet.Numerics.Distributions return ((_scale * _scale * _scale * SpecialFunctions.Gamma(1.0 + (3.0 / _shape))) - (3.0 * sigma2 * mu) - (mu * mu * mu)) / sigma3; } } + #endregion #region IContinuousDistribution implementation @@ -263,10 +237,7 @@ namespace MathNet.Numerics.Distributions /// public double Median { - get - { - return _scale * Math.Pow(Constants.Ln2, 1.0 / _shape); - } + get { return _scale * Math.Pow(Constants.Ln2, 1.0 / _shape); } } /// @@ -340,13 +311,29 @@ namespace MathNet.Numerics.Distributions return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, _shape) * _scalePowShapeInv); } + #endregion + + /// + /// Generates one sample from the Weibull distribution. This method doesn't perform + /// any parameter checks. + /// + /// The random number generator to use. + /// The shape of the Weibull distribution. + /// The scale of the Weibull distribution. + /// A sample from a Weibull distributed random variable. + internal static double SampleUnchecked(Random rnd, double shape, double scale) + { + var x = rnd.NextDouble(); + return scale * Math.Pow(-Math.Log(x), 1.0 / shape); + } + /// /// Generates a sample from the Weibull distribution. /// /// a sample from the distribution. public double Sample() { - return SampleWeibull(RandomSource, _shape, _scale); + return SampleUnchecked(RandomSource, _shape, _scale); } /// @@ -357,12 +344,10 @@ namespace MathNet.Numerics.Distributions { while (true) { - yield return SampleWeibull(RandomSource, _shape, _scale); + yield return SampleUnchecked(RandomSource, _shape, _scale); } } - #endregion - /// /// Generates a sample from the Weibull distribution. /// @@ -377,7 +362,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - return SampleWeibull(rng, shape, scale); + return SampleUnchecked(rng, shape, scale); } /// @@ -396,22 +381,8 @@ namespace MathNet.Numerics.Distributions while (true) { - yield return SampleWeibull(rng, shape, scale); + yield return SampleUnchecked(rng, shape, scale); } } - - /// - /// Generates one sample from the Weibull distribution. This method doesn't perform - /// any parameter checks. - /// - /// The random number generator to use. - /// The shape of the Weibull distribution. - /// The scale of the Weibull distribution. - /// A sample from a Weibull distributed random variable. - internal static double SampleWeibull(Random rnd, double shape, double scale) - { - var x = rnd.NextDouble(); - return scale * Math.Pow(-Math.Log(x), 1.0 / shape); - } } } diff --git a/src/Numerics/Distributions/Multivariate/MatrixNormal.cs b/src/Numerics/Distributions/Multivariate/MatrixNormal.cs index 0070d456..7957ce8d 100644 --- a/src/Numerics/Distributions/Multivariate/MatrixNormal.cs +++ b/src/Numerics/Distributions/Multivariate/MatrixNormal.cs @@ -326,7 +326,7 @@ namespace MathNet.Numerics.Distributions var v = new DenseVector(count, 0.0); for (var d = 0; d < count; d += 2) { - var sample = Normal.SampleBoxMuller(rnd); + var sample = Normal.SampleUncheckedBoxMuller(rnd); v[d] = sample.Item1; if (d + 1 < count) {