Browse Source

Simple sampling implementation

Also corrected some comments.
Added more tests
truncatednormal
BenHewins 11 years ago
parent
commit
656f235e31
  1. 40
      src/Numerics/Distributions/TruncatedNormal.cs
  2. 4
      src/UnitTests/DistributionTests/CommonDistributionTests.cs
  3. 57
      src/UnitTests/DistributionTests/Continuous/TruncatedNormalTests.cs

40
src/Numerics/Distributions/TruncatedNormal.cs

@ -72,9 +72,9 @@ namespace MathNet.Numerics.Distributions {
readonly double _cumulativeDensityWithinBounds;
/// <summary>
/// Initializes a new instance of the TruncatedNormal class with a particular mean, standard deviation, lower bound, and upper bound. The distribution will
/// be initialized with the default <seealso cref="System.Random"/> random number generator. The mean and standard deviation are that of the untruncated
/// normal distribution.
/// Initializes a new instance of the TruncatedNormal class. The distribution will
/// be initialized with the default <seealso cref="System.Random"/> random number generator. The mean
/// and standard deviation are that of the untruncated normal distribution.
/// </summary>
/// <param name="mean">The mean (μ) of the untruncated distribution.</param>
/// <param name="stddev">The standard deviation (σ) of the untruncated distribution. Range: σ > 0.</param>
@ -88,22 +88,26 @@ namespace MathNet.Numerics.Distributions {
}
/// <summary>
/// Initializes a new instance of the Normal class with a particular mean and standard deviation. The distribution will
/// be initialized with the default <seealso cref="System.Random"/> random number generator.
/// Initializes a new instance of the TruncatedNormal class. The distribution will
/// be initialized with the provided <seealso cref="System.Random"/> random number generator.
/// </summary>
/// <param name="mean">The mean (μ) of the normal distribution.</param>
/// <param name="stddev">The standard deviation (σ) of the normal distribution. Range: σ > 0.</param>
/// <param name="untruncatedMean">The mean (μ) of the untruncated normal distribution.</param>
/// <param name="untruncatedStdDev">The standard deviation (σ) of the untruncated normal distribution. Range: σ > 0.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples.</param>
public TruncatedNormal(double mean, double stddev, System.Random randomSource, double lowerBound = double.NegativeInfinity, double upperBound = double.PositiveInfinity)
/// <param name="lowerBound">The inclusive lower bound of the truncated distribution. Default is double.NegativeInfinity.</param>
/// <param name="upperBound">The inclusive upper bound of the truncated distribution. Must be larger than <paramref name="lowerBound"/>.
/// Default is double.PositiveInfinity.</param>
public TruncatedNormal(double untruncatedMean, double untruncatedStdDev, System.Random randomSource, double lowerBound = double.NegativeInfinity, double upperBound = double.PositiveInfinity)
{
if (!IsValidParameterSet(mean, stddev, lowerBound, upperBound))
if (!IsValidParameterSet(untruncatedMean, untruncatedStdDev, lowerBound, upperBound))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
_random = randomSource ?? SystemRandomSource.Default;
_mu = mean;
_sigma = stddev;
_mu = untruncatedMean;
_sigma = untruncatedStdDev;
_lowerBound = lowerBound;
_upperBound = upperBound;
_alpha = (_lowerBound - _mu) / _sigma;
@ -279,22 +283,26 @@ namespace MathNet.Numerics.Distributions {
return _standardNormal.DensityLn((x - _mu) / _sigma) - Math.Log(_sigma) - Math.Log(_cumulativeDensityWithinBounds);
}
//TODO: implement sampling, use method described by Mazet here: http://miv.u-strasbg.fr/mazet/rtnorm/
// see implementations listed on that page for examples.
public double Sample()
{
throw new NotImplementedException();
//TODO: implement sampling more efficiently/accurately, use method described by Mazet here: http://miv.u-strasbg.fr/mazet/rtnorm/
// see implementations listed on that page for examples.
return InverseCumulativeDistribution(RandomSource.NextDouble());
}
public void Samples(double[] values)
{
throw new NotImplementedException();
for(int i = 0; i < values.Length; i++) {
values[i] = Sample();
}
}
public IEnumerable<double> Samples()
{
throw new NotImplementedException();
while (true) {
yield return Sample();
}
}
/// <summary>

4
src/UnitTests/DistributionTests/CommonDistributionTests.cs

@ -89,8 +89,8 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
new StudentT(0.0, 1.0, 5.0),
new Triangular(0, 1, 0.7),
new Weibull(1.0, 1.0),
new TruncatedNormal(0, 1.0, -5.0, 5.0), //Finite
new TruncatedNormal(0, 1.0, -5.0), //Semi-finite
new TruncatedNormal(0, 1.0, -1.0, 1.5), //Finite
new TruncatedNormal(0, 1.0, -0.5), //Semi-finite
};
[Test]

57
src/UnitTests/DistributionTests/Continuous/TruncatedNormalTests.cs

@ -110,7 +110,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
[TestCase(Double.PositiveInfinity, 1.0)]
public void ValidateCumulativeNoBounds(double x, double p) {
var truncatedNormal = new TruncatedNormal(5.0, 2.0);
AssertHelpers.AlmostEqualRelative(p, truncatedNormal.CumulativeDistribution(x), 14);
AssertHelpers.AlmostEqualRelative(p, truncatedNormal.CumulativeDistribution(x), 14);
}
/// <summary>
@ -156,6 +156,61 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
}
}
/// <summary>
/// Validate density when only one bound is are specified.
/// </summary>
/// <param name="mean">Mean value.</param>
/// <param name="sdev">Standard deviation value.</param>
[TestCase(10.0, 0.1, -5.0)]
[TestCase(-5.0, 1.0, 3.0)]
[TestCase(0.0, 10.0, -10.0)]
[TestCase(10.0, 100.0, 15.0)]
[TestCase(-5.0, Double.PositiveInfinity, -5.0)]
public void ValidateDensitySemiFinite(double mean, double sdev, double lowerBound) {
var truncatedNormal = new TruncatedNormal(mean, sdev, lowerBound);
var normal = new Normal(mean, sdev);
for (var i = 0; i < 11; i++) {
var x = i - 5.0;
double density;
if(x < lowerBound) {
density = 0d;
} else {
var d = (mean - x) / sdev;
var pdf = Math.Exp(-0.5 * d * d) / (sdev * Constants.Sqrt2Pi);
density = pdf / (1.0 - normal.CumulativeDistribution(lowerBound));
}
AssertHelpers.AlmostEqualRelative(density, truncatedNormal.Density(x), 14);
}
}
/// <summary>
/// Validate density when both bounds are specified.
/// </summary>
/// <param name="mean">Mean value.</param>
/// <param name="sdev">Standard deviation value.</param>
[TestCase(10.0, 0.1, -5.0, 5.0)]
[TestCase(-5.0, 1.0, double.NegativeInfinity, -5.0)]
[TestCase(0.0, 10.0, -10.0, 15.0)]
[TestCase(10.0, 100.0, 15.0, 100.0)]
[TestCase(-5.0, Double.PositiveInfinity, -5.0, 0.0)]
public void ValidateDensityFinite(double mean, double sdev, double lowerBound, double upperBound) {
var truncatedNormal = new TruncatedNormal(mean, sdev, lowerBound, upperBound);
var normal = new Normal(mean, sdev);
for (var i = 0; i < 11; i++) {
var x = i - 5.0;
double density;
if (x < lowerBound || upperBound < x) {
density = 0d;
} else {
var d = (mean - x) / sdev;
var pdf = Math.Exp(-0.5 * d * d) / (sdev * Constants.Sqrt2Pi);
density = pdf / (normal.CumulativeDistribution(upperBound) - normal.CumulativeDistribution(lowerBound));
}
AssertHelpers.AlmostEqualRelative(density, truncatedNormal.Density(x), 14);
}
}
/// <summary>
/// Validate density log when no bounds are specified. Uses same
/// test cases as the Normal distribution as should be equivalent in this case.

Loading…
Cancel
Save