Added preliminary KernelDensityEstimator

10 years ago · ff148b209d
2 changed files with 392 additions and 0 deletions
--- a/src/Numerics/Statistics/KernelDensityEstimator.cs
+++ b/src/Numerics/Statistics/KernelDensityEstimator.cs
@ -0,0 +1,200 @@
+// <copyright file="KernelDensityEstimator.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using MathNet.Numerics.Distributions;
+using MathNet.Numerics.Threading;
+
+namespace MathNet.Numerics.Statistics
+{
+  /// <summary>
+  /// An enum of the methods the <see cref="MathNet.Numerics.Statistics.KernelDensityEstimator"/> supports
+  /// for automatic bandwidth selection.
+  /// </summary>
+  public enum KDEBandwidthSelectionMethod
+  {
+    /// <summary>
+    /// TBD
+    /// </summary>
+    SilvermansRuleOfThumb,
+    /// <summary>
+    /// TBD
+    /// </summary>
+    SolveTheEquation
+  }
+
+  /// <summary>
+  /// The <see cref="MathNet.Numerics.Statistics.KernelDensityEstimator"/> supports several predefined Kernels.
+  /// Note that you can set your own custom kernel by setting <see cref="MathNet.Numerics.Statistics.KernelDensityEstimator.Kernel"/>
+  /// </summary>
+  public enum KDEKernelType
+  {
+    /// <summary>
+    /// A Gaussian kernel (PDF of Normal distribution with mean 0 and variance 1).
+    /// This kernel is the default.
+    /// </summary>
+    Gaussian,
+
+    /// <summary>
+    /// Epanechnikov Kernel
+    /// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0
+    /// </summary>
+    Epanechnikov,
+
+    /// <summary>
+    /// Uniform Kernel
+    /// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0
+    /// </summary>
+    Uniform,
+
+    /// <summary>
+    /// Triangular Kernel
+    /// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0
+    /// </summary>
+    Triangular,
+
+    /// <summary>
+    /// A custom kernel can be set by property <see cref="MathNet.Numerics.Statistics.KernelDensityEstimator.Kernel"/>
+    /// </summary>
+    Custom
+  }
+
+  /// <summary>
+  /// 
+  /// </summary>
+  public class KernelDensityEstimator
+  {
+    public KernelDensityEstimator(IList<double> samples)
+    {
+      _samples = samples;
+      KernelType = KDEKernelType.Gaussian;
+    }
+
+    public double EstimateDensity(double x)
+    {
+      var n = Samples.Count;
+      var estimate = CommonParallel.Aggregate(0, n,
+        i =>
+        {
+          var s = Samples[i];
+          return Kernel((x - s) / Bandwidth);
+        },
+        (a, b) => a + b,
+        0d) / (n * Bandwidth);
+
+      return estimate;
+    }
+
+    private readonly IList<double> _samples;
+    public IList<double> Samples
+    {
+      get { return _samples; }
+    }
+
+    private double _bandwidth = 1;
+    public double Bandwidth
+    {
+      get
+      {
+        return _bandwidth;
+      }
+      set
+      {
+        if (value <= 0)
+        {
+          throw new ArgumentException("The bandwidth must be a positive number!");
+        }
+        _bandwidth = value;
+      }
+    }
+
+    private KDEKernelType _kernelType;
+    public KDEKernelType KernelType
+    {
+      get { return _kernelType; }
+      set
+      {
+        switch (value)
+        {
+          case KDEKernelType.Gaussian:
+            {
+              Kernel = x => Normal.PDF(0.0, 1.0, x);
+              _kernelType = KDEKernelType.Gaussian;
+            }
+            break;
+          case KDEKernelType.Epanechnikov:
+            {
+              Kernel = x => Math.Abs(x) <= 1.0 ? 0.75 * (1 - x * x) : 0.0;
+              _kernelType = KDEKernelType.Epanechnikov;
+            }
+            break;
+          case KDEKernelType.Uniform:
+            {
+              Kernel = x => ContinuousUniform.PDF(-1.0, 1.0, x);
+              _kernelType = KDEKernelType.Uniform;
+            }
+            break;
+          case KDEKernelType.Triangular:
+            {
+              Kernel = x => Triangular.PDF(-1.0, 1.0, 0.0, x);
+              _kernelType = KDEKernelType.Triangular;
+            }
+            break;
+          case KDEKernelType.Custom:
+            throw new ArgumentException("In order to set a custom Kernel, property Kernel must be set directly.");
+        }
+      }
+    }
+
+    private Func<double, double> _kernel;
+    /// <summary>
+    /// Sets or Gets the Kernel used for the density estimate.
+    /// Setting the Kernel changes the <see cref="KernelDensityEstimator.KernelType"/> to <see cref="KDEKernelType.Custom"/>
+    /// A Kernel is a real function with Integral 1. Typically, it is also positive and symmetric about 0.
+    /// Note that none of these properties are checked.
+    /// </summary>
+    public Func<double, double> Kernel
+    {
+      get { return _kernel; }
+      set
+      {
+        _kernel = value;
+        _kernelType = KDEKernelType.Custom;
+      }
+    }
+
+    public double SelectBandwidth(KDEBandwidthSelectionMethod bandwidthSelectionMethod)
+    {
+      throw new NotImplementedException();
+    }
+  }
+}
--- a/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs
+++ b/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs
@ -0,0 +1,192 @@
+// <copyright file="KernelDensityEstimatorTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using NUnit.Framework;
+using MathNet.Numerics.Statistics;
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests
+{
+  /// <summary>
+  /// Kernel Density Estimator tests.
+  /// </summary>
+  [TestFixture, Category("Statistics")]
+  public class KernelDensityEstimatorTests
+  {
+    private readonly double[] _testData =
+    {
+      0.899822328897223,
+     -0.300111005615676,
+      1.029365712103099,
+     -0.345065971567321,
+      1.012801864262980,
+      0.629334584931419,
+     -0.213015082641055,
+     -0.865697308360524,
+     -1.043108301337627,
+     -0.270068812648099
+    };
+
+    [Test]
+    public void KDETestGaussianKernelBandwidth1()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+
+      Assert.AreEqual(KDEKernelType.Gaussian, kde.KernelType);
+      Assert.AreEqual(1.0d, kde.Bandwidth);
+      AssertHelpers.AlmostEqualRelative(0.398942280401433, kde.Kernel(0), 10); //Density of standard normal distribution at 0
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0.004115405028907, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.310485907659139, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.099698581377801, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestTriangularKernelBandwidth1()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.KernelType = KDEKernelType.Triangular;
+
+      Assert.AreEqual(KDEKernelType.Triangular, kde.KernelType);
+      Assert.AreEqual(1.0d, kde.Bandwidth);
+      Assert.AreEqual(1.0d, kde.Kernel(0)); //Density of standard normal distribution at 0
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.347688490533868, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.004216757636608, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestUniformKernelBandwidth1()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.KernelType = KDEKernelType.Uniform;
+
+      Assert.AreEqual(KDEKernelType.Uniform, kde.KernelType);
+      Assert.AreEqual(1.0d, kde.Bandwidth);
+      Assert.AreEqual(0.5d, kde.Kernel(0));
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.35, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.1, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestEpanechnikovKernelBandwidth1()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.KernelType = KDEKernelType.Epanechnikov;
+
+      Assert.AreEqual(KDEKernelType.Epanechnikov, kde.KernelType);
+      Assert.AreEqual(1.0d, kde.Bandwidth);
+      Assert.AreEqual(0.75d, kde.Kernel(0));
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.353803214812608, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.006248168996717, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestGaussianKernelBandwidth0p5()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.Bandwidth = 0.5d;
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(5.311490430807364e-007, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.369994803886827, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.032447347007482, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestGaussianKernelBandwidth2()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.Bandwidth = 2.0d;
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0.046875864115900, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.186580447512078, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.123339405007761, estimate, 10);
+    }
+
+    [Test]
+    public void KDETestCustomKernelBandwidth1()
+    {
+      var kde = new KernelDensityEstimator(_testData);
+      kde.Bandwidth = 1.0d;
+      kde.Kernel = x => 0.5d * Math.Exp(-Math.Abs(x)); //Picard-Kernel
+
+      Assert.AreEqual(KDEKernelType.Custom, kde.KernelType);
+      Assert.AreEqual(1.0d, kde.Bandwidth);
+      Assert.AreEqual(0.5d, kde.Kernel(0));
+
+      var estimate = kde.EstimateDensity(-3.5);
+      AssertHelpers.AlmostEqualRelative(0.018396636706009, estimate, 10);
+
+      estimate = kde.EstimateDensity(0);
+      AssertHelpers.AlmostEqualRelative(0.272675897096678, estimate, 10);
+
+      estimate = kde.EstimateDensity(2);
+      AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10);
+    }
+  }
+}