diff --git a/src/Examples/ContinuousDistributions/TriangularDistribution.cs b/src/Examples/ContinuousDistributions/TriangularDistribution.cs new file mode 100644 index 00000000..bd4c78a0 --- /dev/null +++ b/src/Examples/ContinuousDistributions/TriangularDistribution.cs @@ -0,0 +1,145 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 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. +// + +using System; +using MathNet.Numerics.Distributions; + +namespace Examples.TriangularExamples +{ + /// + /// ContinuousUniform distribution example + /// + public class TriangularDistribution : IExample + { + /// + /// Gets the name of this example + /// + /// + public string Name + { + get + { + return "Triangular distribution"; + } + } + + /// + /// Gets the description of this example + /// + public string Description + { + get + { + return "Triangular distribution properties and samples generating examples"; + } + } + + /// + /// Run example + /// + /// Triangular distribution + public void Run() + { + // 1. Initialize + var triangular = new Triangular(0, 1, 0.3); + Console.WriteLine(@"1. Initialize the new instance of the Triangular distribution class with parameters Lower = {0}, Upper = {1}, Mode = {2}", triangular.LowerBound, triangular.UpperBound, triangular.Mode); + Console.WriteLine(); + + // 2. Distributuion properties: + Console.WriteLine(@"2. {0} distributuion properties:", triangular); + + // Cumulative distribution function + Console.WriteLine(@"{0} - Сumulative distribution at location '0.3'", triangular.CumulativeDistribution(0.3).ToString(" #0.00000;-#0.00000")); + + // Probability density + Console.WriteLine(@"{0} - Probability density at location '0.3'", triangular.Density(0.3).ToString(" #0.00000;-#0.00000")); + + // Log probability density + Console.WriteLine(@"{0} - Log probability density at location '0.3'", triangular.DensityLn(0.3).ToString(" #0.00000;-#0.00000")); + + // Entropy + Console.WriteLine(@"{0} - Entropy", triangular.Entropy.ToString(" #0.00000;-#0.00000")); + + // Largest element in the domain + Console.WriteLine(@"{0} - Largest element in the domain", triangular.Maximum.ToString(" #0.00000;-#0.00000")); + + // Smallest element in the domain + Console.WriteLine(@"{0} - Smallest element in the domain", triangular.Minimum.ToString(" #0.00000;-#0.00000")); + + // Mean + Console.WriteLine(@"{0} - Mean", triangular.Mean.ToString(" #0.00000;-#0.00000")); + + // Median + Console.WriteLine(@"{0} - Median", triangular.Median.ToString(" #0.00000;-#0.00000")); + + // Mode + Console.WriteLine(@"{0} - Mode", triangular.Mode.ToString(" #0.00000;-#0.00000")); + + // Variance + Console.WriteLine(@"{0} - Variance", triangular.Variance.ToString(" #0.00000;-#0.00000")); + + // Standard deviation + Console.WriteLine(@"{0} - Standard deviation", triangular.StdDev.ToString(" #0.00000;-#0.00000")); + + // Skewness + Console.WriteLine(@"{0} - Skewness", triangular.Skewness.ToString(" #0.00000;-#0.00000")); + Console.WriteLine(); + + // 10 samples + Console.WriteLine(@"3. Generate 10 samples of the Triangular distribution"); + for (var i = 0; i < 10; i++) + { + Console.Write(triangular.Sample().ToString("N05") + @" "); + } + + Console.WriteLine(); + Console.WriteLine(); + + // 10000 samples with starting parameters + Console.WriteLine(@"4. Generate 100000 samples of the Triangular({0}, {1}, {2}) distribution and display histogram", triangular.LowerBound, triangular.UpperBound, triangular.Mode); + var data = new double[100000]; + for (var i = 0; i < data.Length; i++) + { + data[i] = triangular.Sample(); + } + + ConsoleHelper.DisplayHistogram(data); + Console.WriteLine(); + + // 10000 with different parameters + triangular.UpperBound = 10; + triangular.Mode = 8; + triangular.LowerBound = 2; + Console.WriteLine(@"4. Generate 100000 samples of the Triangular({0}, {1}, {2}) distribution and display histogram", triangular.LowerBound, triangular.UpperBound, triangular.Mode); + for (var i = 0; i < data.Length; i++) + { + data[i] = triangular.Sample(); + } + + ConsoleHelper.DisplayHistogram(data); + } + } +} diff --git a/src/Examples/Examples.csproj b/src/Examples/Examples.csproj index 3b594d8b..376d6fdd 100644 --- a/src/Examples/Examples.csproj +++ b/src/Examples/Examples.csproj @@ -88,6 +88,7 @@ + diff --git a/src/Numerics/Distributions/Triangular.cs b/src/Numerics/Distributions/Triangular.cs new file mode 100644 index 00000000..4fea433b --- /dev/null +++ b/src/Numerics/Distributions/Triangular.cs @@ -0,0 +1,420 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 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. +// + +using System; +using System.Collections.Generic; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics.Distributions +{ + /// + /// Triangular distribution. + /// For details, see Wikipedia - Triangular 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 whether all the incoming parameters 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 Triangular : IContinuousDistribution + { + System.Random _random; + + double _lower; + double _upper; + double _mode; + + /// + /// Initializes a new instance of the Triangular class with the given lower bound, upper bound and mode. + /// + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// If the upper bound is smaller than the mode or if the mode is smaller than the lower bound. + public Triangular(double lower, double upper, double mode) + { + _random = MersenneTwister.Default; + SetParameters(lower, upper, mode); + } + + /// + /// Initializes a new instance of the Triangular class with the given lower bound, upper bound and mode. + /// + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// The random number generator which is used to draw random samples. + /// If the upper bound is smaller than the mode or if the mode is smaller than the lower bound. + public Triangular(double lower, double upper, double mode, System.Random randomSource) + { + _random = randomSource ?? MersenneTwister.Default; + SetParameters(lower, upper, mode); + } + + /// + /// A string representation of the distribution. + /// + /// a string representation of the distribution. + public override string ToString() + { + return "Triangular(Lower = " + _lower + ", Upper = " + _upper + ", Mode = " + _mode + ")"; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// When the parameters are out of range. + void SetParameters(double lower, double upper, double mode) + { + if (upper < mode || mode < lower || Double.IsNaN(upper) || Double.IsNaN(lower) || Double.IsNaN(mode)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _lower = lower; + _upper = upper; + _mode = mode; + } + + /// + /// Gets or sets the lower bound of the distribution. + /// + public double LowerBound + { + get { return _lower; } + set { SetParameters(value, _upper, _mode); } + } + + /// + /// Gets or sets the upper bound of the distribution. + /// + public double UpperBound + { + get { return _upper; } + set { SetParameters(_lower, value, _mode); } + } + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public System.Random RandomSource + { + get { return _random; } + set { _random = value ?? MersenneTwister.Default; } + } + + /// + /// Gets the mean of the distribution. + /// + public double Mean + { + get { return (_lower + _upper + _mode) / 3.0; } + } + + /// + /// Gets the variance of the distribution. + /// + public double Variance + { + get + { + var a = _lower; + var b = _upper; + var c = _mode; + return (a * a + b * b + c * c - a * b - a * c - b * c) / 18.0; + } + } + + /// + /// Gets the standard deviation of the distribution. + /// + public double StdDev + { + get { return Math.Sqrt(Variance); } + } + + /// + /// Gets the entropy of the distribution. + /// + /// + public double Entropy + { + get { return 0.5 + Math.Log((_upper - _lower) / 2); } + } + + /// + /// Gets the skewness of the distribution. + /// + public double Skewness + { + get + { + var a = _lower; + var b = _upper; + var c = _mode; + var q = Math.Sqrt(2) * (a + b - 2 * c) * (2 * a - b - c) * (a - 2 * b + c); + var d = 5 * Math.Pow(a * a + b * b + c * c - a * b - a * c - b * c, 3.0 / 2); + return q / d; + } + } + + /// + /// Gets or sets the mode of the distribution. + /// + public double Mode + { + get { return _mode; } + set { SetParameters(_lower, _upper, value); } + } + + /// + /// Gets the median of the distribution. + /// + /// + public double Median + { + get + { + var a = _lower; + var b = _upper; + var c = _mode; + return c >= (a + b) / 2 + ? a + Math.Sqrt((b - a) * (c - a) / 2) + : b - Math.Sqrt((b - a) * (b - c) / 2); + } + } + + /// + /// Gets the minimum of the distribution. + /// + public double Minimum + { + get { return _lower; } + } + + /// + /// Gets the maximum of the distribution. + /// + public double Maximum + { + get { return _upper; } + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location at which to compute the density. + /// the density at . + /// + public double Density(double x) + { + return PDF(_lower, _upper, _mode, x); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location at which to compute the log density. + /// the log density at . + /// + public double DensityLn(double x) + { + return PDFLn(_lower, _upper, _mode, x); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public double CumulativeDistribution(double x) + { + return CDF(_lower, _upper, _mode, x); + } + + /// + /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution + /// at the given probability. This is also known as the quantile or percent point function. + /// + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . + /// + public double InverseCumulativeDistribution(double p) + { + return InvCDF(_lower, _upper, _mode, p); + } + + /// + /// Generates a sample from the Triangular distribution. + /// + /// a sample from the distribution. + public double Sample() + { + return Sample(_random, _lower, _upper, _mode); + } + + /// + /// Generates a sequence of samples from the Triangular distribution. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return Samples(_random, _lower, _upper, _mode); + } + + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double lower, double upper, double mode, double x) + { + if (upper < mode) throw new ArgumentOutOfRangeException("upper", Resources.InvalidDistributionParameters); + if (mode < lower) throw new ArgumentOutOfRangeException("lower", Resources.InvalidDistributionParameters); // TODO: Is "lower" the appropriate argument here? + + var a = lower; + var b = upper; + var c = mode; + + if (a <= x && x <= c) return 2 * (x - a) / ((b - a) * (c - a)); + if (c < x & x <= b) return 2 * (b - x) / ((b - a) * (b - c)); + return 0; + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double lower, double upper, double mode, double x) + { + return Math.Log(PDF(lower, upper, mode, x)); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// the cumulative distribution at location . + /// + public static double CDF(double lower, double upper, double mode, double x) + { + if (upper < mode) throw new ArgumentOutOfRangeException("upper", Resources.InvalidDistributionParameters); + if (mode < lower) throw new ArgumentOutOfRangeException("lower", Resources.InvalidDistributionParameters); // TODO: Is "lower" the appropriate argument here? + + var a = lower; + var b = upper; + var c = mode; + + if (x < a) return 0; + if (a <= x && x <= c) return (x - a) * (x - a) / ((b - a) * (c - a)); + if (c < x & x <= b) return 1 - (b - x) * (b - x) / ((b - a) * (b - c)); + return 1; + } + + /// + /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution + /// at the given probability. This is also known as the quantile or percent point function. + /// + /// The location at which to compute the inverse cumulative density. + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// the inverse cumulative density at . + /// + public static double InvCDF(double lower, double upper, double mode, double p) + { + if (upper < mode) throw new ArgumentOutOfRangeException("upper", Resources.InvalidDistributionParameters); + if (mode < lower) throw new ArgumentOutOfRangeException("lower", Resources.InvalidDistributionParameters); // TODO: Is "lower" the appropriate argument here? + + var a = lower; + var b = upper; + var c = mode; + + if (p <= 0) return lower; + // Taken from http://www.ntrand.com/triangular-distribution/ + if (p < (c - a) / (b - a)) return a + Math.Sqrt(p * (c - a) * (b - a)); + if (p < 1) return b - Math.Sqrt((1 - p) * (b - c) * (b - a)); + return upper; + } + + /// + /// Generates a sample from the Triangular distribution. + /// + /// The random number generator to use. + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// a sample from the distribution. + public double Sample(System.Random rnd, double lower, double upper, double mode) + { + var a = lower; + var b = upper; + var c = mode; + var u = rnd.NextDouble(); + + return u < (c - a) / (b - a) + ? a + Math.Sqrt(u * (b - a) * (c - a)) + : b - Math.Sqrt((1 - u) * (b - a) * (b - c)); ; + } + + /// + /// Generates a sequence of samples from the Triangular distribution. + /// + /// The random number generator to use. + /// Lower bound. Range: lower ≤ mode ≤ upper + /// Upper bound. Range: lower ≤ mode ≤ upper + /// Mode (most frequent value). Range: lower ≤ mode ≤ upper + /// a sequence of samples from the distribution. + public IEnumerable Samples(System.Random rnd, double lower, double upper, double mode) + { + while (true) + { + yield return Sample(rnd, lower, upper, mode); + } + } + + } +} \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index f253ffc0..862e0bb9 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -87,6 +87,7 @@ + diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index 91380f5c..8fcb6862 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/src/UnitTests/DistributionTests/CommonDistributionTests.cs @@ -122,6 +122,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests new Normal(0.0, 1.0), new Weibull(1.0, 1.0), new LogNormal(1.0, 1.0), + new Triangular(0, 1, 0.7), new StudentT(0.0, 1.0, 5.0) }; }