diff --git a/src/Numerics/Distributions/InverseGamma.cs b/src/Numerics/Distributions/InverseGamma.cs index 35bbc784..de538751 100644 --- a/src/Numerics/Distributions/InverseGamma.cs +++ b/src/Numerics/Distributions/InverseGamma.cs @@ -30,6 +30,7 @@ using System; using System.Collections.Generic; +using System.Linq; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Distributions @@ -84,17 +85,6 @@ namespace MathNet.Numerics.Distributions return "InverseGamma(α = " + _shape + ", β = " + _scale + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The shape (α) of the distribution. Range: α > 0. - /// The scale (β) of the distribution. Range: β > 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double shape, double scale) - { - return shape > 0.0 && scale > 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -103,7 +93,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + if (shape <= 0.0 || scale <= 0.0 || Double.IsNaN(shape) || Double.IsNaN(scale)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -241,14 +231,10 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the density. /// the density at . + /// public double Density(double x) { - if (x >= 0.0) - { - return Math.Pow(_scale, _shape)*Math.Pow(x, -_shape - 1.0)*Math.Exp(-_scale/x)/SpecialFunctions.Gamma(_shape); - } - - return 0.0; + return x < 0.0 ? 0.0 : Math.Pow(_scale, _shape)*Math.Pow(x, -_shape - 1.0)*Math.Exp(-_scale/x)/SpecialFunctions.Gamma(_shape); } /// @@ -256,6 +242,7 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the log density. /// the log density at . + /// public double DensityLn(double x) { return Math.Log(Density(x)); @@ -266,42 +253,71 @@ namespace MathNet.Numerics.Distributions /// /// The location at which to compute the cumulative distribution function. /// the cumulative distribution at location . + /// public double CumulativeDistribution(double x) { return SpecialFunctions.GammaUpperRegularized(_shape, _scale/x); } /// - /// Samples the distribution. + /// Draws a random sample from the distribution. + /// + /// A random number from this distribution. + public double Sample() + { + return 1.0/Gamma.Sample(_random, _shape, _scale); + } + + /// + /// Generates a sequence of samples from the Cauchy distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return Gamma.Samples(_random, _shape, _scale).Select(z => 1.0/z); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. /// - /// The random number generator to use. /// The shape (α) of the distribution. Range: α > 0. /// The scale (β) of the distribution. Range: β > 0. - /// a random number from the distribution. - static double SampleUnchecked(System.Random rnd, double shape, double scale) + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double shape, double scale, double x) { - return 1.0/Gamma.Sample(rnd, shape, scale); + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return x < 0.0 ? 0.0 : Math.Pow(scale, shape)*Math.Pow(x, -shape - 1.0)*Math.Exp(-scale/x)/SpecialFunctions.Gamma(shape); } /// - /// Draws a random sample from the distribution. + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). /// - /// A random number from this distribution. - public double Sample() + /// The shape (α) of the distribution. Range: α > 0. + /// The scale (β) of the distribution. Range: β > 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double shape, double scale, double x) { - return SampleUnchecked(_random, _shape, _scale); + return Math.Log(PDF(shape, scale, x)); } /// - /// Generates a sequence of samples from the Cauchy distribution. + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). /// - /// a sequence of samples from the distribution. - public IEnumerable Samples() + /// The location at which to compute the cumulative distribution function. + /// The shape (α) of the distribution. Range: α > 0. + /// The scale (β) of the distribution. Range: β > 0. + /// the cumulative distribution at location . + /// + public static double CDF(double shape, double scale, double x) { - while (true) - { - yield return SampleUnchecked(_random, _shape, _scale); - } + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + return SpecialFunctions.GammaUpperRegularized(shape, scale/x); } /// @@ -313,12 +329,9 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - return SampleUnchecked(rnd, shape, scale); + return 1.0/Gamma.Sample(rnd, shape, scale); } /// @@ -330,15 +343,9 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - while (true) - { - yield return SampleUnchecked(rnd, shape, scale); - } + return Gamma.Samples(rnd, shape, scale).Select(z => 1.0/z); } } } diff --git a/src/UnitTests/DistributionTests/Continuous/InverseGammaTests.cs b/src/UnitTests/DistributionTests/Continuous/InverseGammaTests.cs index 3de4c931..41170c5c 100644 --- a/src/UnitTests/DistributionTests/Continuous/InverseGammaTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/InverseGammaTests.cs @@ -262,14 +262,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous public void ValidateDensity(double a, double b, double x) { var n = new InverseGamma(a, b); - if (x >= 0) - { - Assert.AreEqual(Math.Pow(b, a) * Math.Pow(x, -a - 1.0) * Math.Exp(-b / x) / SpecialFunctions.Gamma(a), n.Density(x)); - } - else - { - Assert.AreEqual(0.0, n.Density(x)); - } + double expected = Math.Pow(b, a)*Math.Pow(x, -a - 1.0)*Math.Exp(-b/x)/SpecialFunctions.Gamma(a); + Assert.AreEqual(expected, n.Density(x)); + Assert.AreEqual(expected, InverseGamma.PDF(a, b, x)); } /// @@ -290,7 +285,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous public void ValidateDensityLn(double a, double b, double x) { var n = new InverseGamma(a, b); - Assert.AreEqual(Math.Log(n.Density(x)), n.DensityLn(x)); + double expected = Math.Log(Math.Pow(b, a)*Math.Pow(x, -a - 1.0)*Math.Exp(-b/x)/SpecialFunctions.Gamma(a)); + Assert.AreEqual(expected, n.DensityLn(x)); + Assert.AreEqual(expected, InverseGamma.PDFLn(a, b, x)); } /// @@ -332,7 +329,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous public void ValidateCumulativeDistribution(double a, double b, double x) { var n = new InverseGamma(a, b); - Assert.AreEqual(SpecialFunctions.GammaUpperRegularized(a, b / x), n.CumulativeDistribution(x)); + double expected = SpecialFunctions.GammaUpperRegularized(a, b/x); + Assert.AreEqual(expected, n.CumulativeDistribution(x)); + Assert.AreEqual(expected, InverseGamma.CDF(a, b, x)); } } }