diff --git a/src/Numerics/Statistics/KernelDensity.cs b/src/Numerics/Statistics/KernelDensity.cs new file mode 100644 index 00000000..06c91bb4 --- /dev/null +++ b/src/Numerics/Statistics/KernelDensity.cs @@ -0,0 +1,133 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2018 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. +// + +using System; +using System.Collections.Generic; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.Threading; + +namespace MathNet.Numerics.Statistics +{ + /// + /// Kernel density estimation (KDE). + /// + public static class KernelDensity + { + /// + /// Estimate the probability density function of a random variable. + /// + /// + /// The routine assumes that the provided kernel is well defined, i.e. a real non-negative function that integrates to 1. + /// + public static double Estimate(double x, double bandwidth, IList samples, Func kernel) + { + if (bandwidth <= 0) + { + throw new ArgumentException("The bandwidth must be a positive number!"); + } + + var n = samples.Count; + var estimate = CommonParallel.Aggregate(0, n, + i => kernel((x - samples[i]) / bandwidth), + (a, b) => a + b, + 0d) / (n * bandwidth); + + return estimate; + } + + /// + /// Estimate the probability density function of a random variable with a Gaussian kernel. + /// + public static double EstimateGaussian(double x, double bandwidth, IList samples) + { + return Estimate(x, bandwidth, samples, GaussianKernel); + } + + /// + /// Estimate the probability density function of a random variable with an Epanechnikov kernel. + /// The Epanechnikov kernel is optimal in a mean square error sense. + /// + public static double EstimateEpanechnikov(double x, double bandwidth, IList samples) + { + return Estimate(x, bandwidth, samples, EpanechnikovKernel); + } + + /// + /// Estimate the probability density function of a random variable with a uniform kernel. + /// + public static double EstimateUniform(double x, double bandwidth, IList samples) + { + return Estimate(x, bandwidth, samples, UniformKernel); + } + + /// + /// Estimate the probability density function of a random variable with a triangular kernel. + /// + public static double EstimateTriangular(double x, double bandwidth, IList samples) + { + return Estimate(x, bandwidth, samples, TriangularKernel); + } + + /// + /// A Gaussian kernel (PDF of Normal distribution with mean 0 and variance 1). + /// This kernel is the default. + /// + public static double GaussianKernel(double x) + { + return Normal.PDF(0.0, 1.0, x); + } + + /// + /// Epanechnikov Kernel: + /// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0 + /// + public static double EpanechnikovKernel(double x) + { + return Math.Abs(x) <= 1.0 ? 0.75 * (1 - x * x) : 0.0; + } + + /// + /// Uniform Kernel: + /// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0 + /// + public static double UniformKernel(double x) + { + return ContinuousUniform.PDF(-1.0, 1.0, x); + } + + /// + /// Triangular Kernel: + /// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0 + /// + public static double TriangularKernel(double x) + { + return Triangular.PDF(-1.0, 1.0, 0.0, x); + } + } +} diff --git a/src/Numerics/Statistics/KernelDensityEstimator.cs b/src/Numerics/Statistics/KernelDensityEstimator.cs deleted file mode 100644 index fd38641d..00000000 --- a/src/Numerics/Statistics/KernelDensityEstimator.cs +++ /dev/null @@ -1,204 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// -// Copyright (c) 2009-2018 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. -// - -using System; -using System.Collections.Generic; -using MathNet.Numerics.Distributions; -using MathNet.Numerics.Threading; - -namespace MathNet.Numerics.Statistics -{ - /// - /// An enum of the methods the supports - /// for automatic bandwidth selection. - /// - public enum KDEBandwidthSelectionMethod - { - /// - /// TBD - /// - SilvermansRuleOfThumb, - - /// - /// TBD - /// - SolveTheEquation - } - - /// - /// The supports several predefined Kernels. - /// Note that you can set your own custom kernel by setting - /// - public enum KDEKernelType - { - /// - /// A Gaussian kernel (PDF of Normal distribution with mean 0 and variance 1). - /// This kernel is the default. - /// - Gaussian, - - /// - /// Epanechnikov Kernel - /// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0 - /// - Epanechnikov, - - /// - /// Uniform Kernel - /// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0 - /// - Uniform, - - /// - /// Triangular Kernel - /// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0 - /// - Triangular, - - /// - /// A custom kernel can be set by property - /// - Custom - } - - /// - /// Kernel density estimation - /// - public class KernelDensityEstimator - { - public KernelDensityEstimator(IList 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 _samples; - - public IList 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 _kernel; - - /// - /// Sets or Gets the Kernel used for the density estimate. - /// Setting the Kernel changes the to - /// 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. - /// - public Func Kernel - { - get { return _kernel; } - set - { - _kernel = value; - _kernelType = KDEKernelType.Custom; - } - } - - public double SelectBandwidth(KDEBandwidthSelectionMethod bandwidthSelectionMethod) - { - throw new NotImplementedException(); - } - } -} diff --git a/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs b/src/UnitTests/StatisticsTests/KernelDensityTests.cs similarity index 60% rename from src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs rename to src/UnitTests/StatisticsTests/KernelDensityTests.cs index eee46c35..6689d7b9 100644 --- a/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs +++ b/src/UnitTests/StatisticsTests/KernelDensityTests.cs @@ -37,7 +37,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests /// Kernel Density Estimator tests. /// [TestFixture, Category("Statistics")] - public class KernelDensityEstimatorTests + public class KernelDensityTests { private readonly double[] _testData = { @@ -56,132 +56,103 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests [Test] public void KDETestGaussianKernelBandwidth1() { - var kde = new KernelDensityEstimator(_testData); + //Density of standard normal distribution at 0 + AssertHelpers.AlmostEqualRelative(0.398942280401433, KernelDensity.GaussianKernel(0), 10); - 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); + var estimate = KernelDensity.EstimateGaussian(-3.5d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.004115405028907, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateGaussian(0.0d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.310485907659139, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateGaussian(2.0d, 1.0d, _testData); 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 + Assert.AreEqual(1.0d, KernelDensity.TriangularKernel(0)); - var estimate = kde.EstimateDensity(-3.5); + var estimate = KernelDensity.EstimateTriangular(-3.5d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateTriangular(0.0d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.347688490533868, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateTriangular(2.0d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.004216757636608, estimate, 10); } [Test] public void KDETestUniformKernelBandwidth1() { - var kde = new KernelDensityEstimator(_testData); - kde.KernelType = KDEKernelType.Uniform; + Assert.AreEqual(0.5d, KernelDensity.UniformKernel(0)); - 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); + var estimate = KernelDensity.EstimateUniform(-3.5d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateUniform(0.0d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.35, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateUniform(2.0d, 1.0d, _testData); 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)); + Assert.AreEqual(0.75d, KernelDensity.EpanechnikovKernel(0)); - var estimate = kde.EstimateDensity(-3.5); + var estimate = KernelDensity.EstimateEpanechnikov(-3.5d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateEpanechnikov(0.0d, 1.0d, _testData); AssertHelpers.AlmostEqualRelative(0.353803214812608, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateEpanechnikov(2.0d, 1.0d, _testData); 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); + var estimate = KernelDensity.EstimateGaussian(-3.5d, 0.5d, _testData); AssertHelpers.AlmostEqualRelative(5.311490430807364e-007, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateGaussian(0.0d, 0.5d, _testData); AssertHelpers.AlmostEqualRelative(0.369994803886827, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateGaussian(2.0d, 0.5d, _testData); 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); + var estimate = KernelDensity.EstimateGaussian(-3.5d, 2.0d, _testData); AssertHelpers.AlmostEqualRelative(0.046875864115900, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.EstimateGaussian(0.0d, 2.0d, _testData); AssertHelpers.AlmostEqualRelative(0.186580447512078, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.EstimateGaussian(2.0d, 2.0d, _testData); 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)); + double Kernel(double x) => 0.5d * Math.Exp(-Math.Abs(x)); + Assert.AreEqual(0.5d, Kernel(0)); - var estimate = kde.EstimateDensity(-3.5); + var estimate = KernelDensity.Estimate(-3.5d, 1.0d, _testData, Kernel); AssertHelpers.AlmostEqualRelative(0.018396636706009, estimate, 10); - estimate = kde.EstimateDensity(0); + estimate = KernelDensity.Estimate(0.0d, 1.0d, _testData, Kernel); AssertHelpers.AlmostEqualRelative(0.272675897096678, estimate, 10); - estimate = kde.EstimateDensity(2); + estimate = KernelDensity.Estimate(2.0d, 1.0d, _testData, Kernel); AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10); } }