非参数估计—核密度估计（Kernel Density Estimation）

本文所有概念仅适用于非参数估计。

1. 核函数定义

定义一个新的核函数用于非参数估计，必须满足下面三个条件：

• 归一性（Normalization）：
  $${\int }_{-\infty }^{+\infty }K(u)du=1$$

• 偶函数对称性（Even-function Symmetry）：
  $$K(-u)=K(u)$$

• 非负性（Non-negative）:
  $$K(u)\ge 0$$

2. 核密度估计的直观理解

为什么要进行核密度估计？

拿到一堆数据，如果我们要看他们分布，最先想到的是采用频数直方图的形式，统计每个区间的样本分布情况。

但是直方图容易受到组数和组距的影响，产生较大的差异，如下图所示。

总结，直方图用于密度估计存在一些缺点：1）不光滑；2）容易受带宽和数据范围的影响。

上述绘图代码如下：

import matplotlib.pyplot as plt
import numpy as np

def generate_gauss_data(mu, std, size=(1024,1)):
    '''
    生成均值mu和标准差std的正态分布数据
    :param mu:
    :param std:
    :param size: 生成数据的大小
    :return:
    '''
    data = np.random.normal(loc=mu, scale=std, size=size)
    return data

data = generate_gauss_data(mu=0.5, std=1, size=(1024, 1))

#plt.hist(data, bins=30, density=False, alpha=0.5, label='Histogram (bins=30)')
plt.hist(data, bins=15, density=False, alpha=0.5, label='Histogram (bins=15)')
plt.grid(linestyle='--', linewidth=1, alpha=0.6)
plt.xlabel('x')
plt.ylabel('Frequency')
plt.legend()
plt.show()

什么是核密度估计？

通俗解释就是估计每个观测点的核密度，然后对所有的观测点密度加权平均。

In nonparametric statistics, a kernel is a weighting function used in non-parametric estimation techniques.

绘图代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 生成一些随机数据点
data = np.random.normal(0, 1, size=5)

# 定义带宽
h = 0.3

# 定义核函数（这里用高斯核）
def gaussian_kernel(x, xi, bandwidth):
    return norm.pdf((x - xi) / bandwidth) / bandwidth

# 绘图范围
x = np.linspace(-4, 4, 1000)

# 单个数据点的核密度分布示意图
plt.figure(figsize=(12, 8))
for xi in data:
    plt.plot(x, gaussian_kernel(x, xi, h), alpha=0.4)

# 计算核密度估计（加权平均）
kde = np.mean([gaussian_kernel(x, xi, h) for xi in data], axis=0)

# 绘制核密度估计曲线
plt.plot(x, kde, color='black', linewidth=1, label='Kernel Density Estimation')

# 绘制数据点位置
plt.scatter(data, np.zeros_like(data), marker='x', color='red', label='Data Points')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

3. 核密度计算公式

定义带宽核函数（Scaled Kernel）

带宽核函数$K_h$是由标准的核函数$K$缩放得到的，定义如下：
$$K_h(x-X_i)=\frac{1}{h}K(\frac{x-X_i}{h})$$
式中，$x$是当前我们希望估计密度的位置，$X_i$是第$i$个样本点，$h$是带宽（控制核函数宽度和平滑程度）。

为什么要定义带宽核函数？

因为原始的核函数无法调节带宽，所以需要引入该参数进行平滑（smoothing parameter）。

证明：${\int }_{-\infty }^{+\infty }{K}_h(x-X_i)du=1$。

$$\begin{array}{rl}{\int }_{-\infty }^{+\infty }{K}_h(x-X_i)dx& =\frac{1}{h}{\int }_{-\infty }^{+\infty }K(\frac{x-X_i}{h})dx\\ & \stackrel{\underset{\mathrm{u}=(\mathrm{x}-{\mathrm{X}}_{\mathrm{i}})/h}{}}{=}\frac{1}{h}{\int }_{-\infty }^{+\infty }K(u)d(uh+X_i)\\ & =\frac{1}{h}\times h\times {\int }_{-\infty }^{+\infty }K(u)du\\ & ={\int }_{-\infty }^{+\infty }K(u)du=1\end{array}$$

注：上面引入了$u=(x-X_i)/h$进行变量代换，则有$dx=d(uh+X_i)=hd(u)$。

核函数的一般形式，$K^{\ast }(u)=\lambda K(\lambda u)$，其中$\lambda$是一个大于零的常数。

核密度估计公式

这里所述的核函数相当于是一个权重函数（weighting function），即以数据点$X_i$为中心的权重函数：

• 对于距离观测样本点$X_i$近的区域赋予大的权重
• 对于距离观测样本$X_i$远的区域赋予小的权重

即$K(|x-X_i|)$，又因为核函数是个偶函数（对称函数），则可以写为$K(x-X_i)$。

引入权重缩放之后，即为
$$\frac{1}{h}K(\frac{x-X_i}{h})$$
然后，计算每个观测点附近的核密度，在求均值

$$\frac{1}{n}\sum _{i=1}^n\frac{1}{h}K(\frac{x-X_i}{h})$$

最终的核密度估计公式如下：
$$\boxed{\hat{f}(x)=\frac{1}{n}\sum _{i=1}^n{K}_h(x-X_i)=\frac{1}{nh}\sum _{i=1}^{n}K(\frac{x-X_i}{h})}$$
采用高斯核函数：
$$K(u)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{u^2}{2}}$$
将高斯核函数带入上面的核密度估计公式，则有
$$\boxed{\hat{f}(x)=\frac{1}{nh}\sum _{i=1}^n\frac{1}{\sqrt{2\pi }}{e}^{-\frac{(x-X_i{)}^2}{2h^2}}}$$
式中，$x$是当前我们希望估计密度的置，$X_i$是第$i$个样本点，$h$是带宽（控制核函数宽度和平滑程度），$n$是样本数量。对于高斯函数而言[2]，一般根据经验带宽$h=1.06\times \sigma \times n^{-1/5}$，其中$\sigma$是数据的标准差。

代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, gaussian_kde
np.random.seed(42)


# 观测数据X
data = np.random.normal(0, 1, size=100)

# 计算带宽（Silverman经验法则）
n = len(data)
sigma = np.std(data, ddof=1)
h = 1.06 * sigma * n ** (-1 / 5)

def gaussian_kernel(u):
    return (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * u ** 2)

def my_kde(x, data, bandwidth):
    n = len(data)
    kde_values = np.zeros_like(x)
    for xi in data:
        kde_values += gaussian_kernel((x - xi) / bandwidth)
    return kde_values / (n * bandwidth)


x_plot = np.linspace(-4, 4, 1000) # x轴范围
kde_manual = my_kde(x_plot, data, h)

# 内置函数计算核密度估计
kde_builtin = gaussian_kde(data, bw_method=h / sigma)
kde_builtin_values = kde_builtin(x_plot)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x_plot, kde_manual, label='my_kde', linewidth=2, color='red')
plt.plot(x_plot, kde_builtin_values, label='Built-in KDE', linewidth=2, linestyle='--', color='blue')
plt.hist(data, bins=20, density=True, alpha=0.3, color='grey', label='Histogram')
plt.xlabel('x')
plt.ylabel('Density')
plt.title('Kernel Density Estimation (KDE)')
plt.legend()
plt.grid(True)
plt.show()

写在最后

周末花了半天+一晚上形成此文，瞟了下右下角的时间，不知不觉已经是第二天的凌晨。由于时间仓促，上述只是一个浅显的理解，还有很多理解不彻底或者理解错误的地方，欢迎批评指正！

本文暂且定义为一个初稿，后续继续修改，并透彻理解。

参考资料：

1. Kernel (statistics) | Wikipedia
2. Kernel density estimation | Wikipedia
3. 核密度估计（kernel density estimation） | 三点水
4. 什么是核密度估计？如何感性认识？ - 慧航的回答 - 知乎
5. 什么是核密度估计？如何感性认识？ - 浮生六记的回答 - 知乎