Distributions: maximum-likelihood estimation for normal, log-normal

src/Numerics/Distributions/Continuous/LogNormal.cs

@@ -24,11 +24,14 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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;
 
     /// <summary>
     /// Implements the univariate Log-Normal distribution. For details about this distribution, see 
@@ -83,6 +86,28 @@ namespace MathNet.Numerics.Distributions
             SetParameters(mu, sigma);
         }
 
+        /// <summary>
+        /// Constructs a log-normal distribution with the desired mean and variance. The distribution will
+        /// be initialized with the default <seealso cref="System.Random"/> random number generator.
+        /// </summary>
+        /// <param name="mean">The mean of the log-normal distribution.</param>
+        /// <param name="var">The variance of the log-normal distribution.</param>
+        /// <returns>a log-normal distribution.</returns>
+        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));
+        }
+
+        /// <summary>
+        /// Estimates the normal distribution parameters from sample data with maximum-likelihood.
+        /// </summary>
+        public static LogNormal Estimate(IEnumerable<double> samples)
+        {
+            var muSigma2 = samples.Select(s => Math.Log(s)).MeanVariance();
+            return new LogNormal(muSigma2.Item1, Math.Sqrt(muSigma2.Item2));
+        }
+
         /// <summary>
         /// A string representation of the distribution.
         /// </summary>

src/Numerics/Distributions/Continuous/Normal.cs

@@ -24,11 +24,13 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Statistics;
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Distributions
 {
     using System;
-    using System.Collections.Generic;
-    using Properties;
 
     /// <summary>
     /// 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));
         }
 
+        /// <summary>
+        /// Estimates the normal distribution parameters from sample data with maximum-likelihood.
+        /// </summary>
+        public static Normal Estimate(IEnumerable<double> samples)
+        {
+            var meanVariance = samples.MeanVariance();
+            return new Normal(meanVariance.Item1, Math.Sqrt(meanVariance.Item2));
+        }
+
         /// <summary>
         /// A string representation of the distribution.
         /// </summary>
@@ -201,8 +212,6 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
-        #region IDistribution implementation
-
         /// <summary>
         /// Gets or sets the random number generator which is used to draw random samples.
         /// </summary>
@@ -263,10 +272,6 @@ namespace MathNet.Numerics.Distributions
             get { return 0.0; }
         }
 
-        #endregion
-
-        #region IContinuousDistribution implementation
-
         /// <summary>
         /// Gets the mode of the normal distribution.
         /// </summary>
@@ -300,7 +305,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the density of the normal distribution.
+        /// Computes the density of the normal distribution (PDF), i.e. dP(X &lt;= x)/dx.
         /// </summary>
         /// <param name="mean">The mean of the normal distribution.</param>
         /// <param name="sdev">The standard deviation of the normal distribution.</param>
@@ -313,7 +318,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the log density of the normal distribution.
+        /// Computes the log density of the normal distribution (lnPDF), i.e. ln(dP(X &lt;= x)/dx).
         /// </summary>
         /// <param name="mean">The mean of the normal distribution.</param>
         /// <param name="sdev">The standard deviation of the normal distribution.</param>
@@ -326,7 +331,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the density of the normal distribution.
+        /// Computes the density of the normal distribution (PDF), i.e. dP(X &lt;= x)/dx.
         /// </summary>
         /// <param name="x">The location at which to compute the density.</param>
         /// <returns>the density at <paramref name="x"/>.</returns>
@@ -336,7 +341,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the log density of the normal distribution.
+        /// Computes the log density of the normal distribution (lnPDF), i.e. ln(dP(X &lt;= x)/dx).
         /// </summary>
         /// <param name="x">The location at which to compute the log density.</param>
         /// <returns>the log density at <paramref name="x"/>.</returns>
@@ -346,7 +351,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the cumulative distribution function of the normal distribution.
+        /// Computes the cumulative distribution function (CDF) of the normal distribution, i.e. P(X &lt;= x).
         /// </summary>
         /// <param name="mean">The mean of the normal distribution.</param>
         /// <param name="sdev">The standard deviation of the normal distribution.</param>
@@ -358,7 +363,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Computes the cumulative distribution function of the normal distribution.
+        /// Computes the cumulative distribution function (CDF) of the normal distribution, i.e. P(X &lt;= x).
         /// </summary>
         /// <param name="x">The location at which to compute the cumulative density.</param>
         /// <returns>the cumulative density at <paramref name="x"/>.</returns>
@@ -367,8 +372,6 @@ namespace MathNet.Numerics.Distributions
             return CumulativeDistribution(_mean, _stdDev, x);
         }
 
-        #endregion
-
         /// <summary>
         /// Computes the inverse cumulative distribution function of the normal distribution.
         /// </summary>

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);
         }
+
+        /// <summary>
+        /// Can estimate distribution parameters.
+        /// </summary>
+        [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);
+        }
     }
 }

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);
         }
+
+        /// <summary>
+        /// Can estimate distribution parameters.
+        /// </summary>
+        [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);
+        }
     }
 }