ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Statistics: Kernel Density Estimation: Cleanup

build
Christoph Ruegg 9 years ago
parent
ff148b209d
commit
da1630cd78
2 changed files with 282 additions and 282 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 274
      src/Numerics/Statistics/KernelDensityEstimator.cs
  2. 290
      src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs

274
src/Numerics/Statistics/KernelDensityEstimator.cs
View File

@ -3,7 +3,7 @@
// http://numerics.mathdotnet.com // http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics // 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 // Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation // obtaining a copy of this software and associated documentation
@ -29,172 +29,176 @@
using System; using System;
using System.Collections.Generic; using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.Distributions; using MathNet.Numerics.Distributions;
using MathNet.Numerics.Threading; using MathNet.Numerics.Threading;
namespace MathNet.Numerics.Statistics 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> /// <summary>
/// TBD /// An enum of the methods the <see cref="KernelDensityEstimator"/> supports
/// for automatic bandwidth selection.
/// </summary> /// </summary>
SilvermansRuleOfThumb, public enum KDEBandwidthSelectionMethod
/// <summary> {
/// TBD /// <summary>
/// </summary> /// TBD
SolveTheEquation /// </summary>
} SilvermansRuleOfThumb,
/// <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> /// <summary>
/// Epanechnikov Kernel /// TBD
/// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0 /// </summary>
/// </summary> SolveTheEquation
Epanechnikov, }
/// <summary> /// <summary>
/// Uniform Kernel /// The <see cref="KernelDensityEstimator"/> supports several predefined Kernels.
/// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0 /// Note that you can set your own custom kernel by setting <see cref="KernelDensityEstimator.Kernel"/>
/// </summary> /// </summary>
Uniform, 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> /// <summary>
/// Triangular Kernel /// Epanechnikov Kernel
/// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0 /// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0
/// </summary> /// </summary>
Triangular, 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="KernelDensityEstimator.Kernel"/>
/// </summary>
Custom
}
/// <summary> /// <summary>
/// A custom kernel can be set by property <see cref="MathNet.Numerics.Statistics.KernelDensityEstimator.Kernel"/> /// Kernel density estimation
/// </summary> /// </summary>
Custom public class KernelDensityEstimator
}
/// <summary>
///
/// </summary>
public class KernelDensityEstimator
{
public KernelDensityEstimator(IList<double> samples)
{ {
_samples = samples; public KernelDensityEstimator(IList<double> samples)
KernelType = KDEKernelType.Gaussian; {
} _samples = samples;
KernelType = KDEKernelType.Gaussian;
}
public double EstimateDensity(double x)
{
var n = Samples.Count;
var estimate = CommonParallel.Aggregate(0, n, public double EstimateDensity(double x)
i =>
{ {
var s = Samples[i]; var n = Samples.Count;
return Kernel((x - s) / Bandwidth); var estimate = CommonParallel.Aggregate(0, n,
}, i =>
(a, b) => a + b, {
0d) / (n * Bandwidth); var s = Samples[i];
return Kernel((x - s) / Bandwidth);
},
(a, b) => a + b,
0d) / (n * Bandwidth);
return estimate; return estimate;
} }
private readonly IList<double> _samples; private readonly IList<double> _samples;
public IList<double> Samples
{
get { return _samples; }
}
private double _bandwidth = 1; public IList<double> Samples
public double Bandwidth
{
get
{
return _bandwidth;
}
set
{
if (value <= 0)
{ {
throw new ArgumentException("The bandwidth must be a positive number!"); get { return _samples; }
} }
_bandwidth = value;
}
}
private KDEKernelType _kernelType; private double _bandwidth = 1;
public KDEKernelType KernelType
{ public double Bandwidth
get { return _kernelType; }
set
{
switch (value)
{ {
case KDEKernelType.Gaussian: get { return _bandwidth; }
{ set
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; if (value <= 0)
_kernelType = KDEKernelType.Epanechnikov; {
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); switch (value)
_kernelType = KDEKernelType.Uniform; {
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<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 = x => Triangular.PDF(-1.0, 1.0, 0.0, x); _kernel = value;
_kernelType = KDEKernelType.Triangular; _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<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) public double SelectBandwidth(KDEBandwidthSelectionMethod bandwidthSelectionMethod)
{ {
throw new NotImplementedException(); throw new NotImplementedException();
}
} }
}
} }

290
src/UnitTests/StatisticsTests/KernelDensityEstimatorTests.cs
View File

@ -3,7 +3,7 @@
// http://numerics.mathdotnet.com // http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics // 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 // Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation // obtaining a copy of this software and associated documentation
@ -28,165 +28,161 @@
// </copyright> // </copyright>
using System; using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using NUnit.Framework; using NUnit.Framework;
using MathNet.Numerics.Statistics; using MathNet.Numerics.Statistics;
namespace MathNet.Numerics.UnitTests.StatisticsTests namespace MathNet.Numerics.UnitTests.StatisticsTests
{ {
/// <summary> /// <summary>
/// Kernel Density Estimator tests. /// Kernel Density Estimator tests.
/// </summary> /// </summary>
[TestFixture, Category("Statistics")] [TestFixture, Category("Statistics")]
public class KernelDensityEstimatorTests public class KernelDensityEstimatorTests
{
private readonly double[] _testData =
{ {
0.899822328897223, private readonly double[] _testData =
-0.300111005615676, {
1.029365712103099, 0.899822328897223,
-0.345065971567321, -0.300111005615676,
1.012801864262980, 1.029365712103099,
0.629334584931419, -0.345065971567321,
-0.213015082641055, 1.012801864262980,
-0.865697308360524, 0.629334584931419,
-1.043108301337627, -0.213015082641055,
-0.270068812648099 -0.865697308360524,
}; -1.043108301337627,
-0.270068812648099
[Test] };
public void KDETestGaussianKernelBandwidth1()
{ [Test]
var kde = new KernelDensityEstimator(_testData); public void KDETestGaussianKernelBandwidth1()
{
Assert.AreEqual(KDEKernelType.Gaussian, kde.KernelType); var kde = new KernelDensityEstimator(_testData);
Assert.AreEqual(1.0d, kde.Bandwidth);
AssertHelpers.AlmostEqualRelative(0.398942280401433, kde.Kernel(0), 10); //Density of standard normal distribution at 0 Assert.AreEqual(KDEKernelType.Gaussian, kde.KernelType);
Assert.AreEqual(1.0d, kde.Bandwidth);
var estimate = kde.EstimateDensity(-3.5); AssertHelpers.AlmostEqualRelative(0.398942280401433, kde.Kernel(0), 10); //Density of standard normal distribution at 0
AssertHelpers.AlmostEqualRelative(0.004115405028907, estimate, 10);
var estimate = kde.EstimateDensity(-3.5);
estimate = kde.EstimateDensity(0); AssertHelpers.AlmostEqualRelative(0.004115405028907, estimate, 10);
AssertHelpers.AlmostEqualRelative(0.310485907659139, estimate, 10);
estimate = kde.EstimateDensity(0);
estimate = kde.EstimateDensity(2); AssertHelpers.AlmostEqualRelative(0.310485907659139, estimate, 10);
AssertHelpers.AlmostEqualRelative(0.099698581377801, estimate, 10);
} estimate = kde.EstimateDensity(2);
AssertHelpers.AlmostEqualRelative(0.099698581377801, estimate, 10);
[Test] }
public void KDETestTriangularKernelBandwidth1()
{ [Test]
var kde = new KernelDensityEstimator(_testData); public void KDETestTriangularKernelBandwidth1()
kde.KernelType = KDEKernelType.Triangular; {
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); var estimate = kde.EstimateDensity(-3.5);
Assert.AreEqual(1.0d, kde.Bandwidth); AssertHelpers.AlmostEqualRelative(0.046875864115900, estimate, 10);
Assert.AreEqual(1.0d, kde.Kernel(0)); //Density of standard normal distribution at 0
var estimate = kde.EstimateDensity(-3.5); estimate = kde.EstimateDensity(0);
AssertHelpers.AlmostEqualRelative(0, estimate, 10); AssertHelpers.AlmostEqualRelative(0.186580447512078, estimate, 10);
estimate = kde.EstimateDensity(0); estimate = kde.EstimateDensity(2);
AssertHelpers.AlmostEqualRelative(0.347688490533868, estimate, 10); AssertHelpers.AlmostEqualRelative(0.123339405007761, 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); [Test]
Assert.AreEqual(1.0d, kde.Bandwidth); public void KDETestCustomKernelBandwidth1()
Assert.AreEqual(0.75d, kde.Kernel(0)); {
var kde = new KernelDensityEstimator(_testData);
var estimate = kde.EstimateDensity(-3.5); kde.Bandwidth = 1.0d;
AssertHelpers.AlmostEqualRelative(0, estimate, 10); kde.Kernel = x => 0.5d * Math.Exp(-Math.Abs(x)); //Picard-Kernel
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(KDEKernelType.Custom, kde.KernelType);
Assert.AreEqual(1.0d, kde.Bandwidth); Assert.AreEqual(1.0d, kde.Bandwidth);
Assert.AreEqual(0.5d, kde.Kernel(0)); Assert.AreEqual(0.5d, kde.Kernel(0));
var estimate = kde.EstimateDensity(-3.5); var estimate = kde.EstimateDensity(-3.5);
AssertHelpers.AlmostEqualRelative(0.018396636706009, estimate, 10); AssertHelpers.AlmostEqualRelative(0.018396636706009, estimate, 10);
estimate = kde.EstimateDensity(0); estimate = kde.EstimateDensity(0);
AssertHelpers.AlmostEqualRelative(0.272675897096678, estimate, 10); AssertHelpers.AlmostEqualRelative(0.272675897096678, estimate, 10);
estimate = kde.EstimateDensity(2); estimate = kde.EstimateDensity(2);
AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10); AssertHelpers.AlmostEqualRelative(0.092580285110347, estimate, 10);
}
} }
}
} }

Write Preview
Loading…
Cancel
Save