diff --git a/src/Numerics/Distributions/TruncatedNormal.cs b/src/Numerics/Distributions/TruncatedNormal.cs index b8419593..72b33b58 100644 --- a/src/Numerics/Distributions/TruncatedNormal.cs +++ b/src/Numerics/Distributions/TruncatedNormal.cs @@ -72,9 +72,9 @@ namespace MathNet.Numerics.Distributions { readonly double _cumulativeDensityWithinBounds; /// - /// 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 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 random number generator. The mean + /// and standard deviation are that of the untruncated normal distribution. /// /// The mean (μ) of the untruncated distribution. /// The standard deviation (σ) of the untruncated distribution. Range: σ > 0. @@ -88,22 +88,26 @@ namespace MathNet.Numerics.Distributions { } /// - /// Initializes a new instance of the Normal class with a particular mean and standard deviation. The distribution will - /// be initialized with the default random number generator. + /// Initializes a new instance of the TruncatedNormal class. The distribution will + /// be initialized with the provided random number generator. /// - /// The mean (μ) of the normal distribution. - /// The standard deviation (σ) of the normal distribution. Range: σ > 0. + /// The mean (μ) of the untruncated normal distribution. + /// The standard deviation (σ) of the untruncated normal distribution. Range: σ > 0. /// The random number generator which is used to draw random samples. - public TruncatedNormal(double mean, double stddev, System.Random randomSource, double lowerBound = double.NegativeInfinity, double upperBound = double.PositiveInfinity) + /// The inclusive lower bound of the truncated distribution. Default is double.NegativeInfinity. + /// The inclusive upper bound of the truncated distribution. Must be larger than . + /// Default is double.PositiveInfinity. + + 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 Samples() { - throw new NotImplementedException(); + while (true) { + yield return Sample(); + } } /// diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index e0991dd6..82fcb5f4 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/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] diff --git a/src/UnitTests/DistributionTests/Continuous/TruncatedNormalTests.cs b/src/UnitTests/DistributionTests/Continuous/TruncatedNormalTests.cs index 7605c165..faf48cb6 100644 --- a/src/UnitTests/DistributionTests/Continuous/TruncatedNormalTests.cs +++ b/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); } /// @@ -156,6 +156,61 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous } } + /// + /// Validate density when only one bound is are specified. + /// + /// Mean value. + /// Standard deviation value. + [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); + } + } + + /// + /// Validate density when both bounds are specified. + /// + /// Mean value. + /// Standard deviation value. + [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); + } + } + + /// /// Validate density log when no bounds are specified. Uses same /// test cases as the Normal distribution as should be equivalent in this case.