-
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】韦来生.数理统计
转载：https://blog.csdn.net/weixin_41503009/article/details/109396849