概率密度函数的核估计
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
from scipy import stats
from typing import *
核密度估计（kernel density estimation）
核密度估计法是一种通过某个（连续的）概率分布的样本来估计这个概率分布的密度函数的方法。
说到用样本来估计概率密度，最基础的就应该是“直方图”了。我们可以把直方图看作是一个几乎处处连续的函数，用这样一个连续的函数作为未知概率分布的近似。对样本点
，取分点
，直方图这样一个连续函数：
当样本数量趋于无穷并且划分区间长度趋于0时，是几乎处处收敛与原概率分布的密度函数的。
以下代码生成了100个标准正态分布随机数并画出了它们的直方图
sample = np.random.randn(100)
sns.distplot(sample, kde=False, bins=15, norm_hist=True)
plt.show()
但是用直方图来近似有一个问题，就是它不够光滑，同时，如果分布在两侧或一侧有重尾，这时直方图的部分小区间可能只有很少点甚至没有点，部分小区间则集中了过多的点，等距概率直方图就不能很好地反映分布密度形状。
核密度估计是一种比较平滑地估计未知分布概率密度的方法。
我们可以针对每一个
，用
来估计
（其中
表示集合的元素个数）
即：
如果把上面的区间改为左开右闭区间
, 就有：
，
是经验分布函数。
即
是对经验分布函数用差分近似估计
导数的结果。
这种估计叫做「Rosenblatt 直方图估计」
设函数
Rosenblatt 直方图估计可以写成
这里的
叫做核函数。
def kernel_density(K, sample, h):
"""
K: density function, h: bandwidth
返回样本的核密度估计函数
"""
sample = np.array(sample)
f = lambda y: np.mean(np.vectorize(K)((y - sample) / h)) / h
return np.vectorize(f)
y = np.arange(-2, 2, 0.1)
func = kernel_density(lambda x: int(abs(x) <= 1/2), sample, h=1)
plt.figure(figsize=(15, 4.5))
for i, h in enumerate([0.2, 0.5, 1, 3]):
plt.subplot(1, 4, i + 1)
func = kernel_density(lambda x: int(abs(x) <= 1/2), sample, h=h)
plt.plot(y, func(y))
plt.show()
上图是用Rosenblatt直方图方法估计的标准正态分布样本点的概率密度。
注意到：
- 当
很小时，密度估计很不光滑，在每个
处有一个尖锐的峰而没有观测值的地方密度估计值非常小，估计偏差小而方差大。
- 当
较大时，估计比较光滑，估计偏差大而方差小。
这个
我们一般叫做带宽（bandwidth），它的选取需要平衡偏差和方差。渐近地取
, 核密度估计的均方误差为
。
除了Rosenblatt直方图估计，还有一些其它的核函数：
比如说高斯核函数，用它来估计就具有非常好的光滑性。sns.displot函数的kde=True就会使用高斯核密度估计来拟合样本！
Gauss 核密度估计
K = lambda x: 1 * np.exp(-x**2/2) / np.sqrt(2 * np.pi)
plt.figure(figsize=(15, 4.5))
for i, h in enumerate([0.2, 0.6, 1, 3]):
plt.subplot(1, 4, i + 1)
func = kernel_density(K, sample, h=h)
plt.plot(y, func(y))
plt.show()
下图是标准的概率分布，可以看到，选取比较合适的bandwidth，高斯核密度估计能够很好地近似原分布！
plt.plot(y, stats.norm.pdf(y)); plt.show()
二次曲线核
K = lambda x: 3 / 4 * (1 - x**2) * (abs(x) <= 1)
func = kernel_density(K, sample, h=1)
plt.figure(figsize=(15, 4.5))
for i, h in enumerate([0.2, 0.6, 1, 3]):
plt.subplot(1, 4, i + 1)
func = kernel_density(K, sample, h=h)
plt.plot(y, func(y))
plt.show()
关于厚尾分布
sample = np.random.exponential(size=100)
sns.distplot(sample, norm_hist=True, kde=False)
<matplotlib.axes._subplots.AxesSubplot at 0x7f90b57333c8>
关于指数分布这种厚尾分布，直方图显得很无能为力，但是核密度估计法的效果是非常稳定的！
「以下是真实分布与核密度方法近似的比较：」
y = np.arange(-1, 7, 0.05)
plt.plot(y, stats.expon.pdf(y))
[<matplotlib.lines.Line2D at 0x7f90cb621128>]
K = lambda x: 3 / 4 * (1 - x**2) * (abs(x) <= 1)
func = kernel_density(K, sample, h=1)
plt.plot(y, func(y))
[<matplotlib.lines.Line2D at 0x7f90b5808ef0>]
可以看到核密度估计能够把分布的“尾巴”给近似出来！
参考：【1】韦来生.数理统计