核密度估计KDE原理及Python实现（PDF｜CDF｜Kernel｜KDE）

PDF和CDF：CDF vs. PDF: What is the Difference?
KDE：Kernel Density Estimator explained step by step
本文只是从理解的角度介绍KDE，没有涉及复杂的公式推导

1.PDF和CDF

概率密度函数（Probability Density Function，PDF）和累积分布函数（Cumulative Distribution Function，CDF）都是在处理连续随机变量统计问题中非常重要的概念，两个函数都提供了对概率的分析，但是有不同的目的和数据分布视角。

1.1 概率密度函数PDF

PDF 代表概率密度函数，它是统计学中的一个重要概念，用于理解与连续随机变量相关的概率。它是一条平滑的曲线，显示了在一定范围内不同结果出现的可能性。

比如，对于某个城市某天的问题，PDF可以显示出温度落在某些特定范围的可能性，例如在20度-25度之间的可能性。

PDF不会给出特定值的概率，而是给出变量落在特定值附近一小段区间的概率，区间段的PDF曲线面积表示变量落在该范围内的概率。要找到单个值的概率，需要计算该点的 PDF 积分，也就是找到该特定值处的曲线下面积。

1.2 累计分布函数CDF

CDF代表累积分布函数，是对PDF的补充，提供了与随机变量相关的概率累积视图，CDF代表了一个随机变量小于等于某个特定值的概率。

CDF从负值对应的0开始，组建随着随机变量取值的增加变成1，对于离散随机变量，CDF 逐步增加，与每个可能结果的概率相对应。对于连续随机变量，CDF 平稳增加，并反映不同间隔内的组合概率

CDF	PDF
一个随机变量小于等于一个特定值的概率	一个随机变量落在一个特定值周围小区间的可能性
累积概率的总量	概率密度的瞬时变化率
提供特定值或者特定范围的概率	不提供特定值概率，但是提供每个值周围的概率密度
非递减函数	非负的曲线，曲线下面积总和为1

• PDF：随机变量在某个特定值附近的分布情况，$f(x)$
• CDF：随机变量在某个值之前的总体分布情况，$F(x)$

CDF是PDF的积分
$$F(x)={\int }_{-\infty }^{x}f(t)dt$$

1.3 正态分布举例

标准正态分布Standard Normal Distribution是均值为0，标准差为1的正态分布，也称为高斯分布Gaussian distribution，下面是正态分布的PDF，高斯分布的PDF为$\phi (x)$，CDF为$\varphi (z)$


下图为高斯分布的PDF和CDF表示

2.核密度估计KDE

当给定一组数据时，我们想要知道关于数据分布的信息，也就是概率密度函数PDF，一些数据样本比较容易拟合成常见的密度函数，比如正态分布，泊松分布等等，这时可以采用最大似然方法（the maximum likelihood）来拟合密度函数。
但是一些情况下，数据分布是没有规律的，不能直接对应上常见的PDFs上。这种情况下，Kernel Density Estimator (KDE) 提供了一种对数据分布比较合理和乐观的表达。

2.1 Kernel函数

Kernel：用来决定每个数据点周围密度分布的概率密度函数

Kernel Density Estimator：所有点Kernel的组合

比如乐高积木和汽车模型的关系，汽车模型是由很多块乐高积木组成的，Kernel就是乐高积木，KDE就是汽车模型，只是此处的乐高积木只有一种类型，也就是Kernel函数

首先定义一个数据点的Kernel函数，此处使用正态分布（即高斯分布表示），均值为0，方差为1
$$K(x)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}}$$

那么对于第$i$个数据点，其Kernel函数为

$$K(x-x_i)$$

为了让曲线变窄或者变宽，可以添加一个常数$h$（kernel bandwidth）

$$K(\frac{x-x_i}{h})$$
由于正态分布下的曲线面积要保持为1，h越大，曲线越宽，面积越大，因此函数需要再除以h
$$\frac{1}{h}K(\frac{x-x_i}{h})$$

如下图表示了不同h值下的高斯分布曲线图，h越大，越宽也越矮，保持面积和为1

2.2 Kernel bandwidth

在Kernel函数中有一个比较关键的参数，$h$（kernel bandwidth），通常有两种设置带宽的方法：Scott原则和Silverman原则

其中$n$是样本数据点数，${\sigma }_d$是从样本数据中得出的标准差。一些第三方数据分析库，比如Scipy，Seaborn，SKlearn等提供的KDE模块都能设置带框选择方法。本文在2.3中实现的代码为自设带宽。

2.3 Kernel Density Estimator

在确定了第一个数据点的Kernel函数后，同样的对于第二个数据点的kernel函数也有
$$\frac{1}{h}K(\frac{x-x_2}{h})$$

把两个数据点的kernel加起来，但是需要注意的是，每个数据点的kernel函数区域面积都是1，合起来的函数区域面积也要为1，因此需要取平均值

$$f(x)=\frac{1}{2}(\frac{1}{h}K(\frac{x-x_1}{h})+\frac{1}{h}K(\frac{x-x_2}{h}))$$

如下图是两个数据点的kernel函数组合成KDE，可以理解为x1的kernel曲线和x2的kernel曲线取均值变成KDE

同样的思路，n个数据点的KDE函数如下，K是核函数，h是带宽，n是样本数量

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

kernel函数的选择有很多种，不同的选择对概率分布函数的估计效果也不同，此处选择的高斯kernel

下面用一个简单的数据样本显示不同$h$值下的KDE函数

import numpy as np
import matplotlib.pyplot as plt

"""
（1）定义数据：原始数据+依照定义间隔生成均匀分布的x
"""
dataset = np.array([1.33, 0.3, 0.97, 1.1, 0.1, 1.4, 0.4])  # 数据点
x_range = np.linspace(np.min(dataset) - 0.3, np.max(dataset) + 0.3, num=600)  # kernel函数的x轴
sample_nb = len(dataset)

"""
（2）设置参数：不同的bandwidth
"""
bandwidth = [0.3, 0.1, 0.03]
color_list = ['goldenrod', 'black', 'maroon']
alpha_list = [0.8, 1, 0.8]
width_list = [1.7, 2.5, 1.7]

"""
（3）定义函数：kernel函数，KDE绘图函数
"""
# kernel函数
def kernel_func(x):
    return np.exp(-x**2/2)/np.sqrt(2*np.pi)

# 绘制每个bandwidth下的KDE
def plot(h, c, a, w):
    total_sum = 0
    # 遍历每个数据点，计算kernel函数下的y值总和
    for index, x in enumerate(dataset):
        total_sum += kernel_func((x_range - x) / h)
        plt.annotate(r'$x_{}$'.format(index+1), xy=[x, 0.13], horizontalalignment='center', fontsize=18)
    # 取平均得到组合后的KDE的y值
    y_range = total_sum / (h * sample_nb)
    # 绘制KDE曲线
    plt.plot(x_range, y_range, color=c, alpha=a, linewidth=w, label=f'{h}')
    # 标注每个数据点
    plt.plot(dataset, np.zeros_like(dataset), 's', markersize=8, color='black')

"""
（4）结果展示
"""
for dot_h, dot_c, dot_a, dot_w in zip(bandwidth, color_list, alpha_list, width_list):
    plot(dot_h, dot_c, dot_a, dot_w)

plt.xlabel('$x$', fontsize=22)
plt.ylabel('$f(x)$', fontsize=22)
plt.legend(fontsize=14, shadow=True, title='$h$', title_fontsize=16)
plt.show()

2.4 KDE with Python Seaborn

seaborn中kdeplot使用了更加复杂和优化的KDE计算过程，比如加入了内部优化技术，权重平滑和边界处理等，并且Seaborn的bw_adjust参数调整的是Seaborn内部计算的带宽因子，不代表带宽，因此下面的图和上面自己实现的并不是完全对应

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

dataset = np.array([1.33, 0.3, 0.97, 1.1, 0.1, 1.4, 0.4])  # 数据点
sns.set()
fig, ax = plt.subplots(figsize=(10, 4))
sns.kdeplot(ax=ax, data=dataset, bw_adjust=0.3)
ax.plot(dataset, np.zeros_like(dataset)+0.05, 's', markersize=8, color='black')
for index, x in enumerate(dataset):
    ax.annotate(r'$x_{}$'.format(index + 1), xy=[x, 0.13], horizontalalignment='center', fontsize=18)
plt.show()