From 31583dcdd5949f1d5c2a7e5decd644c1f7fa7018 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 15 Aug 2013 09:12:56 +0200 Subject: [PATCH] Distributions: maximum-likelihood estimation for normal, log-normal --- .../Distributions/Continuous/LogNormal.cs | 29 +++++++++++++-- .../Distributions/Continuous/Normal.cs | 35 ++++++++++--------- .../Continuous/LogNormalTests.cs | 17 +++++++++ .../Continuous/NormalTests.cs | 17 +++++++++ 4 files changed, 80 insertions(+), 18 deletions(-) diff --git a/src/Numerics/Distributions/Continuous/LogNormal.cs b/src/Numerics/Distributions/Continuous/LogNormal.cs index d5a88a74..87332318 100644 --- a/src/Numerics/Distributions/Continuous/LogNormal.cs +++ b/src/Numerics/Distributions/Continuous/LogNormal.cs @@ -24,11 +24,14 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Linq; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Statistics; +using System.Collections.Generic; + namespace MathNet.Numerics.Distributions { using System; - using System.Collections.Generic; - using Properties; /// /// Implements the univariate Log-Normal distribution. For details about this distribution, see @@ -83,6 +86,28 @@ namespace MathNet.Numerics.Distributions SetParameters(mu, sigma); } + /// + /// Constructs a log-normal distribution with the desired mean and variance. The distribution will + /// be initialized with the default random number generator. + /// + /// The mean of the log-normal distribution. + /// The variance of the log-normal distribution. + /// a log-normal distribution. + public static LogNormal WithMeanVariance(double mean, double var) + { + var sigma2 = Math.Log(var / (mean * mean) + 1.0); + return new LogNormal(Math.Log(mean) - sigma2 / 2.0, Math.Sqrt(sigma2)); + } + + /// + /// Estimates the normal distribution parameters from sample data with maximum-likelihood. + /// + public static LogNormal Estimate(IEnumerable samples) + { + var muSigma2 = samples.Select(s => Math.Log(s)).MeanVariance(); + return new LogNormal(muSigma2.Item1, Math.Sqrt(muSigma2.Item2)); + } + /// /// A string representation of the distribution. /// diff --git a/src/Numerics/Distributions/Continuous/Normal.cs b/src/Numerics/Distributions/Continuous/Normal.cs index ba145633..b0d88988 100644 --- a/src/Numerics/Distributions/Continuous/Normal.cs +++ b/src/Numerics/Distributions/Continuous/Normal.cs @@ -24,11 +24,13 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.Properties; +using MathNet.Numerics.Statistics; +using System.Collections.Generic; + namespace MathNet.Numerics.Distributions { using System; - using System.Collections.Generic; - using Properties; /// /// Implements the univariate Normal (or Gaussian) distribution. For details about this distribution, see @@ -138,6 +140,15 @@ namespace MathNet.Numerics.Distributions return new Normal(mean, 1.0/Math.Sqrt(precision)); } + /// + /// Estimates the normal distribution parameters from sample data with maximum-likelihood. + /// + public static Normal Estimate(IEnumerable samples) + { + var meanVariance = samples.MeanVariance(); + return new Normal(meanVariance.Item1, Math.Sqrt(meanVariance.Item2)); + } + /// /// A string representation of the distribution. /// @@ -201,8 +212,6 @@ namespace MathNet.Numerics.Distributions } } - #region IDistribution implementation - /// /// Gets or sets the random number generator which is used to draw random samples. /// @@ -263,10 +272,6 @@ namespace MathNet.Numerics.Distributions get { return 0.0; } } - #endregion - - #region IContinuousDistribution implementation - /// /// Gets the mode of the normal distribution. /// @@ -300,7 +305,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the density of the normal distribution. + /// Computes the density of the normal distribution (PDF), i.e. dP(X <= x)/dx. /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. @@ -313,7 +318,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the log density of the normal distribution. + /// Computes the log density of the normal distribution (lnPDF), i.e. ln(dP(X <= x)/dx). /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. @@ -326,7 +331,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the density of the normal distribution. + /// Computes the density of the normal distribution (PDF), i.e. dP(X <= x)/dx. /// /// The location at which to compute the density. /// the density at . @@ -336,7 +341,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the log density of the normal distribution. + /// Computes the log density of the normal distribution (lnPDF), i.e. ln(dP(X <= x)/dx). /// /// The location at which to compute the log density. /// the log density at . @@ -346,7 +351,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the cumulative distribution function of the normal distribution. + /// Computes the cumulative distribution function (CDF) of the normal distribution, i.e. P(X <= x). /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. @@ -358,7 +363,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Computes the cumulative distribution function of the normal distribution. + /// Computes the cumulative distribution function (CDF) of the normal distribution, i.e. P(X <= x). /// /// The location at which to compute the cumulative density. /// the cumulative density at . @@ -367,8 +372,6 @@ namespace MathNet.Numerics.Distributions return CumulativeDistribution(_mean, _stdDev, x); } - #endregion - /// /// Computes the inverse cumulative distribution function of the normal distribution. /// diff --git a/src/UnitTests/DistributionTests/Continuous/LogNormalTests.cs b/src/UnitTests/DistributionTests/Continuous/LogNormalTests.cs index 76e6480c..6a0d12c6 100644 --- a/src/UnitTests/DistributionTests/Continuous/LogNormalTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/LogNormalTests.cs @@ -517,5 +517,22 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous var n = new LogNormal(mu, sigma); AssertHelpers.AlmostEqual(f, n.CumulativeDistribution(x), 8); } + + /// + /// Can estimate distribution parameters. + /// + [TestCase(0.0, 0.0)] + [TestCase(10.0, 0.1)] + [TestCase(-5.0, 1.0)] + [TestCase(0.0, 5.0)] + [TestCase(10.0, 50.0)] + public void CanEstimateParameters(double mu, double sigma) + { + var original = new LogNormal(mu, sigma, new Random(100)); + var estimated = LogNormal.Estimate(original.Samples().Take(10000)); + + AssertHelpers.AlmostEqual(mu, estimated.Mu, 2); + AssertHelpers.AlmostEqual(sigma, estimated.Sigma, 2); + } } } diff --git a/src/UnitTests/DistributionTests/Continuous/NormalTests.cs b/src/UnitTests/DistributionTests/Continuous/NormalTests.cs index d1a6a8ad..d1bfa5c9 100644 --- a/src/UnitTests/DistributionTests/Continuous/NormalTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/NormalTests.cs @@ -491,5 +491,22 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous var n = Normal.WithMeanStdDev(5.0, 2.0); AssertHelpers.AlmostEqual(x, n.InverseCumulativeDistribution(f), 15); } + + /// + /// Can estimate distribution parameters. + /// + [TestCase(0.0, 0.0)] + [TestCase(10.0, 0.1)] + [TestCase(-5.0, 1.0)] + [TestCase(0.0, 5.0)] + [TestCase(10.0, 50.0)] + public void CanEstimateParameters(double mean, double stddev) + { + var original = new Normal(mean, stddev, new Random(100)); + var estimated = Normal.Estimate(original.Samples().Take(10000)); + + AssertHelpers.AlmostEqual(mean, estimated.Mean, 2); + AssertHelpers.AlmostEqual(stddev, estimated.StdDev, 2); + } } }