Browse Source

Merge pull request #184 from Superbest/triangular-distribution

Triangular distribution
optimization-3
Christoph Ruegg 13 years ago
parent
commit
566ec8d695
  1. 145
      src/Examples/ContinuousDistributions/TriangularDistribution.cs
  2. 1
      src/Examples/Examples.csproj
  3. 420
      src/Numerics/Distributions/Triangular.cs
  4. 1
      src/Numerics/Numerics.csproj
  5. 1
      src/UnitTests/DistributionTests/CommonDistributionTests.cs

145
src/Examples/ContinuousDistributions/TriangularDistribution.cs

@ -0,0 +1,145 @@
// <copyright file="ContinuousUniformDistribution.cs" company="Math.NET">
// 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.
// </copyright>
using System;
using MathNet.Numerics.Distributions;
namespace Examples.TriangularExamples
{
/// <summary>
/// ContinuousUniform distribution example
/// </summary>
public class TriangularDistribution : IExample
{
/// <summary>
/// Gets the name of this example
/// </summary>
/// <seealso cref="http://reference.wolfram.com/mathematica/ref/TriangularDistribution.html"/>
public string Name
{
get
{
return "Triangular distribution";
}
}
/// <summary>
/// Gets the description of this example
/// </summary>
public string Description
{
get
{
return "Triangular distribution properties and samples generating examples";
}
}
/// <summary>
/// Run example
/// </summary>
/// <a href="https://en.wikipedia.org/wiki/Triangular_distribution">Triangular distribution</a>
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);
}
}
}

1
src/Examples/Examples.csproj

@ -88,6 +88,7 @@
<Compile Include="ContinuousDistributions\CauchyDistribution.cs" /> <Compile Include="ContinuousDistributions\CauchyDistribution.cs" />
<Compile Include="ContinuousDistributions\ChiDistribution.cs" /> <Compile Include="ContinuousDistributions\ChiDistribution.cs" />
<Compile Include="ContinuousDistributions\ChiSquareDistribution.cs" /> <Compile Include="ContinuousDistributions\ChiSquareDistribution.cs" />
<Compile Include="ContinuousDistributions\TriangularDistribution.cs" />
<Compile Include="ContinuousDistributions\ContinuousUniformDistribution.cs" /> <Compile Include="ContinuousDistributions\ContinuousUniformDistribution.cs" />
<Compile Include="ContinuousDistributions\WeibullDistribution.cs" /> <Compile Include="ContinuousDistributions\WeibullDistribution.cs" />
<Compile Include="ContinuousDistributions\StudentTDistribution.cs" /> <Compile Include="ContinuousDistributions\StudentTDistribution.cs" />

420
src/Numerics/Distributions/Triangular.cs

@ -0,0 +1,420 @@
// <copyright file="TriangularUniform.cs" company="Math.NET">
// 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.
// </copyright>
using System;
using System.Collections.Generic;
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
namespace MathNet.Numerics.Distributions
{
/// <summary>
/// Triangular distribution.
/// For details, see <a href="https://en.wikipedia.org/wiki/Triangular_distribution">Wikipedia - Triangular distribution</a>.
/// </summary>
/// <remarks><para>The distribution will use the <see cref="System.Random"/> by default.
/// Users can get/set the random number generator by using the <see cref="RandomSource"/> property.</para>
/// <para>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 <c>false</c>, all parameter checks can be turned off.</para></remarks>
public class Triangular : IContinuousDistribution
{
System.Random _random;
double _lower;
double _upper;
double _mode;
/// <summary>
/// Initializes a new instance of the Triangular class with the given lower bound, upper bound and mode.
/// </summary>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <exception cref="ArgumentException">If the upper bound is smaller than the mode or if the mode is smaller than the lower bound.</exception>
public Triangular(double lower, double upper, double mode)
{
_random = MersenneTwister.Default;
SetParameters(lower, upper, mode);
}
/// <summary>
/// Initializes a new instance of the Triangular class with the given lower bound, upper bound and mode.
/// </summary>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <param name="randomSource">The random number generator which is used to draw random samples.</param>
/// <exception cref="ArgumentException">If the upper bound is smaller than the mode or if the mode is smaller than the lower bound.</exception>
public Triangular(double lower, double upper, double mode, System.Random randomSource)
{
_random = randomSource ?? MersenneTwister.Default;
SetParameters(lower, upper, mode);
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return "Triangular(Lower = " + _lower + ", Upper = " + _upper + ", Mode = " + _mode + ")";
}
/// <summary>
/// Sets the parameters of the distribution after checking their validity.
/// </summary>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <exception cref="ArgumentOutOfRangeException">When the parameters are out of range.</exception>
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;
}
/// <summary>
/// Gets or sets the lower bound of the distribution.
/// </summary>
public double LowerBound
{
get { return _lower; }
set { SetParameters(value, _upper, _mode); }
}
/// <summary>
/// Gets or sets the upper bound of the distribution.
/// </summary>
public double UpperBound
{
get { return _upper; }
set { SetParameters(_lower, value, _mode); }
}
/// <summary>
/// Gets or sets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? MersenneTwister.Default; }
}
/// <summary>
/// Gets the mean of the distribution.
/// </summary>
public double Mean
{
get { return (_lower + _upper + _mode) / 3.0; }
}
/// <summary>
/// Gets the variance of the distribution.
/// </summary>
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;
}
}
/// <summary>
/// Gets the standard deviation of the distribution.
/// </summary>
public double StdDev
{
get { return Math.Sqrt(Variance); }
}
/// <summary>
/// Gets the entropy of the distribution.
/// </summary>
/// <value></value>
public double Entropy
{
get { return 0.5 + Math.Log((_upper - _lower) / 2); }
}
/// <summary>
/// Gets the skewness of the distribution.
/// </summary>
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;
}
}
/// <summary>
/// Gets or sets the mode of the distribution.
/// </summary>
public double Mode
{
get { return _mode; }
set { SetParameters(_lower, _upper, value); }
}
/// <summary>
/// Gets the median of the distribution.
/// </summary>
/// <value></value>
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);
}
}
/// <summary>
/// Gets the minimum of the distribution.
/// </summary>
public double Minimum
{
get { return _lower; }
}
/// <summary>
/// Gets the maximum of the distribution.
/// </summary>
public double Maximum
{
get { return _upper; }
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return PDF(_lower, _upper, _mode, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return PDFLn(_lower, _upper, _mode, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return CDF(_lower, _upper, _mode, x);
}
/// <summary>
/// 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.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative density.</param>
/// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
/// <seealso cref="InvCDF"/>
public double InverseCumulativeDistribution(double p)
{
return InvCDF(_lower, _upper, _mode, p);
}
/// <summary>
/// Generates a sample from the <c>Triangular</c> distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return Sample(_random, _lower, _upper, _mode);
}
/// <summary>
/// Generates a sequence of samples from the <c>Triangular</c> distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
return Samples(_random, _lower, _upper, _mode);
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
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;
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double lower, double upper, double mode, double x)
{
return Math.Log(PDF(lower, upper, mode, x));
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
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;
}
/// <summary>
/// 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.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative density.</param>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
/// <seealso cref="InverseCumulativeDistribution"/>
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;
}
/// <summary>
/// Generates a sample from the <c>Triangular</c> distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <returns>a sample from the distribution.</returns>
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)); ;
}
/// <summary>
/// Generates a sequence of samples from the <c>Triangular</c> distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="lower">Lower bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="upper">Upper bound. Range: lower ≤ mode ≤ upper</param>
/// <param name="mode">Mode (most frequent value). Range: lower ≤ mode ≤ upper</param>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples(System.Random rnd, double lower, double upper, double mode)
{
while (true)
{
yield return Sample(rnd, lower, upper, mode);
}
}
}
}

1
src/Numerics/Numerics.csproj

@ -87,6 +87,7 @@
<Reference Include="System.Xml" /> <Reference Include="System.Xml" />
</ItemGroup> </ItemGroup>
<ItemGroup> <ItemGroup>
<Compile Include="Distributions\Triangular.cs" />
<Compile Include="Euclid.cs" /> <Compile Include="Euclid.cs" />
<Compile Include="Generate.cs" /> <Compile Include="Generate.cs" />
<Compile Include="GoodnessOfFit.cs" /> <Compile Include="GoodnessOfFit.cs" />

1
src/UnitTests/DistributionTests/CommonDistributionTests.cs

@ -122,6 +122,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
new Normal(0.0, 1.0), new Normal(0.0, 1.0),
new Weibull(1.0, 1.0), new Weibull(1.0, 1.0),
new LogNormal(1.0, 1.0), new LogNormal(1.0, 1.0),
new Triangular(0, 1, 0.7),
new StudentT(0.0, 1.0, 5.0) new StudentT(0.0, 1.0, 5.0)
}; };
} }

Loading…
Cancel
Save