From 05269f5c7a2228ca90a437ef22dc2ff7ec27e5d1 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Thu, 9 Jul 2009 06:21:15 +0800 Subject: [PATCH] Continued implementation of the normal distribution. Signed-off-by: jvangael --- .../Continuous/NormalTests.cs | 110 +++++++++++ .../Distributions/Continuous/Normal.cs | 175 ++++++++++++++++-- src/Managed/Distributions/IDistribution.cs | 2 +- 3 files changed, 267 insertions(+), 20 deletions(-) diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs index 24bdbc97..c825d36a 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs @@ -215,5 +215,115 @@ namespace MathNet.Numerics.UnitTests var n = new Normal(); n.Mean = mean; } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void ValidateEntropy(double sdev) + { + var n = new Normal(1.0, sdev); + AssertEx.AreEqual(MathNet.Numerics.Constants.LogSqrt2PiE + Math.Log(n.StdDev), n.Entropy); + } + + [Test] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void ValidateSkewness(double sdev) + { + var n = new Normal(1.0, sdev); + AssertEx.AreEqual(0.0, n.Skewness); + } + + [Test] + [Row(Double.NegativeInfinity)] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void ValidateMode(double mean) + { + var n = new Normal(mean, 1.0); + AssertEx.AreEqual(mean, n.Mode); + } + + [Test] + [Row(Double.NegativeInfinity)] + [Row(-0.0)] + [Row(0.0)] + [Row(0.1)] + [Row(1.0)] + [Row(10.0)] + [Row(Double.PositiveInfinity)] + public void ValidateMedian(double mean) + { + var n = new Normal(mean, 1.0); + AssertEx.AreEqual(mean, n.Median); + } + + [Test] + public void ValidateMinimum() + { + var n = new Normal(); + AssertEx.AreEqual(System.Double.NegativeInfinity, n.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var n = new Normal(); + AssertEx.AreEqual(System.Double.PositiveInfinity, n.Maximum); + } + + [Test] + [Row(0.0, 0.0)] + [Row(0.0, 0.1)] + [Row(0.0, 1.0)] + [Row(0.0, 10.0)] + [Row(10.0, 1.0)] + [Row(-5.0, 100.0)] + [Row(0.0, Double.PositiveInfinity)] + public void ValidateDensity(double mean, double sdev) + { + var n = Normal.WithMeanStdDev(mean, sdev); + for(int i = 0; i < 11; i++) + { + double x = i - 5.0; + double d = (mean - x)/sdev; + double pdf = Math.Exp(-0.5*d*d)/(sdev*Constants.Sqrt2Pi); + AssertEx.AreEqual(pdf, n.Density(x)); + } + } + + [Test] + [Row(0.0, 0.0)] + [Row(0.0, 0.1)] + [Row(0.0, 1.0)] + [Row(0.0, 10.0)] + [Row(10.0, 1.0)] + [Row(-5.0, 100.0)] + [Row(0.0, Double.PositiveInfinity)] + public void ValidateDensityLn(double mean, double sdev) + { + var n = Normal.WithMeanStdDev(mean, sdev); + for (int i = 0; i < 11; i++) + { + double x = i - 5.0; + double d = (mean - x) / sdev; + double pdfln = -0.5 * d * d - Math.Log(sdev) - Constants.LogSqrt2Pi; + AssertEx.AreEqual(pdfln, n.DensityLn(x)); + } + } + + test samplers } } diff --git a/src/Managed/Distributions/Continuous/Normal.cs b/src/Managed/Distributions/Continuous/Normal.cs index 3b434838..8110a000 100644 --- a/src/Managed/Distributions/Continuous/Normal.cs +++ b/src/Managed/Distributions/Continuous/Normal.cs @@ -32,7 +32,8 @@ namespace MathNet.Numerics.Distributions using System.Collections.Generic; /// - /// Implements the univariate Normal (or Gaussian) distribution. + /// Implements the univariate Normal (or Gaussian) distribution. For details about this distribution, see + /// Wikipedia - Normal distribution. /// public class Normal : IContinuousDistribution { @@ -43,24 +44,28 @@ namespace MathNet.Numerics.Distributions /// /// Constructs a standard normal distribution. This is a normal distribution with mean 0.0 - /// and standard deviation 1.0. + /// and standard deviation 1.0. The distribution will + /// be initialized with the default random number generator. /// public Normal() : this(0.0, 1.0) { } /// - /// Construct a normal distribution with a particular mean and standard deviation. + /// Construct a normal distribution with a particular mean and standard deviation. The distribution will + /// be initialized with the default random number generator. /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. public Normal(double mean, double stddev) { SetParameters(mean, stddev); + RandomSource = new Random(); } /// - /// Constructs a normal distribution from a mean and standard deviation. + /// Constructs a normal distribution from a mean and standard deviation. The distribution will + /// be initialized with the default random number generator. /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. @@ -70,7 +75,8 @@ namespace MathNet.Numerics.Distributions } /// - /// Constructs a normal distribution from a mean and variance. + /// Constructs a normal distribution from a mean and variance. The distribution will + /// be initialized with the default random number generator. /// /// The mean of the normal distribution. /// The variance of the normal distribution. @@ -80,7 +86,8 @@ namespace MathNet.Numerics.Distributions } /// - /// Constructs a normal distribution from a mean and precision. + /// Constructs a normal distribution from a mean and precision. The distribution will + /// be initialized with the default random number generator. /// /// The mean of the normal distribution. /// The precision of the normal distribution. @@ -141,6 +148,9 @@ namespace MathNet.Numerics.Distributions } } + /// + /// The precision of the normal distribution. + /// public double Precision { get { return 1.0 / (mStdDev * mStdDev); } @@ -148,38 +158,117 @@ namespace MathNet.Numerics.Distributions } #region IDistribution implementation - public Random RandomNumberGenerator { get; set; } + /// + /// The random number generator which is used to draw random samples. + /// + public Random RandomSource { get; set; } + + /// + /// The mean of the normal distribution. + /// public double Mean { get { return mMean; } set { SetParameters(value, mStdDev); } } + + /// + /// The variance of the normal distribution. + /// public double Variance { get { return mStdDev * mStdDev; } set { SetParameters(mMean, value); } } + + /// + /// The standard deviation of the normal distribution. + /// public double StdDev { get { return mStdDev; } set { SetParameters(mMean, value); } } - public double Entropy { get { throw new NotImplementedException(); } } - public double Skewness { get { throw new NotImplementedException(); } } + + /// + /// The entropy of the normal distribution. + /// + public double Entropy { get { return Math.Log(mStdDev) + Constants.LogSqrt2PiE; } } + + /// + /// The skewness of the normal distribution. + /// + public double Skewness { get { return 0.0; } } #endregion #region IContinuousDistribution implementation - public double Mode { get { throw new NotImplementedException(); } } - public double Median { get { throw new NotImplementedException(); } } - public double Minimum { get { throw new NotImplementedException(); } } - public double Maximum { get { throw new NotImplementedException(); } } - public double Density(double x) { throw new NotImplementedException(); } - public double DensityLn(double x) { throw new NotImplementedException(); } + + /// + /// The mode of the normal distribution. + /// + public double Mode { get { return mMean; } } + + /// + /// The median of the normal distribution. + /// + public double Median { get { return mMean; } } + + /// + /// The minimum of the normal distribution. + /// + public double Minimum { get { return System.Double.NegativeInfinity; } } + + /// + /// The maximum of the normal distribution. + /// + public double Maximum { get { return System.Double.PositiveInfinity; } } + + /// + /// Computes the density of the normal distribution. + /// + /// The location at which to compute the density. + public double Density(double x) + { + double d = (x - mMean) / mStdDev; + return Math.Exp(-0.5*d*d) / (Constants.Sqrt2Pi*mStdDev); + } + + /// + /// Computes the log density of the normal distribution. + /// + /// The location at which to compute the log density. + public double DensityLn(double x) + { + double d = (x - mMean) / mStdDev; + return -0.5 * d * d - Math.Log(mStdDev) - Constants.LogSqrt2Pi; + } + public double CumulativeDistribution(double x) { throw new NotImplementedException(); } - public double Sample() { throw new NotImplementedException(); } - public IEnumerable Samples() { throw new NotImplementedException(); } + /// + /// Generates a sample from the normal distribution using the Box-Muller algorithm. + /// + public double Sample() + { + double r2; + return mMean + mStdDev * SampleBoxMuller(rng, r2); + } + + /// + /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. + /// + public IEnumerable Samples() + { + double r2; + + while (true) + { + double r1 = SampleBoxMuller(rng, r2); + yield return mean + stddev * r1; + yield return mean + stddev * r2; + } + } #endregion public double InverseCumulativeDistribution(double p) @@ -187,7 +276,55 @@ namespace MathNet.Numerics.Distributions throw new NotImplementedException(); } - public static double Sample(System.Random rng, double mean, double stddev) { throw new NotImplementedException(); } - public static IEnumerable Samples(System.Random rng, double mean, double stddev) { throw new NotImplementedException(); } + /// + /// Generates a sample from the normal distribution using the Box-Muller algorithm. + /// + /// The random number generator to use. + /// The mean of the normal distribution from which to generate samples. + /// The standard deviation of the normal distribution from which to generate samples. + public static double Sample(System.Random rng, double mean, double stddev) + { + double r2; + return mean + stddev * SampleBoxMuller(rng, r2); + } + + /// + /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. + /// + /// The random number generator to use. + /// The mean of the normal distribution from which to generate samples. + /// The standard deviation of the normal distribution from which to generate samples. + public static IEnumerable Samples(System.Random rng, double mean, double stddev) + { + double r2; + + while(true) + { + double r1 = SampleBoxMuller(rng, r2); + yield return mean + stddev * r1; + yield return mean + stddev * r2; + } + } + + /// + /// Samples a pair of standard normal distributed random variables using the Box-Muller algorithm. + /// + /// The random number generator to use. + /// The second random number. + internal static double SampleBoxMuller(System.Random rnd, out double r2) + { + double v1 = 2.0 * rnd.NextDouble() - 1.0; + double v2 = 2.0 * rnd.NextDouble() - 1.0; + double r = v1 * v1 + v2 * v2; + while (r >= 1.0 || r == 0.0) + { + v1 = 2.0 * rnd.NextDouble() - 1.0; + v2 = 2.0 * rnd.NextDouble() - 1.0; + r = v1 * v1 + v2 * v2; + } + double fac = System.Math.Sqrt(-2.0 * System.Math.Log(r) / r); + r2 = v2 * fac; + return v1 * fac; + } } } diff --git a/src/Managed/Distributions/IDistribution.cs b/src/Managed/Distributions/IDistribution.cs index 7ad434e6..86392e48 100644 --- a/src/Managed/Distributions/IDistribution.cs +++ b/src/Managed/Distributions/IDistribution.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.Distributions /// /// Gets or sets the random number generator which is used to generate random samples from the distribution. /// - Random RandomNumberGenerator { get; set; } + Random RandomSource { get; set; } /// /// The mean of the distribution.