1 核密度估计
对给定的一维各自独立的从小到大排列的样本\(x_1,x_2,…,x_n\),核密度估计函数定义如下
(1-1)
$$
f(x)=\dfrac{1}{nh}\sum_{i=1}^nK(\dfrac{x-x_i}{h})
$$
公式里\(h>0\)为窗宽,\(K(x)\)为核函数。
为了方便表示，设第\(i\)个概率密度贡献函数为
(1-2)
$$
f_i(x)=K(\dfrac{x-x_i}{h})
$$
根据累积分布函数定义，设第\(i\)个累积分布贡献函数为
(1-3)
$$
c_i(x)=\int_{-\infty}^x f_i(t)dt
$$
因此最终累积分布函数可表示为
(1-4)
$$c(x)=\int_{-\infty}^x f(t)dt=\dfrac{1}{nh}\sum_{i=1}^nc_i(x)$$
2 最佳窗宽
在实际数据处理中，选择合适的核函数以及窗宽，才能更好地还原真实的数据分布。
2.1 最小化目标函数
为了获取最佳窗宽,本质上需要假定某种误差模型(最小化目标函数),然后找到一个最佳\(h\)使得这个误差模型的值最小，选择核函数以及求窗宽本质是一个纯优化问题。
对于这类优化问题，比较常见的是建立如下的积分误差模型
(2-1)
$$
\begin{cases}
\text{概率密度积分平方误差}:&e(h)=\int (f(x)-\bar{f}(x))^2dx
\\
\text{概率密度积分均方误差}:&e(h)=E[\int (f(x)-\bar{f}(x))^2dx]
\end{cases}
$$
公式里\(\bar{f}(x)\)为概率密度真值模型，很遗憾地是，我们本身就不知道这个模型长什么样，很难获得比较真的值，因此很多算法都在对\(\bar{f}(x)\)提出各种近似，进而最小化误差函数。
这里，我们需要理解，概率密度函数本质反应的是某个点累积概率变化程度，因此它本质反应的是累积概率的一阶导（从定义上就可以看出），这时可以先直接得到累积概率密度值，然后使用某种算法获得累积概率各个点的一阶导值，这样就可以得到近似的\(\bar{f}(x)\)，不过往往得到的值噪声较大。
因此，这里不使用公式(2-1)里的任意一个误差模型。我们换一种思路构建最小化目标函数。由于公式(1-4)已经获得了假设模型的累积分布函数\(c(x)\)，因此可以构建如下最小二乘误差
(2-2)
$$
e(h)=\sum_{i=1}^n(c(x_i)-\bar{c}(x_i))^2
$$
公式里\(\bar{c}(x)\)为真值累积分布值，当样本足够多时，可用如下方式获取
(2-3)
$$
\bar{c}(x_i)\approx \dfrac{i}{n}
$$
注意, 公式(2-3)要求样本从小到大排序且无重复,有重复需要稍微处理一下。
这里需要注意的是，之所以选择累积分布作为误差值，因为概率密度函数"真值"\(\bar{f}(x)\)往往是对累积分布值再加工后获取，既然是再加工，就可能存在真值损失，这时直接使用累积分布值作为直接度量，可以一定程度减少真值损失带来的误差！
因此，使用(2-2)最小化目标的重点就变成了对核函数进行积分。积分得到累计分布函数后，就可以使用(2-2)套用对应求解器优化出最佳窗宽。
注意，如果直接使用样本数据点作为每个核函数里的\(x_i\)，直接使用(1-4)的累计概率函数去进行逼近时，可以发现，对于无重复的样本数据，最优解会出现在窗宽接近0的位置，这种不是预期内的，这时需要单独再特殊处理。因此实际可以选择部分点作为核函数点，再进行最优带宽求解。
2.2 核函数
这里列举常见的核函数。由于我们有采样数据，因此可以获得采样数据的均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)，如果我们采样数据足够多，理论上核密度估计模型的期望\(\mu\)与方差\(\sigma^2\)与采样数据的均值、方差越接近越好，因此在核函数列举过程中，会罗列出模型对应的期望、方差。
2.2.1 三角核函数
对于参数 \(a>0\), 这里定义如下三角形核函数
(2-4)
$$
K(x)=\begin{cases}
0&,x<-a
\\
\\
\dfrac{1}{a}(1-\dfrac{1}{a}|x|)&,-a\leq x \leq a
\\
\\
0&,x>a
\end{cases}
$$
根据(1-3)定义，可以得到累积分布贡献函数
(2-5)
$$
c_i(x)=\begin{cases}
0&,x<x_i-ah
\\
\\
\dfrac{1}{2a^2h}(x_i-ah-x)^2&,x_i-ah\leq x\leq x_i
\\
\\
\dfrac{1}{2}h+\dfrac{1}{2a^2h}(x-x_i)(2ah+x_i-x)&,x_i<x\leq x_i+ah
\\
\\
h&,x>x_i+ah
\end{cases}
$$
根据期望与方差定义，可以很容易得到模型期望\(\mu\)、方差\(\sigma^2\)与样本数据均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)的关系
(2-6)
$$
\begin{cases}
\mu&=\bar{\mu}
\\
\sigma^2&=\bar{\sigma}^2+\dfrac{a^2}{6}h^2
\end{cases}
$$
特别地，在Matlab里的三角函数核取\(a=\sqrt{6}\), 在维基百科上的三角核取\(a=1\)。
2.2.2 抛物线核函数
对于参数\(a>0\)，这里定义如下抛物线核函数
(2-7)
$$
K(x)=\begin{cases}
0&,x<-a
\\
\\
\dfrac{3}{4a}(1-\dfrac{x^2}{a^2})&,-a\leq x\leq a
\\
\\
0&,x>a
\end{cases}
$$
根据(1-3)定义，可以得到累积分布贡献函数
(2-8)
$$
c_i(x)=\begin{cases}
0&,x<x_i-ah
\\
\\
\dfrac{1}{4a^3h^2}(ah+x-x_i)^2(2ah-x+x_i)&,x_i-ah\leq x \leq x_i+ah
\\
\\
h&,x>x_i+ah
\end{cases}
$$
根据期望与方差定义，可以很容易得到模型期望\(\mu\)、方差\(\sigma^2\)与样本数据均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)的关系
(2-9)
$$
\begin{cases}
\mu&=\bar{\mu}
\\
\\
\sigma^2&=\bar{\sigma}^2+\dfrac{a^2}{5}h^2
\end{cases}
$$
特别地，当\(a=1\)时，就得到了常见的Epanechnikov核。
2.2.3 高斯核函数
对于参数\(a>0\),这里定义如下高斯核函数
(2-10)
$$
K(x)=\dfrac{1}{\sqrt{2\pi}a}\exp(-\dfrac{1}{2a^2}(x-b)^2)
$$
根据(1-3)定义，可以得到累积分布贡献函数
(2-11)
$$
c_i(x)=\dfrac{h}{2}(1-\mathbf{erf}(\dfrac{bh+x_i-x}{\sqrt{2}ah}))
$$
根据期望定义，可以得到模型期望\(\mu\)与样本均值直接的关系
(2-12)
$$
\mu=\bar{\mu}-hb
$$
很多时候，我们希望模型的期望逼近样本的均值，因此这时核函数设置\(b=0\)，此时模型期望\(\mu\)、方差\(\sigma^2\)与样本数据均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)有如下关系
(2-13)
$$
\begin{cases}
\mu&=\bar{\mu}
\\
\sigma^2&=\bar{\sigma}^2+a^2h^2
\end{cases}
$$
特别地，当\(a=1,b=0\)时，就得到了常见的高斯核。
2.2.4 余弦核函数
对于参数\(a>0\),这里定义如下余弦核函数
(2-14)
$$
K(x)=\begin{cases}
0&,x<-a
\\
\\
\dfrac{\pi}{4a}\cos(\dfrac{\pi}{2a}x)&,-a\leq x \leq a
\\
\\
0&,x>a
\end{cases}
$$
根据(1-3)定义，可以得到累积分布贡献函数
(2-15)
$$
c_i(x)=\begin{cases}
0&,x<x_i-ah
\\
\\
\dfrac{h}{2}(1+\sin(\dfrac{\pi(x-x_i)}{2ah}))&,x_i-ah\leq x \leq x_i+ah
\\
\\
h&,x>x_i+ah
\end{cases}
$$
根据期望与方差定义，可以得到模型期望\(\mu\)、方差\(\sigma^2\)与样本数据均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)的关系
(2-16)
$$
\begin{cases}
\mu&=\bar{\mu}
\\
\\
\sigma^2&=\bar{\sigma}^2+(1-\dfrac{8}{\pi^2})a^2h^2
\end{cases}
$$
特别地，当\(a=1\)时，就得到常见的余弦核函数。
2.2.5 逻辑核函数
这里定义核函数
(2-17)
$$
K(x)=\dfrac{1}{2+\exp(x)+\exp(-x)}
$$
根据(1-3)定义，可以得到累积分布贡献函数
(2-18)
$$
c_i(x)=h-\dfrac{h}{1+\exp(\dfrac{x-x_i}{h})}
$$
根据期望与方差定义，可以得到模型期望\(\mu\)、方差\(\sigma^2\)与样本数据均值\(\bar{\mu}\)、方差\(\bar{\sigma}^2\)的关系
(2-19)
$$
\begin{cases}
\mu&=\bar{\mu}
\\
\\
\sigma^2&=\bar{\sigma}^2+\dfrac{\pi^2}{3}h^2
\end{cases}
$$
3 经验窗宽
生活中，大部分模型逼近正态分布，因此如果不想麻烦地去优化窗宽的话，可以认为数据符合正态分布，这时根据Silverman提出的拇指准则(一种经验准则)，可使用如下方式计算窗宽
(3-1)
$$
\begin{cases}
\mu&=\dfrac{1}{n}\sum_{i=1}^nx_i
\\
\\
\sigma&=\sqrt{\dfrac{\sum_{i=1}^n(x_i-\mu)^2}{n-1}}
\\
\\
h&=(\dfrac{4}{3n})^{1/5}\sigma
\end{cases}
$$
4 后话
在核密度估计里，为了加速计算，往往会对求解域进行截断，具体可以参看三角核函数、抛物线核函数，在实际数据处理中需要充分利用这个优势。同样地，高斯核函数本质也可以截断，可以根据实际情况进行处理。