From da1630cd7832fc8d71911a2ea228cba239c7009c Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 7 Jan 2018 16:04:41 +0100 Subject: [PATCH] Statistics: Kernel Density Estimation: Cleanup --- .../Statistics/KernelDensityEstimator.cs | 274 +++++++++-------- .../KernelDensityEstimatorTests.cs | 290 +++++++++--------- 2 files changed, 282 insertions(+), 282 deletions(-) diff --git a/src/Numerics/Statistics/KernelDensityEstimator.cs b/src/Numerics/Statistics/KernelDensityEstimator.cs index 1947b985..fd38641d 100644 --- a/src/Numerics/Statistics/KernelDensityEstimator.cs +++ b/src/Numerics/Statistics/KernelDensityEstimator.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// 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 @@ -29,172 +29,176 @@ using System; using System.Collections.Generic; -using System.Linq; -using System.Text; 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 + /// An enum of the methods the supports + /// for automatic bandwidth selection. /// - 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, + public enum KDEBandwidthSelectionMethod + { + /// + /// TBD + /// + SilvermansRuleOfThumb, - /// - /// Epanechnikov Kernel - /// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0 - /// - Epanechnikov, + /// + /// TBD + /// + SolveTheEquation + } /// - /// Uniform Kernel - /// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0 + /// The supports several predefined Kernels. + /// Note that you can set your own custom kernel by setting /// - Uniform, + public enum KDEKernelType + { + /// + /// A Gaussian kernel (PDF of Normal distribution with mean 0 and variance 1). + /// This kernel is the default. + /// + Gaussian, - /// - /// Triangular Kernel - /// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0 - /// - Triangular, + /// + /// 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 + } /// - /// A custom kernel can be set by property + /// Kernel density estimation /// - Custom - } - - /// - /// - /// - public class KernelDensityEstimator - { - public KernelDensityEstimator(IList samples) + public class KernelDensityEstimator { - _samples = samples; - KernelType = KDEKernelType.Gaussian; - } + 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 => + + + + public double EstimateDensity(double x) { - var s = Samples[i]; - return Kernel((x - s) / Bandwidth); - }, - (a, b) => a + b, - 0d) / (n * Bandwidth); + 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; - } + return estimate; + } - private readonly IList _samples; - public IList Samples - { - get { return _samples; } - } + private readonly IList _samples; - private double _bandwidth = 1; - public double Bandwidth - { - get - { - return _bandwidth; - } - set - { - if (value <= 0) + public IList Samples { - throw new ArgumentException("The bandwidth must be a positive number!"); + get { return _samples; } } - _bandwidth = value; - } - } - private KDEKernelType _kernelType; - public KDEKernelType KernelType - { - get { return _kernelType; } - set - { - switch (value) + private double _bandwidth = 1; + + public double Bandwidth { - case KDEKernelType.Gaussian: - { - Kernel = x => Normal.PDF(0.0, 1.0, x); - _kernelType = KDEKernelType.Gaussian; - } - break; - case KDEKernelType.Epanechnikov: + get { return _bandwidth; } + set { - Kernel = x => Math.Abs(x) <= 1.0 ? 0.75 * (1 - x * x) : 0.0; - _kernelType = KDEKernelType.Epanechnikov; + if (value <= 0) + { + throw new ArgumentException("The bandwidth must be a positive number!"); + } + + _bandwidth = value; } - break; - case KDEKernelType.Uniform: + } + + private KDEKernelType _kernelType; + + public KDEKernelType KernelType + { + get { return _kernelType; } + set { - Kernel = x => ContinuousUniform.PDF(-1.0, 1.0, x); - _kernelType = KDEKernelType.Uniform; + 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."); + } } - break; - case KDEKernelType.Triangular: + } + + 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 = x => Triangular.PDF(-1.0, 1.0, 0.0, x); - _kernelType = KDEKernelType.Triangular; + _kernel = value; + _kernelType = KDEKernelType.Custom; } - 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(); + public double SelectBandwidth(KDEBandwidthSelectionMethod bandwidthSelectionMethod) + { + throw new NotImplementedException(); + } } - } } diff --git a/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs b/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs index 3ed56a21..eee46c35 100644 --- a/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs +++ b/src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2016 Math.NET +// 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 @@ -28,165 +28,161 @@ // 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 { - /// - /// Kernel Density Estimator tests. - /// - [TestFixture, Category("Statistics")] - public class KernelDensityEstimatorTests - { - private readonly double[] _testData = + /// + /// Kernel Density Estimator tests. + /// + [TestFixture, Category("Statistics")] + public class KernelDensityEstimatorTests { - 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; + 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; - 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.046875864115900, estimate, 10); - var estimate = kde.EstimateDensity(-3.5); - AssertHelpers.AlmostEqualRelative(0, estimate, 10); + estimate = kde.EstimateDensity(0); + AssertHelpers.AlmostEqualRelative(0.186580447512078, 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; + estimate = kde.EstimateDensity(2); + AssertHelpers.AlmostEqualRelative(0.123339405007761, estimate, 10); + } - 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 + [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)); + 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); + 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(0); + AssertHelpers.AlmostEqualRelative(0.272675897096678, estimate, 10); - estimate = kde.EstimateDensity(2); - AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10); + estimate = kde.EstimateDensity(2); + AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10); + } } - } }