// <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-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.
|
// </copyright>
|
|
using System;
|
using System.Collections.Generic;
|
using IStation.Numerics.Distributions;
|
using IStation.Numerics.Threading;
|
|
namespace IStation.Numerics.Statistics
|
{
|
/// <summary>
|
/// Kernel density estimation (KDE).
|
/// </summary>
|
public static class KernelDensity
|
{
|
/// <summary>
|
/// Estimate the probability density function of a random variable.
|
/// </summary>
|
/// <remarks>
|
/// The routine assumes that the provided kernel is well defined, i.e. a real non-negative function that integrates to 1.
|
/// </remarks>
|
public static double Estimate(double x, double bandwidth, IList<double> samples, Func<double, double> 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;
|
}
|
|
/// <summary>
|
/// Estimate the probability density function of a random variable with a Gaussian kernel.
|
/// </summary>
|
public static double EstimateGaussian(double x, double bandwidth, IList<double> samples)
|
{
|
return Estimate(x, bandwidth, samples, GaussianKernel);
|
}
|
|
/// <summary>
|
/// Estimate the probability density function of a random variable with an Epanechnikov kernel.
|
/// The Epanechnikov kernel is optimal in a mean square error sense.
|
/// </summary>
|
public static double EstimateEpanechnikov(double x, double bandwidth, IList<double> samples)
|
{
|
return Estimate(x, bandwidth, samples, EpanechnikovKernel);
|
}
|
|
/// <summary>
|
/// Estimate the probability density function of a random variable with a uniform kernel.
|
/// </summary>
|
public static double EstimateUniform(double x, double bandwidth, IList<double> samples)
|
{
|
return Estimate(x, bandwidth, samples, UniformKernel);
|
}
|
|
/// <summary>
|
/// Estimate the probability density function of a random variable with a triangular kernel.
|
/// </summary>
|
public static double EstimateTriangular(double x, double bandwidth, IList<double> samples)
|
{
|
return Estimate(x, bandwidth, samples, TriangularKernel);
|
}
|
|
/// <summary>
|
/// A Gaussian kernel (PDF of Normal distribution with mean 0 and variance 1).
|
/// This kernel is the default.
|
/// </summary>
|
public static double GaussianKernel(double x)
|
{
|
return Normal.PDF(0.0, 1.0, x);
|
}
|
|
/// <summary>
|
/// Epanechnikov Kernel:
|
/// x => Math.Abs(x) <= 1.0 ? 3.0/4.0(1.0-x^2) : 0.0
|
/// </summary>
|
public static double EpanechnikovKernel(double x)
|
{
|
return Math.Abs(x) <= 1.0 ? 0.75 * (1 - x * x) : 0.0;
|
}
|
|
/// <summary>
|
/// Uniform Kernel:
|
/// x => Math.Abs(x) <= 1.0 ? 1.0/2.0 : 0.0
|
/// </summary>
|
public static double UniformKernel(double x)
|
{
|
return ContinuousUniform.PDF(-1.0, 1.0, x);
|
}
|
|
/// <summary>
|
/// Triangular Kernel:
|
/// x => Math.Abs(x) <= 1.0 ? (1.0-Math.Abs(x)) : 0.0
|
/// </summary>
|
public static double TriangularKernel(double x)
|
{
|
return Triangular.PDF(-1.0, 1.0, 0.0, x);
|
}
|
}
|
}
|
