核密度估计Kernel Density Estimation(KDE)

前言

已知一个服从密度函数p(x)的随机变量x的N个观测x1,x2,…,xn，但不知p(x),这里p(x)称为期望学到的密度函数，试求p(x)。这个问题称为密度估计问题。为求样本集的密度函数，有两种方法：
1. 参数估计方法
简单来讲，即假定样本集符合某一概率分布，然后根据样本集拟合该分布中的参数，例如：似然估计，混合高斯，贝叶斯估计等，由于参数估计方法中需要加入主观的先验知识，往往很难拟合出与真实分布的模型；
2. 非参数估计方法
和参数估计不同，非参数估计并不加入任何先验知识，而是根据数据本身的特点、性质来拟合分布，这样能比参数估计方法得出更好的模型。核密度估计就是非参数估计中的一种，由Rosenblatt (1955)和Emanuel Parzen(1962)提出，又名Parzen窗（Parzen window）。Ruppert和Cline基于数据集密度函数聚类算法提出修订的核密度估计方法。

核密度估计

直方图是密度估计的一种非参数方法，虽然直方图法直观简单，但是由于样本数据有限，因此导致得到的p ̂(x)间断不连续，与生活常识不符。为了使p ̂(x)连续，每一个观测样本对密度的影响也应该是连续的，其对密度的影响力应该随着距离的增加而平滑减小。由此得到核密度估计公式：
p(x)=1N∑k=1N1hK(x−xkh)\text{p}(x)=\frac{1}{N}\sum\limits_{k=1}^{N}{\frac{1}{h}K(\frac{x-{{x}_{k}}}{h})} p(x)=N1k=1∑Nh1K(hx−xk)
其中，参数h称为带宽，K(x)称为核函数，K(x)满足条件：
K(x)≥0,∫K(x)dx=1K(x)\ge 0,\int{K(x)dx=1}K(x)≥0,∫K(x)dx=1
∫xK(x)dx=0\int{xK(x)dx=0}∫xK(x)dx=0
∫x2K(x)dx>0\int{{{x}^{2}}K(x)dx>0}∫x2K(x)dx>0
可见，核函数K(x)并不是唯一的，因此，核密度估计公式也不是唯一的。我们要做的是选取合适的带宽h和合适的核函数K(x)。

核函数

从支持向量机、meansift都接触过核函数，应该说核函数是一种理论概念，但每种核函数的功能都是不一样的，这里的核函数有Uniform,Triangular, Epanechnikov, Quartic, Triweight, Gaussian, Cosine等。
这些核函数的图像大致如下图：

此处对一些常用的核函数作以下说明：
Epanechnikov核：
K(x)=34(1−x2)I(∣x∣≤1)K(x)=\frac{3}{4}(1-{{x}^{2}})I(|x|\le 1)K(x)=43(1−x2)I(∣x∣≤1)
高斯核：
K(x)=12πexp⁡(−x22)K(x)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{{{x}^{2}}}{2})K(x)=2π1exp(−2x2)

如果你想了解更多关于核函数的内容，请点击这里

带宽h的选择

在核函数确定之后，比如上面选择的高斯核，那么高斯核的方差，也就是h（也叫带宽，也叫窗口，我们这里说的邻域）应该选择多大呢？不同的带宽会导致最后的拟合结果差别很大。同时上面也提到过，理论上h->0的，但h太小，邻域中参与拟合的点就会过少。那么借助机器学习的理论，我们当然可以使用交叉验证选择最好的h。另外，也有一个理论的推导给你选择h提供一些信息。
在样本集给定的情况下，我们只能对样本点的概率密度进行计算，那拟合过后的概率密度应该核计算的值更加接近才好，基于这一点，我们定义一个误差函数，然后最小化该误差函数便能为h的选择提供一个大致的方向。选择均平方积分误差函数(mean intergrated squared error)，该函数的定义是：
MISE(h)=E∫(^f(x)−f(x))2dxMISE(h)=E\int{{{\hat{(}f(x)-f(x))}^{2}}dx}MISE(h)=E∫(^f(x)−f(x))2dx
在weak assumptions下，MISE(h)=AMISE(h)+o(1nh+h4)MISE(h)=AMISE(h)+o(\frac{1}{nh}+{{h}^{4}})MISE(h)=AMISE(h)+o(nh1+h4)，其中AMISE为渐进的MISE，而AMISE为：
AMISE(h)=R(k)nh+14m2(K)2h4R(f′′)AMISE(h)=\frac{R(k)}{nh}+\frac{1}{4}{{m}_{2}}{{(K)}^{2}}{{h}^{4}}R({{f}^{''}})AMISE(h)=nhR(k)+41m2(K)2h4R(f′′)，其中：
R(k)=∫K(x)2dx,m2(k)=∫x2K(x)dxR(k)=\int{K{{(x)}^{2}}dx},{{m}_{2}}(k)=\int{{{x}^{2}}K(x)dx}R(k)=∫K(x)2dx,m2(k)=∫x2K(x)dx
最小化MISE(h)等价于最小化AMISE(h)，求导，令导数为0有：
∂∂hAMISE(h)=−R(k)nh2+m2(K)2h3R(f′′)=0\frac{\partial }{\partial h}AMISE(h)=-\frac{R(k)}{n{{h}^{2}}}+{{m}_{2}}{{(K)}^{2}}{{h}^{3}}R({{f}^{''}})=0∂h∂AMISE(h)=−nh2R(k)+m2(K)2h3R(f′′)=0
得：
hAMISE=R(K)15m2(K)25R(f′′)15n15{{h}_{AMISE}}=\frac{R{{(K)}^{\frac{1}{5}}}}{{{m}_{2}}{{(K)}^{\frac{2}{5}}}R{{({{f}^{''}})}^{\frac{1}{5}}}{{n}^{\frac{1}{5}}}}hAMISE=m2(K)52R(f′′)51n51R(K)51
当核函数确定之后，h公式里的R、m、f’'都可以确定下来，h便存在解析解。如果带宽不是固定的，其变化取决于估计的位置（balloon estimator）或样本点（逐点估计pointwise estimator），由此可以产生一个非常强大的方法称为自适应或可变带宽核密度估计。

参考

https://blog.csdn.net/unixtch/article/details/78556499
https://www.zhihu.com/question/27301358/answer/105267357?from=profile_answer_card
http://blog.csdn.net/yuanxing14/article/details/41948485