概率密度函数的核估计

2020-11-05 11:39:39 浏览数 (2)
代码语言：javascript复制
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）

核密度估计法是一种通过某个（连续的）概率分布的样本来估计这个概率分布的密度函数的方法。

说到用样本来估计概率密度，最基础的就应该是“直方图”了。我们可以把直方图看作是一个几乎处处连续的函数，用这样一个连续的函数作为未知概率分布的近似。对样本点

y_{1}, y_{2}, ldots, y_{n}

，取分点

t_0 < t_1 < cdots < t_m

，直方图这样一个连续函数：

当样本数量趋于无穷并且划分区间长度趋于0时，是几乎处处收敛与原概率分布的密度函数的。

以下代码生成了100个标准正态分布随机数并画出了它们的直方图

代码语言：javascript复制
sample = np.random.randn(100)
代码语言：javascript复制
sns.distplot(sample, kde=False, bins=15, norm_hist=True)
plt.show()

但是用直方图来近似有一个问题，就是它不够光滑，同时，如果分布在两侧或一侧有重尾，这时直方图的部分小区间可能只有很少点甚至没有点，部分小区间则集中了过多的点，等距概率直方图就不能很好地反映分布密度形状。

核密度估计是一种比较平滑地估计未知分布概率密度的方法。

我们可以针对每一个

y

，用

tilde{p}(y)=frac{#left(left{y_{i}: y_{i} inleft[y-frac{h}{2}, y frac{h}{2}right]right}right) / n}{h}

来估计

p(y)

（其中

#

表示集合的元素个数）

即：

begin{aligned} tilde{p}(y) &=frac{1}{h} frac{1}{n} sum_{i=1}^{n} I_{left[y-frac{h}{2}, y frac{h}{2}right]}left(y_{i}right) \ &=frac{1}{n} sum_{i=1}^{n} frac{1}{h} I_{left[-frac{1}{2}, frac{1}{2}right]}left(frac{y-y_{i}}{h}right), y in(-infty, infty) end{aligned}

如果把上面的区间改为左开右闭区间

(y-frac{h}{2}, y frac{h}{2}]

, 就有：

tilde{p}(y)=frac{F_{n}left(y frac{h}{2}right)-F_{n}left(y-frac{h}{2}right)}{h}, y in(-infty, infty)

，

F_n(x)

是经验分布函数。

即

tilde{p}(y)

是对经验分布函数用差分近似估计

F(x)

导数的结果。

这种估计叫做「Rosenblatt 直方图估计」

设函数

K_{1}(x)=I_{left[-frac{1}{2}, frac{1}{2}right]}(x),

Rosenblatt 直方图估计可以写成

tilde{p}(y)=frac{1}{n} sum_{i=1}^{n} frac{1}{h} K_{1}left(frac{y-y_{i}}{h}right)

这里的

K_1(x)

叫做核函数。

代码语言：javascript复制
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)
代码语言：javascript复制
y = np.arange(-2, 2, 0.1)
代码语言：javascript复制
func = kernel_density(lambda x: int(abs(x) <= 1/2), sample, h=1)
代码语言：javascript复制
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直方图方法估计的标准正态分布样本点的概率密度。

注意到：

  • 当
h

很小时，密度估计很不光滑，在每个

y_i

处有一个尖锐的峰而没有观测值的地方密度估计值非常小，估计偏差小而方差大。

  • 当
h

较大时，估计比较光滑，估计偏差大而方差小。

这个

h

我们一般叫做带宽（bandwidth），它的选取需要平衡偏差和方差。渐近地取

h=O(n^{-1/5})

, 核密度估计的均方误差为

O(n^{−4/5})

。

除了Rosenblatt直方图估计，还有一些其它的核函数：

begin{aligned} K_{2}(x) &=(1-|x|) I_{[-1,1]}(x)(text {三角形核}) \ K_{3}(x) &=frac{3}{4}left(1-x^{2}right) I_{[-1,1]}(x)(text {二次曲线核}) \ K_{4}(x) &=frac{15}{16}left(1-x^{2}right)^{2} I_{[-1,1]}(x)(text {双二次核}) \ K_{5}(x) &=frac{70}{81}left(1-|x|^{3}right)^{3} I_{[-1,1]}(x)(text {双三次核}) \ K_{6}(x) &=frac{1}{sqrt{2 pi}} e^{-frac{x^{2}}{2}}(text { 高斯核 }) end{aligned}

比如说高斯核函数，用它来估计就具有非常好的光滑性。sns.displot函数的kde=True就会使用高斯核密度估计来拟合样本！

Gauss 核密度估计
代码语言：javascript复制
K = lambda x: 1 * np.exp(-x**2/2) / np.sqrt(2 * np.pi)
代码语言：javascript复制
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，高斯核密度估计能够很好地近似原分布！

代码语言：javascript复制
plt.plot(y, stats.norm.pdf(y)); plt.show()
二次曲线核
代码语言：javascript复制
K = lambda x: 3 / 4 * (1 - x**2) * (abs(x) <= 1)
代码语言：javascript复制
func = kernel_density(K, sample, h=1)
代码语言：javascript复制
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()
关于厚尾分布
代码语言：javascript复制
sample = np.random.exponential(size=100)
代码语言：javascript复制
sns.distplot(sample, norm_hist=True, kde=False)
代码语言：javascript复制
<matplotlib.axes._subplots.AxesSubplot at 0x7f90b57333c8>

关于指数分布这种厚尾分布，直方图显得很无能为力，但是核密度估计法的效果是非常稳定的！

「以下是真实分布与核密度方法近似的比较：」

代码语言：javascript复制
y = np.arange(-1, 7, 0.05)
代码语言：javascript复制
plt.plot(y, stats.expon.pdf(y))
代码语言：javascript复制
[<matplotlib.lines.Line2D at 0x7f90cb621128>]
代码语言：javascript复制
K = lambda x: 3 / 4 * (1 - x**2) * (abs(x) <= 1)
func = kernel_density(K, sample, h=1)
代码语言：javascript复制
plt.plot(y, func(y))
代码语言：javascript复制
[<matplotlib.lines.Line2D at 0x7f90b5808ef0>]

可以看到核密度估计能够把分布的“尾巴”给近似出来！

参考：【1】韦来生.数理统计

kernel bandwidth

0 人点赞