diff --git a/src/Numerics/Distributions/Continuous/Gamma.cs b/src/Numerics/Distributions/Continuous/Gamma.cs index d5489377..fa1f9171 100644 --- a/src/Numerics/Distributions/Continuous/Gamma.cs +++ b/src/Numerics/Distributions/Continuous/Gamma.cs @@ -528,9 +528,10 @@ namespace MathNet.Numerics.Distributions } /// - /// Sampling implementation based on: + /// 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. + /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372. + /// This method performs no parameter checks. /// /// The random number generator to use. /// The shape of the Gamma distribution. diff --git a/src/Numerics/Distributions/Continuous/Weibull.cs b/src/Numerics/Distributions/Continuous/Weibull.cs new file mode 100644 index 00000000..50d8eb67 --- /dev/null +++ b/src/Numerics/Distributions/Continuous/Weibull.cs @@ -0,0 +1,413 @@ +// +// 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 Weibull distribution. For details about this distribution, see + /// Wikipedia - Weibull distribution. + /// + /// + /// The Weibull distribution is parametrized by a shape and scale parameter. + /// 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 Weibull : IContinuousDistribution + { + /// + /// Weibull shape parameter. + /// + private double _shape; + + /// + /// Weibull inverse scale parameter. + /// + private double _scale; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Initializes a new instance of the Weibull class. + /// + /// The shape of the Weibull distribution. + /// The inverse scale of the Weibull distribution. + public Weibull(double shape, double scale) + { + SetParameters(shape, scale); + RandomSource = new Random(); + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "Weibull(Shape = " + _shape + ", Scale = " + _scale + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// 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) + { + if (shape <= 0.0 || scale <= 0.0 || Double.IsNaN(shape) || Double.IsNaN(scale)) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// 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) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _shape = shape; + _scale = scale; + } + + /// + /// Gets or sets the shape of the Weibull distribution. + /// + public double Shape + { + get + { + return _shape; + } + + set + { + SetParameters(value, _scale); + } + } + + /// + /// Gets or sets the scale of the Weibull distribution. + /// + public double Scale + { + get + { + return _scale; + } + + 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 + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } + + /// + /// Gets the mean of the Weibull distribution. + /// + public double Mean + { + get + { + return _scale * SpecialFunctions.Gamma(1.0 + 1.0 / _shape); + } + } + + /// + /// Gets the variance of the Weibull distribution. + /// + public double Variance + { + get + { + double mu = this.Mean; + return _scale * _scale * SpecialFunctions.Gamma(1.0 + 2.0 / _shape) - mu * mu; + } + } + + /// + /// Gets the standard deviation of the Weibull distribution. + /// + public double StdDev + { + get + { + return Math.Sqrt(this.Variance); + } + } + + /// + /// Gets the entropy of the Weibull distribution. + /// + public double Entropy + { + get + { + throw new NotImplementedException(); + } + } + + /// + /// Gets the skewness of the Weibull distribution. + /// + public double Skewness + { + get + { + double mu = this.Mean; + double sigma = this.StdDev; + double sigma2 = sigma * sigma; + double sigma3 = sigma2 * sigma; + return (_scale * _scale * _scale * SpecialFunctions.Gamma(1.0 + 3.0 / _shape) + - 3.0 * sigma2 * mu - mu * mu * mu) / sigma3; + } + } + #endregion + + #region IContinuousDistribution implementation + + /// + /// Gets the mode of the Weibull distribution. + /// + public double Mode + { + get + { + if (_shape > 1.0) + { + return _scale * Math.Pow((_shape - 1.0) / _shape, 1.0 / _shape); + } + else + { + return 0.0; + } + } + } + + /// + /// Gets the median of the Weibull distribution. + /// + public double Median + { + get + { + return _scale * Math.Pow(Constants.Ln2, 1.0 / _shape); + } + } + + /// + /// Gets the minimum of the Weibull distribution. + /// + public double Minimum + { + get { return 0.0; } + } + + /// + /// Gets the maximum of the Weibull distribution. + /// + public double Maximum + { + get { return Double.PositiveInfinity; } + } + + /// + /// Computes the density of the Weibull distribution. + /// + /// The location at which to compute the density. + /// the density at . + public double Density(double x) + { + if (x >= 0.0) + { + if (x == 0.0 && _shape == 1.0) + { + return _shape / _scale; + } + else + { + return _shape * Math.Pow(x / _scale, _shape - 1.0) * Math.Exp(-Math.Pow(x / _scale, _shape)) / _scale; + } + } + + return 0.0; + } + + /// + /// Computes the log density of the Weibull distribution. + /// + /// The location at which to compute the log density. + /// the log density at . + public double DensityLn(double x) + { + if (x >= 0.0) + { + if (x == 0.0 && _shape == 1.0) + { + return Math.Log(_shape) - Math.Log(_scale); + } + else + { + return Math.Log(_shape) + (_shape - 1.0) * Math.Log(x / _scale) - Math.Pow(x / _scale, _shape) - Math.Log(_scale); + } + } + + return double.NegativeInfinity; + } + + /// + /// Computes the cumulative distribution function of the Weibull distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + public double CumulativeDistribution(double x) + { + if (x >= 0.0) + { + return 1.0 - Math.Exp(-Math.Pow(x / _scale, _shape)); + } + + return 0.0; + } + + /// + /// Generates a sample from the Weibull distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return SampleWeibull(RandomSource, _shape, _scale); + } + + /// + /// Generates a sequence of samples from the Weibull distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return SampleWeibull(RandomSource, _shape, _scale); + } + } + + #endregion + + /// + /// Generates a sample from the Weibull distribution. + /// + /// The random number generator to use. + /// The shape of the Weibull distribution from which to generate samples. + /// The scale of the Weibull distribution from which to generate samples. + /// a sample from the distribution. + public static double Sample(Random rng, double shape, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return SampleWeibull(rng, shape, scale); + } + + /// + /// Generates a sequence of samples from the Weibull distribution. + /// + /// The random number generator to use. + /// The shape of the Weibull distribution from which to generate samples. + /// The scale of the Weibull distribution from which to generate samples. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(Random rng, double shape, double scale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return SampleWeibull(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(System.Random rnd, double shape, double scale) + { + double x = rnd.NextDouble(); + return scale * Math.Pow(-Math.Log(x), 1.0 / shape); + } + } +} diff --git a/src/Numerics/Distributions/Multivariate/Dirichlet.cs b/src/Numerics/Distributions/Multivariate/Dirichlet.cs index 3f7fa0de..8ecf7f88 100644 --- a/src/Numerics/Distributions/Multivariate/Dirichlet.cs +++ b/src/Numerics/Distributions/Multivariate/Dirichlet.cs @@ -257,7 +257,15 @@ namespace MathNet.Numerics.Distributions double sum = 0.0; for (int i = 0; i < n; i++) { - gv[i] = Gamma.Sample(rnd, alpha[i], 1.0); + if (alpha[i] == 0.0) + { + gv[i] = 0.0; + } + else + { + gv[i] = Gamma.Sample(rnd, alpha[i], 1.0); + } + sum += gv[i]; } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index bc812b2a..3a877efd 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -53,6 +53,7 @@ + diff --git a/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs b/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs new file mode 100644 index 00000000..5c56adef --- /dev/null +++ b/src/UnitTests/DistributionTests/Continuous/WeibullTests.cs @@ -0,0 +1,304 @@ +// +// 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 WeibullTests + { + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + } + + [Test, MultipleAsserts] + [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 CanCreateWeibull(double shape, double scale) + { + var n = new Weibull(shape, scale); + AssertEx.AreEqual(shape, n.Shape); + AssertEx.AreEqual(scale, n.Scale); + } + + [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)] + [Row(0.0, 0.0)] + [Row(0.0, 1.0)] + [Row(1.0, 0.0)] + public void WeibullCreateFailsWithBadParameters(double shape, double scale) + { + var n = new Weibull(shape, scale); + } + + [Test] + public void ValidateToString() + { + var n = new Weibull(1.0, 2.0); + AssertEx.AreEqual("Weibull(Shape = 1, Scale = 2)", n.ToString()); + } + + [Test] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetShape(double shape) + { + var n = new Weibull(1.0, 1.0); + n.Shape = shape; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(-0.0)] + [Row(0.0)] + [Row(-1.0)] + public void SetShapeFailsWithNegativeShape(double shape) + { + var n = new Weibull(1.0, 1.0); + n.Shape = shape; + } + + [Test] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void CanSetScale(double scale) + { + var n = new Weibull(1.0, 1.0); + n.Scale = scale; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(-0.0)] + [Row(0.0)] + [Row(-1.0)] + public void SetScaleFailsWithNegativeScale(double scale) + { + var n = new Weibull(1.0, 1.0); + n.Scale = scale; + } + + [Test] + [Row(1.0, 0.1, 0.1)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 9.5135076986687318362924871772654021925505786260884)] + [Row(10.0, 1.0, 0.95135076986687318362924871772654021925505786260884)] + public void ValidateMean(double shape, double scale, double mean) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(mean, n.Mean, 13); + } + + [Test] + [Row(1.0, 0.1, 0.01)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 1.3100455073468309147154581687505295026863354547057)] + [Row(10.0, 1.0, 0.013100455073468309147154581687505295026863354547057)] + public void ValidateVariance(double shape, double scale, double var) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(var, n.Variance, 13); + } + + [Test] + [Row(1.0, 0.1, 0.1)] + [Row(1.0, 1.0, 1.0)] + [Row(10.0, 10.0, 1.1445721940300799194124723631014002560036613065794)] + [Row(10.0, 1.0, 0.11445721940300799194124723631014002560036613065794)] + public void ValidateStdDev(double shape, double scale, double sdev) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(sdev, n.StdDev, 13); + } + + [Test] + [Row(1.0, 0.1, 2.0)] + [Row(1.0, 1.0, 2.0)] + [Row(10.0, 10.0, -0.63763713390314440916597757156663888653981696212127)] + [Row(10.0, 1.0, -0.63763713390314440916597757156663888653981696212127)] + public void ValidateSkewness(double shape, double scale, double skewness) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(skewness, n.Skewness, 11); + } + + [Test] + [Row(1.0, 0.1, 0.0)] + [Row(1.0, 1.0, 0.0)] + [Row(10.0, 10.0, 9.8951925820621439264623017041980483215553841533709)] + [Row(10.0, 1.0, 0.98951925820621439264623017041980483215553841533709)] + public void ValidateMode(double shape, double scale, double mode) + { + var n = new Weibull(shape, scale); + AssertEx.AreEqual(mode, n.Mode); + } + + [Test] + [Row(1.0, 0.1, 0.069314718055994530941723212145817656807550013436026)] + [Row(1.0, 1.0, 0.69314718055994530941723212145817656807550013436026)] + [Row(10.0, 10.0, 9.6401223546778973665856033763604752124634905617583)] + [Row(10.0, 1.0, 0.96401223546778973665856033763604752124634905617583)] + public void ValidateMedian(double shape, double scale, double median) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(median, n.Median, 13); + } + + [Test] + public void ValidateMinimum() + { + var n = new Weibull(1.0,1.0); + AssertEx.AreEqual(0.0, n.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var n = new Weibull(1.0, 1.0); + AssertEx.AreEqual(System.Double.PositiveInfinity, n.Maximum); + } + + [Test] + [Row(1.0, 0.1, 0.0, 10.0)] + [Row(1.0, 0.1, 1.0, 0.00045399929762484851535591515560550610237918088866565)] + [Row(1.0, 0.1, 10.0, 3.7200759760208359629596958038631183373588922923768e-43)] + [Row(1.0, 1.0, 0.0, 1.0)] + [Row(1.0, 1.0, 1.0, 0.36787944117144232159552377016146086744581113103177)] + [Row(1.0, 1.0, 10.0, 0.000045399929762484851535591515560550610237918088866565)] + [Row(10.0, 10.0, 0.0, 0.0)] + [Row(10.0, 10.0, 1.0, 9.9999999990000000000499999999983333333333750000000e-10)] + [Row(10.0, 10.0, 10.0, 0.36787944117144232159552377016146086744581113103177)] + [Row(10.0, 1.0, 0.0, 0.0)] + [Row(10.0, 1.0, 1.0, 3.6787944117144232159552377016146086744581113103177)] + [Row(10.0, 1.0, 10.0, 0.0)] + public void ValidateDensity(double shape, double scale, double x, double pdf) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(pdf, n.Density(x), 14); + } + + [Test] + [Row(1.0, 0.1, 0.0, 2.3025850929940456840179914546843642076011014886288)] + [Row(1.0, 0.1, 1.0, -7.6974149070059543159820085453156357923988985113712)] + [Row(1.0, 0.1, 10.0, -97.697414907005954315982008545315635792398898511371)] + [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, -20.723265837046411156161923092159277868409913397659)] + [Row(10.0, 10.0, 10.0, -1.0)] + [Row(10.0, 1.0, 0.0, Double.NegativeInfinity)] + [Row(10.0, 1.0, 1.0, 1.3025850929940456840179914546843642076011014886288)] + [Row(10.0, 1.0, 10.0, -9.999999976974149070059543159820085453156357923988985113712e9)] + public void ValidateDensityLn(double shape, double scale, double x, double pdfln) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14); + } + + [Test] + public void CanSampleStatic() + { + var d = Weibull.Sample(new Random(), 1.0, 1.0); + } + + [Test] + public void CanSampleSequenceStatic() + { + var ied = Weibull.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, Ignore("Catastrophic cancellation in one case. Fix this.")] + [Row(1.0, 0.1, 0.0, 0.0)] + [Row(1.0, 0.1, 1.0, 0.99995460007023751514846440848443944938976208191113)] + [Row(1.0, 0.1, 10.0, 0.99999999999999999999999999999999999999999996279924)] + [Row(1.0, 1.0, 0.0, 0.0)] + [Row(1.0, 1.0, 1.0, 0.63212055882855767840447622983853913255418886896823)] + [Row(1.0, 1.0, 10.0, 0.99995460007023751514846440848443944938976208191113)] + [Row(10.0, 10.0, 0.0, 0.0)] + [Row(10.0, 10.0, 1.0, 9.9999999995000000000166666666662500000000083333333e-11)] + [Row(10.0, 10.0, 10.0, 0.63212055882855767840447622983853913255418886896823)] + [Row(10.0, 1.0, 0.0, 0.0)] + [Row(10.0, 1.0, 1.0, 0.63212055882855767840447622983853913255418886896823)] + [Row(10.0, 1.0, 10.0, 1.0)] + public void ValidateCumulativeDistribution(double shape, double scale, double x, double cdf) + { + var n = new Weibull(shape, scale); + AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 15); + } + } +} diff --git a/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs b/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs index c31899a4..17943614 100644 --- a/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs @@ -165,10 +165,17 @@ namespace MathNet.Numerics.UnitTests.DistributionTests } [Test] - public void Sample() + public void CanSampleSymmetricDirichlet() { Dirichlet d = new Dirichlet(1.0, 5); double[] s = d.Sample(); } + + [Test] + public void CanSampleSingularDirichlet() + { + Dirichlet d = new Dirichlet(new double[] {2.0, 1.0, 0.0, 3.0}); + double[] s = d.Sample(); + } } } \ No newline at end of file diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index b9d2d67b..f6cd4957 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -67,6 +67,7 @@ +