From 453023c35beb509ee46ad256de43865c9a76b7d6 Mon Sep 17 00:00:00 2001 From: jvangael Date: Thu, 30 Jul 2009 19:16:28 +0800 Subject: [PATCH] Added Gamma & Beta distribution unit tests. Signed-off-by: Christoph Ruegg --- .../CommonDistributionTests.cs | 83 +++ .../DistributionTests/Continuous/BetaTests.cs | 378 +++++++++++++ .../Continuous/ContinuousUniformTests.cs | 21 +- .../Continuous/GammaTests.cs | 393 ++++++++++++++ .../Continuous/NormalTests.cs | 21 +- .../Managed.UnitTests.csproj | 3 + src/Managed/Distributions/Continuous/Beta.cs | 324 +++++++++++ src/Managed/Distributions/Continuous/Gamma.cs | 504 ++++++++++++++++++ src/Managed/Managed.csproj | 2 + src/Managed/SpecialFunctions.cs | 32 ++ src/Native.UnitTests/Native.UnitTests.csproj | 6 + src/Native/Native.csproj | 6 + 12 files changed, 1733 insertions(+), 40 deletions(-) create mode 100644 src/Managed.UnitTests/DistributionTests/CommonDistributionTests.cs create mode 100644 src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs create mode 100644 src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs create mode 100644 src/Managed/Distributions/Continuous/Beta.cs create mode 100644 src/Managed/Distributions/Continuous/Gamma.cs diff --git a/src/Managed.UnitTests/DistributionTests/CommonDistributionTests.cs b/src/Managed.UnitTests/DistributionTests/CommonDistributionTests.cs new file mode 100644 index 00000000..8adf6d33 --- /dev/null +++ b/src/Managed.UnitTests/DistributionTests/CommonDistributionTests.cs @@ -0,0 +1,83 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.DistributionTests +{ + using System; + using System.Linq; + using MbUnit.Framework; + using MathNet.Numerics.Distributions; + + [TestFixture] + public class CommonDistributionTests + { + private IDistribution[] dists; + + [SetUp] + public void SetupDistributions() + { + dists = new IDistribution[4]; + + dists[0] = new Beta(1.0, 1.0); + dists[1] = new ContinuousUniform(0.0, 1.0); + dists[2] = new Gamma(1.0, 1.0); + dists[3] = new Normal(0.0, 1.0); + } + + [Test] + [Row(0)] + [Row(1)] + [Row(2)] + [Row(3)] + public void CanCreateNormal(int i) + { + Assert.IsNotNull(dists[i].RandomSource); + } + + [Test] + [Row(0)] + [Row(1)] + [Row(2)] + [Row(3)] + public void CanSetRandomSource(int i) + { + dists[i].RandomSource = new Random(); + } + + [Test] + [Row(0)] + [Row(1)] + [Row(2)] + [Row(3)] + [ExpectedException(typeof(ArgumentNullException))] + public void FailSetRandomSourceWithNullReference(int i) + { + dists[i].RandomSource = null; + } + } +} diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs new file mode 100644 index 00000000..ebc129eb --- /dev/null +++ b/src/Managed.UnitTests/DistributionTests/Continuous/BetaTests.cs @@ -0,0 +1,378 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.DistributionTests +{ + using System; + using System.Linq; + using MbUnit.Framework; + using MathNet.Numerics.Distributions; + + [TestFixture] + public class BetaTests + { + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + } + + [Test, MultipleAsserts] + [Row(0.0, 0.0)] + [Row(0.0, 0.1)] + [Row(1.0, 0.0)] + [Row(1.0, 1.0)] + [Row(9.0, 1.0)] + [Row(5.0, 100.0)] + [Row(1.0, Double.PositiveInfinity)] + [Row(Double.PositiveInfinity, 1.0)] + [Row(0.0, Double.PositiveInfinity)] + [Row(Double.PositiveInfinity, 0.0)] + public void CanCreateBeta(double a, double b) + { + var n = new Beta(a, b); + AssertEx.AreEqual(a, n.A); + AssertEx.AreEqual(b, n.B); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(Double.NaN, 1.0)] + [Row(1.0, Double.NaN)] + [Row(Double.NaN, Double.NaN)] + [Row(1.0, -1.0)] + [Row(-1.0, 1.0)] + [Row(-1.0, -1.0)] + public void BetaCreateFailsWithBadParameters(double a, double b) + { + var n = new Beta(a, b); + } + + [Test] + public void ValidateToString() + { + var n = new Beta(1.0, 2.0); + AssertEx.AreEqual("Beta(A = 1, B = 2)", n.ToString()); + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetShapeA(double a) + { + var n = new Beta(1.0, 1.0); + n.A = a; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetShapeAFailsWithNegativeA() + { + var n = new Beta(1.0, 1.0); + n.A = -1.0; + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetShapeB(double b) + { + var n = new Beta(1.0, 1.0); + n.B = b; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetShapeBFailsWithNegativeB() + { + var n = new Beta(1.0, 1.0); + n.B = -1.0; + } + + [Test] + [Row(0.0, 0.0, 0.5)] + [Row(0.0, 0.1, 0.1)] + [Row(1.0, 0.0, 1.0)] + [Row(1.0, 1.0, 0.5)] + [Row(9.0, 1.0, 0.9)] + [Row(5.0, 100.0, 0.047619047619047619047616)] + [Row(1.0, Double.PositiveInfinity, 1.0)] + [Row(Double.PositiveInfinity, 1.0, 0.0)] + [Row(0.0, Double.PositiveInfinity, 1.0)] + [Row(Double.PositiveInfinity, 0.0, 0.0)] + public void ValidateMean(double a, double b, double mean) + { + var n = new Beta(a, b); + AssertEx.AreEqual(mean, n.Mean); + } + + [Test] + [Row(0.0, 0.0, 0.5)] + [Row(0.0, 0.1, 0.1)] + [Row(1.0, 0.0, 1.0)] + [Row(1.0, 1.0, 0.0)] + [Row(9.0, 1.0, -1.3083356884473304939016015849561625204060922267565917)] + [Row(5.0, 100.0, -2.5201623187602743679459255108827601222133603091753153)] + [Row(1.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 0.0)] + [Row(0.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 0.0)] + public void ValidateEntropy(double a, double b, double entropy) + { + var n = new Beta(a, b); + AssertEx.AreEqual(entropy, n.Entropy); + } + + [Test] + [Row(0.0, 0.0, 0.0)] + [Row(0.0, 0.1, 2.0)] + [Row(1.0, 0.0, -2.0)] + [Row(1.0, 1.0, 0.0)] + [Row(9.0, 1.0, -1.4740554623801777107177478829647496373009282424841579)] + [Row(5.0, 100.0, 0.81759410927553430354583159143895018978562196953345572)] + [Row(1.0, Double.PositiveInfinity, 2.0)] + [Row(Double.PositiveInfinity, 1.0, -2.0)] + [Row(0.0, Double.PositiveInfinity, 2.0)] + [Row(Double.PositiveInfinity, 0.0, -2.0)] + public void ValidateSkewness(double a, double b, double skewness) + { + var n = new Beta(a, b); + AssertEx.AreEqual(skewness, n.Skewness); + } + + [Test] + [Row(0.0, 0.0, 0.5)] + [Row(0.0, 0.1, 1.0)] + [Row(1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 0.5)] + [Row(9.0, 1.0, 1.0)] + [Row(5.0, 100.0, 0.038834951456310676243255386452801758423447608947753906)] + [Row(1.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 1.0)] + [Row(0.0, Double.PositiveInfinity, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 1.0)] + public void ValidateMode(double a, double b, double mode) + { + var n = new Beta(a, b); + AssertEx.AreEqual(mode, n.Mode); + } + + [Test] + [ExpectedException(typeof(NotSupportedException))] + [Row(0.0, 0.0)] + [Row(0.0, 0.1)] + [Row(1.0, 0.0)] + [Row(1.0, 1.0)] + [Row(9.0, 1.0)] + [Row(5.0, 100.0)] + [Row(1.0, Double.PositiveInfinity)] + [Row(Double.PositiveInfinity, 1.0)] + [Row(0.0, Double.PositiveInfinity)] + [Row(Double.PositiveInfinity, 0.0)] + public void ValidateMedian(double a, double b) + { + var n = new Beta(a, 1.0); + var m = n.Median; + } + + [Test] + public void ValidateMinimum() + { + var n = new Beta(1.0, 1.0); + AssertEx.AreEqual(0.0, n.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var n = new Beta(1.0, 1.0); + AssertEx.AreEqual(1.0, n.Maximum); + } + + [Test] + public void CanSampleStatic() + { + var d = Beta.Sample(new Random(), 2.0, 3.0); + } + + [Test] + public void CanSampleSequenceStatic() + { + var ied = Beta.Samples(new Random(), 2.0, 3.0); + var arr = ied.Take(5).ToArray(); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleStatic() + { + var d = Beta.Sample(new Random(), 1.0, -1.0); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleSequenceStatic() + { + var ied = Beta.Samples(new Random(), 1.0, -1.0).First(); + } + + [Test] + public void CanSample() + { + var n = new Beta(2.0, 3.0); + var d = n.Sample(); + } + + [Test] + public void CanSampleSequence() + { + var n = new Beta(2.0, 3.0); + var ied = n.Samples(); + var e = ied.Take(5).ToArray(); + } + + [Test] + [Row(0.0, 0.0, 0.0, Double.PositiveInfinity)] + [Row(0.0, 0.0, 0.5, 0.0)] + [Row(0.0, 0.0, 1.0, Double.PositiveInfinity)] + [Row(0.0, 0.1, 0.0, Double.PositiveInfinity)] + [Row(0.0, 0.1, 0.5, 0.0)] + [Row(0.0, 0.1, 1.0, 0.0)] + [Row(1.0, 0.0, 0.0, 0.0)] + [Row(1.0, 0.0, 0.5, 0.0)] + [Row(1.0, 0.0, 1.0, Double.PositiveInfinity)] + [Row(1.0, 1.0, 0.0, 1.0)] + [Row(1.0, 1.0, 0.5, 1.0)] + [Row(1.0, 1.0, 1.0, 1.0)] + [Row(9.0, 1.0, 0.0, 0.0)] + [Row(9.0, 1.0, 0.5, 0.035155378090821160189479427593561667617370600927556366)] + [Row(9.0, 1.0, 1.0, 8.9997767912502170085067334639517869100468738374544298)] + [Row(5.0, 100.0, 0.0, 0.0)] + [Row(5.0, 100.0, 0.5, 1.0881845516040810386311829462908430145307026037926335e-21)] + [Row(5.0, 100.0, 1.0, 0.0)] + [Row(1.0, Double.PositiveInfinity, 0.0, Double.PositiveInfinity)] + [Row(1.0, Double.PositiveInfinity, 0.5, 0.0)] + [Row(1.0, Double.PositiveInfinity, 1.0, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 0.0, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 0.5, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 1.0, Double.PositiveInfinity)] + [Row(0.0, Double.PositiveInfinity, 0.0, Double.PositiveInfinity)] + [Row(0.0, Double.PositiveInfinity, 0.5, 0.0)] + [Row(0.0, Double.PositiveInfinity, 1.0, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 0.0, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 0.5, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 1.0, Double.PositiveInfinity)] + public void ValidateDensity(double a, double b, double x, double pdf) + { + var n = new Beta(a, b); + AssertHelpers.AlmostEqual(pdf, n.Density(x), 15); + } + + [Test] + [Row(0.0, 0.0, 0.0, Double.PositiveInfinity)] + [Row(0.0, 0.0, 0.5, Double.NegativeInfinity)] + [Row(0.0, 0.0, 1.0, Double.PositiveInfinity)] + [Row(0.0, 0.1, 0.0, Double.PositiveInfinity)] + [Row(0.0, 0.1, 0.5, Double.NegativeInfinity)] + [Row(0.0, 0.1, 1.0, Double.NegativeInfinity)] + [Row(1.0, 0.0, 0.0, Double.NegativeInfinity)] + [Row(1.0, 0.0, 0.5, Double.NegativeInfinity)] + [Row(1.0, 0.0, 1.0, Double.PositiveInfinity)] + [Row(1.0, 1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 0.5, 0.0)] + [Row(1.0, 1.0, 1.0, 0.0)] + [Row(9.0, 1.0, 0.0, Double.NegativeInfinity)] + [Row(9.0, 1.0, 0.5, -3.3479528671433430925473664978203611353090199592365404)] + [Row(9.0, 1.0, 1.0, 2.197199776056471990627201569939718235188380979206923)] + [Row(5.0, 100.0, 0.0, Double.NegativeInfinity)] + [Row(5.0, 100.0, 0.5, -51.447830024537682154565870837960406410586196074573801)] + [Row(5.0, 100.0, 1.0, Double.NegativeInfinity)] + [Row(1.0, Double.PositiveInfinity, 0.0, Double.PositiveInfinity)] + [Row(1.0, Double.PositiveInfinity, 0.5, Double.NegativeInfinity)] + [Row(1.0, Double.PositiveInfinity, 1.0, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 1.0, 0.0, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 1.0, 0.5, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 1.0, 1.0, Double.PositiveInfinity)] + [Row(0.0, Double.PositiveInfinity, 0.0, Double.PositiveInfinity)] + [Row(0.0, Double.PositiveInfinity, 0.5, Double.NegativeInfinity)] + [Row(0.0, Double.PositiveInfinity, 1.0, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 0.0, 0.0, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 0.0, 0.5, Double.NegativeInfinity)] + [Row(Double.PositiveInfinity, 0.0, 1.0, Double.PositiveInfinity)] + public void ValidateDensityLn(double a, double b, double x, double pdfln) + { + var n = new Beta(a, b); + AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 15); + } + + [Test] + [Row(0.0, 0.0, 0.0, 0.5)] + [Row(0.0, 0.0, 0.5, 0.5)] + [Row(0.0, 0.0, 1.0, 1.0)] + [Row(0.0, 0.1, 0.0, 1.0)] + [Row(0.0, 0.1, 0.5, 1.0)] + [Row(0.0, 0.1, 1.0, 1.0)] + [Row(1.0, 0.0, 0.0, 0.0)] + [Row(1.0, 0.0, 0.5, 0.0)] + [Row(1.0, 0.0, 1.0, 1.0)] + [Row(1.0, 1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 0.5, 0.5)] + [Row(1.0, 1.0, 1.0, 1.0)] + [Row(9.0, 1.0, 0.0, 0.0)] + [Row(9.0, 1.0, 0.5, 0.00195313)] + [Row(9.0, 1.0, 1.0, 1.0)] + [Row(5.0, 100.0, 0.0, 0.0)] + [Row(5.0, 100.0, 0.5, 1.0)] + [Row(5.0, 100.0, 1.0, 1.0)] + [Row(1.0, Double.PositiveInfinity, 0.0, 1.0)] + [Row(1.0, Double.PositiveInfinity, 0.5, 1.0)] + [Row(1.0, Double.PositiveInfinity, 1.0, 1.0)] + [Row(Double.PositiveInfinity, 1.0, 0.0, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 0.5, 0.0)] + [Row(Double.PositiveInfinity, 1.0, 1.0, 1.0)] + [Row(0.0, Double.PositiveInfinity, 0.0, 1.0)] + [Row(0.0, Double.PositiveInfinity, 0.5, 1.0)] + [Row(0.0, Double.PositiveInfinity, 1.0, 1.0)] + [Row(Double.PositiveInfinity, 0.0, 0.0, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 0.5, 0.0)] + [Row(Double.PositiveInfinity, 0.0, 1.0, 1.0)] + public void ValidateCumulativeDistribution(double a, double b, double x, double cdf) + { + var n = new Beta(a, b); + AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 15); + } + } +} diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/ContinuousUniformTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/ContinuousUniformTests.cs index daf8fbe0..5d0be8b3 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/ContinuousUniformTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/ContinuousUniformTests.cs @@ -83,21 +83,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests AssertEx.AreEqual("ContinuousUniform(Lower = 1, Upper = 2)", n.ToString()); } - [Test] - public void CanGetRandomSource() - { - var n = new ContinuousUniform(); - var rs = n.RandomSource; - Assert.IsNotNull(rs); - } - - [Test] - public void CanSetRandomSource() - { - var n = new ContinuousUniform(); - n.RandomSource = new Random(); - } - [Test] [Row(-10.0)] [Row(-0.0)] @@ -272,11 +257,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests public void CanSampleSequenceStatic() { var ied = ContinuousUniform.Samples(new Random(), 0.0, 1.0); - var e = ied.GetEnumerator(); - e.MoveNext(); - var d = e.Current; - e.MoveNext(); - var g = e.Current; + var arr = ied.Take(5).ToArray(); } [Test] diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs new file mode 100644 index 00000000..ae7557d4 --- /dev/null +++ b/src/Managed.UnitTests/DistributionTests/Continuous/GammaTests.cs @@ -0,0 +1,393 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.DistributionTests +{ + using System; + using System.Linq; + using MbUnit.Framework; + using MathNet.Numerics.Distributions; + + [TestFixture] + public class GammaTests + { + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + } + + [Test, MultipleAsserts] + [Row(0.0, 0.0)] + [Row(1.0, 0.1)] + [Row(1.0, 1.0)] + [Row(10.0, 10.0)] + [Row(10.0, 1.0)] + [Row(10.0, Double.PositiveInfinity)] + public void CanCreateGamma(double shape, double invScale) + { + var n = new Gamma(shape, invScale); + AssertEx.AreEqual(shape, n.Shape); + AssertEx.AreEqual(invScale, n.InvScale); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(Double.NaN, 1.0)] + [Row(1.0, Double.NaN)] + [Row(Double.NaN, Double.NaN)] + [Row(1.0, -1.0)] + [Row(-1.0, 1.0)] + [Row(-1.0, -1.0)] + public void GammaCreateFailsWithBadParameters(double shape, double invScale) + { + var n = new Gamma(shape, invScale); + } + + [Test, MultipleAsserts] + [Row(0.0, 0.0)] + [Row(1.0, 0.1)] + [Row(1.0, 1.0)] + [Row(10.0, 10.0)] + [Row(10.0, 1.0)] + [Row(10.0, Double.PositiveInfinity)] + public void CanCreateGammaWithShapeInvScale(double shape, double invScale) + { + var n = Gamma.WithShapeInvScale(shape, invScale); + AssertEx.AreEqual(shape, n.Shape); + AssertEx.AreEqual(invScale, n.InvScale); + } + + [Test, MultipleAsserts] + [Row(0.0, 0.0)] + [Row(1.0, 0.1)] + [Row(1.0, 1.0)] + [Row(10.0, 10.0)] + [Row(10.0, 1.0)] + [Row(10.0, Double.PositiveInfinity)] + public void CanCreateGammaWithShapeScale(double shape, double scale) + { + var n = Gamma.WithShapeScale(shape, scale); + AssertEx.AreEqual(shape, n.Shape); + AssertEx.AreEqual(scale, n.Scale); + } + + [Test] + public void ValidateToString() + { + var n = new Gamma(1.0, 2.0); + AssertEx.AreEqual("Gamma(Shape = 1, Inverse Scale = 2)", n.ToString()); + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetShape(double shape) + { + var n = new Gamma(1.0, 1.0); + n.Shape = shape; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetShapeFailsWithNegativeShape() + { + var n = new Gamma(1.0, 1.0); + n.Shape = -1.0; + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetScale(double scale) + { + var n = new Gamma(1.0, 1.0); + n.Scale = scale; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetScaleFailsWithNegativeScale() + { + var n = new Gamma(1.0, 1.0); + n.Scale = -1.0; + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetInvScale(double invScale) + { + var n = new Gamma(1.0, 1.0); + n.InvScale = invScale; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetInvScaleFailsWithNegativeInvScale() + { + var n = new Gamma(1.0, 1.0); + n.InvScale = -1.0; + } + + [Test] + [Row(0.0, 0.0, 0.0)] + [Row(1.0, 0.1, 10.0)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 1.0)] + [Row(10.0, 1.0, 10.0)] + [Row(10.0, Double.PositiveInfinity, 0.0)] + public void CanGetMean(double shape, double invScale, double mean) + { + var n = new Gamma(shape, invScale); + AssertEx.AreEqual(mean, n.Mean); + } + + [Test] + [Row(0.0, 0.0, 0.0)] + [Row(1.0, 0.1, 100.0)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 0.1)] + [Row(10.0, 1.0, 10.0)] + [Row(10.0, Double.PositiveInfinity, 0.0)] + public void CanGetVariance(double shape, double invScale, double var) + { + var n = new Gamma(shape, invScale); + AssertEx.AreEqual(var, n.Variance); + } + + [Test] + [Row(0.0, 0.0, 0.0)] + [Row(1.0, 0.1, 10.0)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 0.31622776601683794197697302588502426416723164097476643)] + [Row(10.0, 1.0, 3.1622776601683793319988935444327185337195551393252168)] + [Row(10.0, Double.PositiveInfinity, 0.0)] + public void CanGetStdDev(double shape, double invScale, double sdev) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(sdev, n.StdDev, 15); + } + + [Test] + [Row(0.0, 0.0, Double.PositiveInfinity)] + [Row(1.0, 0.1, 3.3025850929940456285068402234265387271634735938763824)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 0.23346908548693395836262094490967812177376750477943892)] + [Row(10.0, 1.0, 2.5360541784809796423806123995940423293748689934081866)] + [Row(10.0, Double.PositiveInfinity, 0.0)] + public void ValidateEntropy(double shape, double invScale, double entropy) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(entropy, n.Entropy, 15); + } + + [Test] + [Row(0.0, 0.0, Double.PositiveInfinity)] + [Row(1.0, 0.1, 2.0)] + [Row(1.0, 1.0, 2.0)] + [Row(10.0, 10.0, 0.63245553203367586639977870888654370674391102786504337)] + [Row(10.0, 1.0, 0.63245553203367586639977870888654370674391102786504337)] + [Row(10.0, Double.PositiveInfinity, 0.0)] + public void ValidateSkewness(double shape, double invScale, double skewness) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(skewness, n.Skewness, 15); + } + + [Test] + [Row(0.0, 0.0, Double.PositiveInfinity)] + [Row(1.0, 0.1, 0.0)] + [Row(1.0, 1.0, 0.0)] + [Row(10.0, 10.0, 0.9)] + [Row(10.0, 1.0, 9.0)] + [Row(10.0, Double.PositiveInfinity, 10.0)] + public void ValidateMode(double shape, double invScale, double mode) + { + var n = new Gamma(shape, invScale); + AssertEx.AreEqual(mode, n.Mode); + } + + [Test] + [ExpectedException(typeof(NotSupportedException))] + [Row(0.0, 0.0)] + [Row(1.0, 0.1)] + [Row(1.0, 1.0)] + [Row(10.0, 10.0)] + [Row(10.0, 1.0)] + [Row(10.0, Double.PositiveInfinity)] + public void ValidateMedian(double shape, double invScale, double mode) + { + var n = new Gamma(shape, invScale); + var median = n.Median; + } + + [Test] + public void ValidateMinimum() + { + var n = new Gamma(1.0,1.0); + AssertEx.AreEqual(0.0, n.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var n = new Gamma(1.0, 1.0); + AssertEx.AreEqual(System.Double.PositiveInfinity, n.Maximum); + } + + [Test] + [Row(0.0, 0.0, 0.0, 0.0)] + [Row(0.0, 0.0, 1.0, 0.0)] + [Row(0.0, 0.0, 10.0, 0.0)] + [Row(1.0, 0.1, 0.0, 0.10000000000000000555111512312578270211815834045410156)] + [Row(1.0, 0.1, 1.0, 0.099004983374916810660915381324116279472953858297127391)] + [Row(1.0, 0.1, 10.0, 0.036787944117144234201693506390001264039984687455876246)] + [Row(1.0, 1.0, 0.0, 1.0)] + [Row(1.0, 1.0, 1.0, 0.36787944117144232159552377016146086744581113103176804)] + [Row(1.0, 1.0, 10.0, 0.000045399929762484851535591515560550610237918088866564953)] + [Row(10.0, 10.0, 0.0, 0.0)] + [Row(10.0, 10.0, 1.0, 1.2511003572113329898476497894772544708420990097708588)] + [Row(10.0, 10.0, 10.0, 1.0251532120868705806216092933926141802686541811003037e-30)] + [Row(10.0, 1.0, 0.0, 0.0)] + [Row(10.0, 1.0, 1.0, 0.0000010137771196302974029859010421116095333052555418644397)] + [Row(10.0, 1.0, 10.0, 0.12511003572113329898476497894772544708420990097708601)] + [Row(10.0, Double.PositiveInfinity, 0.0, 0.0)] + [Row(10.0, Double.PositiveInfinity, 1.0, 0.0)] + [Row(10.0, Double.PositiveInfinity, 10.0, Double.PositiveInfinity)] + public void ValidateDensity(double shape, double invScale, double x, double pdf) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(pdf, n.Density(x), 15); + } + + [Test] + [Row(0.0, 0.0, 0.0, Double.NegativeInfinity)] + [Row(0.0, 0.0, 1.0, Double.NegativeInfinity)] + [Row(0.0, 0.0, 10.0, Double.NegativeInfinity)] + [Row(1.0, 0.1, 0.0, -2.3025850929940456285068402234265387271634735938763824)] + [Row(1.0, 0.1, 1.0, -2.312585092994045630449730516520562672904829013035318)] + [Row(1.0, 0.1, 10.0, -3.3025850929940456285068402234265387271634735938763824)] + [Row(1.0, 1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 1.0, -1.0)] + [Row(1.0, 1.0, 10.0, -10.0)] + [Row(10.0, 10.0, 0.0, Double.NegativeInfinity)] + [Row(10.0, 10.0, 1.0, 0.22402344985898722897219667227693591172986563062456522)] + [Row(10.0, 10.0, 10.0, -69.052710713194601614865880235563786219860220971716511)] + [Row(10.0, 1.0, 0.0, Double.NegativeInfinity)] + [Row(10.0, 1.0, 1.0, -13.801827480081469611207717874566706164281149255663166)] + [Row(10.0, 1.0, 10.0, -2.0785616431350584550457947824074282958712358580042068)] + [Row(10.0, Double.PositiveInfinity, 0.0, Double.NegativeInfinity)] + [Row(10.0, Double.PositiveInfinity, 1.0, Double.NegativeInfinity)] + [Row(10.0, Double.PositiveInfinity, 10.0, Double.PositiveInfinity)] + public void ValidateDensityLn(double shape, double invScale, double x, double pdfln) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 15); + } + + [Test] + public void CanSampleStatic() + { + var d = Gamma.Sample(new Random(), 1.0, 1.0); + } + + [Test] + public void CanSampleSequenceStatic() + { + var ied = Gamma.Samples(new Random(), 1.0, 1.0); + var arr = ied.Take(5).ToArray(); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleStatic() + { + var d = Normal.Sample(new Random(), 1.0, -1.0); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleSequenceStatic() + { + var ied = Normal.Samples(new Random(), 1.0, -1.0).First(); + } + + [Test] + public void CanSample() + { + var n = new Normal(); + var d = n.Sample(); + } + + [Test] + public void CanSampleSequence() + { + var n = new Normal(); + var ied = n.Samples(); + var e = ied.Take(5).ToArray(); + } + + [Test] + [Row(0.0, 0.0, 0.0, 0.0)] + [Row(0.0, 0.0, 1.0, 0.0)] + [Row(0.0, 0.0, 10.0, 0.0)] + [Row(1.0, 0.1, 0.0, 0.0)] + [Row(1.0, 0.1, 1.0, 0.095162581964040431858607615783064404690935346242622848)] + [Row(1.0, 0.1, 10.0, 0.63212055882855767840447622983853913255418886896823196)] + [Row(1.0, 1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 1.0, 0.63212055882855767840447622983853913255418886896823196)] + [Row(1.0, 1.0, 10.0, 0.99995460007023751514846440848443944938976208191113396)] + [Row(10.0, 10.0, 0.0, 0.0)] + [Row(10.0, 10.0, 1.0, 0.54207028552814779168583514294066541824736464003242184)] + [Row(10.0, 10.0, 10.0, 0.99999999999999999999999999999988746526039157266114706)] + [Row(10.0, 1.0, 0.0, 0.0)] + [Row(10.0, 1.0, 1.0, 0.00000011142547833872067735305068724025236288094949815466035)] + [Row(10.0, 1.0, 10.0, 0.54207028552814779168583514294066541824736464003242184)] + [Row(10.0, Double.PositiveInfinity, 0.0, 0.0)] + [Row(10.0, Double.PositiveInfinity, 1.0, 0.0)] + [Row(10.0, Double.PositiveInfinity, 10.0, 1.0)] + public void ValidateCumulativeDistribution(double shape, double invScale, double x, double cdf) + { + var n = new Gamma(shape, invScale); + AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 15); + } + } +} diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs index 5ba4b03e..82426a7c 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs @@ -128,21 +128,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests AssertEx.AreEqual("Normal(Mean = 1, StdDev = 2)", n.ToString()); } - [Test] - public void CanGetRandomSource() - { - var n = new Normal(); - var rs = n.RandomSource; - Assert.IsNotNull(rs); - } - - [Test] - public void CanSetRandomSource() - { - var n = new Normal(); - n.RandomSource = new Random(); - } - [Test] [Row(-0.0)] [Row(0.0)] @@ -338,11 +323,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests public void CanSampleSequenceStatic() { var ied = Normal.Samples(new Random(), 0.0, 1.0); - var e = ied.GetEnumerator(); - e.MoveNext(); - var d = e.Current; - e.MoveNext(); - var g = e.Current; + var arr = ied.Take(5).ToArray(); } [Test] diff --git a/src/Managed.UnitTests/Managed.UnitTests.csproj b/src/Managed.UnitTests/Managed.UnitTests.csproj index d987d51d..f140563f 100644 --- a/src/Managed.UnitTests/Managed.UnitTests.csproj +++ b/src/Managed.UnitTests/Managed.UnitTests.csproj @@ -62,7 +62,10 @@ + + + diff --git a/src/Managed/Distributions/Continuous/Beta.cs b/src/Managed/Distributions/Continuous/Beta.cs new file mode 100644 index 00000000..c205bd8f --- /dev/null +++ b/src/Managed/Distributions/Continuous/Beta.cs @@ -0,0 +1,324 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + /// + /// Implements the Beta distribution. For details about this distribution, see + /// Wikipedia - Beta distribution. + /// + /// The distribution will use the by default. + /// Users can get/set the random number generator by using the property. + /// The statistics classes will check all the incoming parameters whether they are in the allowed + /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters + /// to false, all parameter checks can be turned off. + public class Beta : IContinuousDistribution + { + /// + /// Beta shape parameter a. + /// + private double _shapeA; + + /// + /// Beta shape parameter b. + /// + private double _shapeB; + + /// + /// Initializes a new instance of the Beta distribution. + /// + /// The a shape parameter of the Beta distribution. + /// The b shape parameter of the Beta distribution. + /// If any of the Beta parameters are negative. + public Beta(double a, double b) + { + SetParameters(a, b); + RandomSource = new Random(); + } + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Beta(A = " + _shapeA + ", B = " + _shapeB + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// 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) + { + if (a < 0.0 || b < 0.0) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// 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) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _shapeA = a; + _shapeB = b; + } + + /// + /// Gets or sets the A shape parameter of the Beta distribution. + /// + public double A + { + get { return _shapeA; } + set { SetParameters(value, _shapeB); } + } + + /// + /// Gets or sets the B shape parameter of the Beta distribution. + /// + public double B + { + get { return _shapeB; } + set { SetParameters(_shapeA, value); } + } + + #region IDistribution implementation + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource { get; set; } + + /// + /// Gets the mean of the Beta distribution. + /// + public double Mean + { + get { return _shapeA / (_shapeA + _shapeB); } + } + + /// + /// Gets the variance of the Beta distribution. + /// + public double Variance + { + get { return (_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0)); } + } + + /// + /// Gets the standard deviation of the Beta distribution. + /// + public double StdDev + { + get { return Math.Sqrt((_shapeA * _shapeB) / ((_shapeA + _shapeB) * (_shapeA + _shapeB) * (_shapeA + _shapeB + 1.0))); } + } + + /// + /// Gets the entropy of the Beta distribution. + /// + public double Entropy + { + get + { + 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); + } + } + + /// + /// Gets the skewness of the Beta distribution. + /// + public double Skewness + { + get + { + return 2.0 * (_shapeB - _shapeA) * Math.Sqrt(_shapeA + _shapeB + 1.0) + / ((_shapeA + _shapeB + 2.0) * Math.Sqrt(_shapeA * _shapeB)); + } + } + #endregion + + #region IContinuousDistribution implementation + + /// + /// Gets the mode of the Beta distribution. + /// + public double Mode + { + get { return (_shapeA - 1) / (_shapeA + _shapeB - 2); } + } + + /// + /// Gets the median of the Beta distribution. + /// + public double Median + { + get { throw new NotSupportedException(); } + } + + /// + /// Gets the minimum of the Beta distribution. + /// + public double Minimum + { + get { return 0.0; } + } + + /// + /// Gets the maximum of the Beta distribution. + /// + public double Maximum + { + get { return 1.0; } + } + + /// + /// Computes the density of the Beta distribution. + /// + /// The location at which to compute the density. + /// the density at . + public double Density(double x) + { + double b = SpecialFunctions.Gamma(_shapeA + _shapeB) / (SpecialFunctions.Gamma(_shapeA) * SpecialFunctions.Gamma(_shapeB)); + return b * Math.Pow(x, _shapeA - 1.0) * Math.Pow(1.0 - x, _shapeB - 1.0); + } + + /// + /// Computes the log density of the Beta distribution. + /// + /// The location at which to compute the log density. + /// the log density at . + public double DensityLn(double x) + { + double b = SpecialFunctions.GammaLn(_shapeA + _shapeB) - SpecialFunctions.GammaLn(_shapeA) - SpecialFunctions.GammaLn(_shapeB); + return b + (_shapeA - 1.0)*Math.Log(x) + (_shapeB - 1.0)*Math.Log(1.0 - x); + } + + /// + /// Computes the cumulative distribution function of the Beta distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + public double CumulativeDistribution(double x) + { + return SpecialFunctions.BetaRegularized(_shapeA, _shapeB, x); + } + + /// + /// Generates a sample from the Beta distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleBeta(RandomSource, _shapeA, _shapeB); + } + + /// + /// Generates a sequence of samples from the Beta distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return SampleBeta(RandomSource, _shapeA, _shapeB); + } + } + #endregion + + /// + /// Generates a sample from the normal distribution using the Box-Muller algorithm. + /// + /// 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) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleBeta(rng, a, b); + } + + /// + /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. + /// + /// 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) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(a, b)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleBeta(rng, 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) + { + double x = Gamma.SampleGamma(rnd, a, 1.0); + double y = Gamma.SampleGamma(rnd, b, 1.0); + return x / (x + y); + } + } +} diff --git a/src/Managed/Distributions/Continuous/Gamma.cs b/src/Managed/Distributions/Continuous/Gamma.cs new file mode 100644 index 00000000..379cb71b --- /dev/null +++ b/src/Managed/Distributions/Continuous/Gamma.cs @@ -0,0 +1,504 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + /// + /// Implements the univariate Gamma distribution. For details about this distribution, see + /// Wikipedia - Gamma distribution. + /// + /// + /// The Gamma distribution is parametrized by a shape and inverse scale parameter. When we want + /// to specify a Gamma distribution which is a point distribution we set the shape parameter to be the + /// location of the point distribution and the inverse scale as positive infinity. + /// Random number generation for the Gamma distribution is based on the algorithm in: + /// "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang + /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372. + /// The distribution will use the by default. + /// Users can get/set the random number generator by using the property. + /// The statistics classes will check all the incoming parameters whether they are in the allowed + /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters + /// to false, all parameter checks can be turned off. + public class Gamma : IContinuousDistribution + { + /// + /// Gamma shape parameter. + /// + private double _shape; + + /// + /// Gamma inverse scale parameter. + /// + private double _invScale; + + /// + /// Initializes a new instance of the Gamma distribution. + /// + /// The shape of the Gamma distribution. + /// The inverse scale of the Gamma distribution. + public Gamma(double shape, double invScale) + { + SetParameters(shape, invScale); + RandomSource = new Random(); + } + + /// + /// Constructs a Gamma distribution from a shape and scale parameter. The distribution will + /// be initialized with the default random number generator. + /// + /// The shape of the Gamma distribution. + /// The scale of the Gamma distribution. + /// a normal distribution. + public static Gamma WithShapeScale(double shape, double scale) + { + return new Gamma(shape, 1.0/scale); + } + + /// + /// Constructs a Gamma distribution from a shape and inverse scale parameter. The distribution will + /// be initialized with the default random number generator. + /// + /// The shape of the Gamma distribution. + /// The inverse scale of the Gamma distribution. + /// a normal distribution. + public static Gamma WithShapeInvScale(double shape, double invScale) + { + return new Gamma(shape, invScale); + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "Gamma(Shape = " + _shape + ", Inverse Scale = " + _invScale + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// 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) + { + if (shape < 0.0 || invScale < 0.0 || Double.IsNaN(shape) || Double.IsNaN(invScale)) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// 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) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, invScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _shape = shape; + _invScale = invScale; + } + + /// + /// Gets or sets the shape of the Gamma distribution. + /// + public double Shape + { + get + { + return _shape; + } + + set + { + SetParameters(value, _invScale); + } + } + + /// + /// Gets or sets the scale of the Gamma distribution. + /// + public double Scale + { + get + { + return 1.0 / _invScale; + } + + set + { + SetParameters(_shape, 1.0/value); + } + } + + /// + /// Gets or sets the inverse scale of the Gamma distribution. + /// + public double InvScale + { + get + { + return _invScale; + } + + set + { + SetParameters(_shape, value); + } + } + + #region IDistribution implementation + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource { get; set; } + + /// + /// Gets the mean of the Gamma distribution. + /// + public double Mean + { + get + { + if (Double.IsPositiveInfinity(_invScale)) + { + return _shape; + } + else + { + return _shape / _invScale; + } + } + } + + /// + /// Gets the variance of the Gamma distribution. + /// + public double Variance + { + get + { + if (Double.IsPositiveInfinity(_invScale)) + { + return 0.0; + } + else + { + return _shape / (_invScale * _invScale); + } + } + } + + /// + /// Gets the standard deviation of the Gamma distribution. + /// + public double StdDev + { + get + { + if (Double.IsPositiveInfinity(_invScale)) + { + return 0.0; + } + else + { + return Math.Sqrt(_shape / (_invScale * _invScale)); + } + } + } + + /// + /// Gets the entropy of the Gamma distribution. + /// + public double Entropy + { + get + { + if (Double.IsPositiveInfinity(_invScale)) + { + return 0.0; + } + else + { + return _shape - Math.Log(_invScale) + SpecialFunctions.GammaLn(_shape) + (1.0 - _shape) * SpecialFunctions.DiGamma(_shape); + } + } + } + + /// + /// Gets the skewness of the Gamma distribution. + /// + public double Skewness + { + get { return 2.0 / Math.Sqrt(_shape); } + } + #endregion + + #region IContinuousDistribution implementation + + /// + /// Gets the mode of the Gamma distribution. + /// + public double Mode + { + get + { + if (Double.IsPositiveInfinity(_invScale)) + { + return _shape; + } + else + { + return (_shape - 1.0) / _invScale; + } + } + } + + /// + /// Gets the median of the Gamma distribution. + /// + public double Median + { + get { throw new NotSupportedException(); } + } + + /// + /// Gets the minimum of the Gamma distribution. + /// + public double Minimum + { + get { return 0.0; } + } + + /// + /// Gets the maximum of the Gamma distribution. + /// + public double Maximum + { + get { return Double.PositiveInfinity; } + } + + /// + /// Computes the density of the Gamma distribution. + /// + /// The location at which to compute the density. + /// the density at . + public double Density(double x) + { + if (Double.IsPositiveInfinity(_invScale)) + { + if (x == _shape) + { + return Double.PositiveInfinity; + } + else + { + return 0.0; + } + } + else + { + return Math.Pow(_invScale, _shape) * Math.Pow(x, _shape - 1.0) * Math.Exp(-_invScale * x) / SpecialFunctions.Gamma(_shape); + } + } + + /// + /// Computes the log density of the Gamma distribution. + /// + /// The location at which to compute the log density. + /// the log density at . + public double DensityLn(double x) + { + if (Double.IsPositiveInfinity(_invScale)) + { + if (x == _shape) + { + return Double.PositiveInfinity; + } + else + { + return Double.NegativeInfinity; + } + } + else + { + return _shape * Math.Log(_invScale) + (_shape - 1.0) * Math.Log(x) - _invScale * x - SpecialFunctions.GammaLn(_shape); + } + } + + /// + /// Computes the cumulative distribution function of the Gamma distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + public double CumulativeDistribution(double x) + { + if (Double.IsPositiveInfinity(_invScale)) + { + if (x >= _shape) + { + return 1.0; + } + else + { + return 0.0; + } + } + else + { + return SpecialFunctions.IncompleteGamma(_shape, x * _invScale, true); + } + } + + /// + /// 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); + } + } + #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 + /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372. + /// + /// The random number generator to use. + /// 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(System.Random rnd, double shape, double invScale) + { + if (Double.IsPositiveInfinity(invScale)) + { + return shape; + } + else + { + double a = shape; + double alphafix = 1.0; + + // Fix when alpha is less than one. + if (shape < 1.0) + { + a = shape + 1.0; + alphafix = System.Math.Pow(rnd.NextDouble(), 1.0 / shape); + } + + double d = a - 1.0 / 3.0; + double c = 1.0 / System.Math.Sqrt(9.0 * d); + while (true) + { + double x = Normal.Sample(rnd, 0.0, 1.0); + double v = 1.0 + c * x; + while (v <= 0.0) + { + x = Normal.Sample(rnd, 0.0, 1.0); + v = 1.0 + c * x; + } + v = v * v * v; + double u = rnd.NextDouble(); + x = x * x; + if (u < 1.0 - 0.0331 * x * x) + { + return alphafix * d * v / invScale; + } + if (System.Math.Log(u) < 0.5 * x + d * (1.0 - v + System.Math.Log(v))) + { + return alphafix * d * v / invScale; + } + } + } + } + } +} diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index e7f69772..069bd8f8 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -48,7 +48,9 @@ + + diff --git a/src/Managed/SpecialFunctions.cs b/src/Managed/SpecialFunctions.cs index d8bcc487..76e1e74e 100644 --- a/src/Managed/SpecialFunctions.cs +++ b/src/Managed/SpecialFunctions.cs @@ -58,5 +58,37 @@ namespace MathNet.Numerics return 0d; } + + + public static double BetaLn(double a, double b) + { + return Double.NaN; + } + + + public static double BetaRegularized(double a, double b, double x) + { + return Double.NaN; + } + + public static double DiGamma(double x) + { + return Double.NaN; + } + + public static double Gamma(double x) + { + return Double.NaN; + } + + public static double GammaLn(double x) + { + return Double.NaN; + } + + public static double IncompleteGamma(double x, double z, bool reg) + { + return Double.NaN; + } } } diff --git a/src/Native.UnitTests/Native.UnitTests.csproj b/src/Native.UnitTests/Native.UnitTests.csproj index ef565cd1..214689f2 100644 --- a/src/Native.UnitTests/Native.UnitTests.csproj +++ b/src/Native.UnitTests/Native.UnitTests.csproj @@ -68,9 +68,15 @@ ComplexTests\ComplexTest.cs + + DistributionTests\Continuous\BetaTests.cs + DistributionTests\Continuous\ContinuousUniformTests.cs + + DistributionTests\Continuous\GammaTests.cs + DistributionTests\Continuous\NormalTests.cs diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index 941b2ce6..d8d5ab04 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -53,9 +53,15 @@ Control.cs + + Distributions\Continuous\Beta.cs + Distributions\Continuous\ContinuousUniform.cs + + Distributions\Continuous\Gamma.cs + Distributions\Continuous\Normal.cs