概率密度函数的核估计

用户3577892

发布于 2020-11-05 11:39:39

1.9K0

发布于 2020-11-05 11:39:39

文章被收录于专栏：数据科学CLUB

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个标准正态分布随机数并画出了它们的直方图

sample = np.random.randn(100)

sns.distplot(sample, kde=False, bins=15, norm_hist=True)
plt.show()

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

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

我们可以针对每一个

，用

\tilde{p}(y)=\frac{\#\left(\left\{y_{i}: y_{i} \in\left[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)

叫做核函数。

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直方图方法估计的标准正态分布样本点的概率密度。

注意到：

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

y_i

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

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

这个

我们一般叫做带宽（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 核密度估计

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】韦来生.数理统计

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-10-30，如有侵权请联系 cloudcommunity@tencent.com 删除

kernel

bandwidth

本文分享自数据科学CLUB 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

kernel

bandwidth

登录后参与评论

0 条评论

热度