CQF预备知识：Python相关库 -- 核密度估计 scipy.stats

文中内容仅限技术学习与代码实践参考，市场存在不确定性，技术分析需谨慎验证，不构成任何投资建议。

核密度估计

在统计学中，一个常见的任务是从一组数据样本中估计随机变量的概率密度函数（PDF），这一任务被称为密度估计。完成这一任务最著名的工具是直方图。直方图是一种用于可视化的有用工具（主要是因为每个人都能理解它），但它并没有高效地利用可用数据。核密度估计（KDE）是用于相同任务的更高效工具。scipy.stats.gaussian_kde 估计器可用于估计单变量和多变量数据的概率密度函数。如果数据是单峰的，它效果最佳。

单变量估计

我们从少量数据开始，以便了解 scipy.stats.gaussian_kde 的工作原理以及带宽选择的不同选项。从概率密度函数中采样的数据以蓝色短划线显示在图的底部（这被称为地毯图）：

>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt

>>> x1 = np.array([-7, -5, 1, 4, 5], dtype=np.float64)
>>> kde1 = stats.gaussian_kde(x1)
>>> kde2 = stats.gaussian_kde(x1, bw_method='silverman')

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)

>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # 地毯图
>>> x_eval = np.linspace(-10, 10, num=200)
>>> ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott 法则")
>>> ax.plot(x_eval, kde2(x_eval), 'r-', label="Silverman 法则")

我们看到 Scott 法则和 Silverman 法则之间几乎没有区别，并且在数据量有限的情况下，带宽选择可能有点太宽了。我们可以定义自己的带宽函数，以获得不太平滑的结果。

>>> def my_kde_bandwidth(obj, fac=1./5):
...     """我们使用 Scott 法则，乘以一个常数因子。"""
...     return np.power(obj.n, -1./(obj.d+4)) * fac

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)

>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # 地毯图
>>> kde3 = stats.gaussian_kde(x1, bw_method=my_kde_bandwidth)
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="带宽更小")

我们看到，如果我们将带宽设置得非常窄，那么得到的概率密度函数（PDF）估计仅仅是每个数据点周围的高斯函数之和。

我们现在来看一个更现实的例子，并看看两种可用的带宽选择规则之间的差异。这些规则已知对（接近）正态分布效果良好，即使对于非正态的单峰分布，它们也能合理地工作。作为非正态分布，我们采用 5 个自由度的 t 分布。

import numpy as np
import matplotlib.pyplot as plt
from scipy import statsrng = np.random.default_rng()
x1 = rng.normal(size=200)  # 随机数据，正态分布
xs = np.linspace(x1.min()-1, x1.max()+1, 200)kde1 = stats.gaussian_kde(x1)
kde2 = stats.gaussian_kde(x1, bw_method='silverman')fig = plt.figure(figsize=(8, 6))ax1 = fig.add_subplot(211)
ax1.plot(x1, np.zeros(x1.shape), 'b+', ms=12)  # 地毯图
ax1.plot(xs, kde1(xs), 'k-', label="Scott 法则")
ax1.plot(xs, kde2(xs), 'b-', label="Silverman 法则")
ax1.plot(xs, stats.norm.pdf(xs), 'r--', label="真实 PDF")ax1.set_xlabel('$x$')
ax1.set_ylabel('密度')
ax1.set_title("正态分布（上）和 t 分布（下）")
ax1.legend(loc=1)x2 = stats.t.rvs(5, size=200, random_state=rng)  # 随机数据，t 分布
xs = np.linspace(x2.min() - 1, x2.max() + 1, 200)kde3 = stats.gaussian_kde(x2)
kde4 = stats.gaussian_kde(x2, bw_method='silverman')ax2 = fig.add_subplot(212)
ax2.plot(x2, np.zeros(x2.shape), 'b+', ms=12)  # 地毯图
ax2.plot(xs, kde3(xs), 'k-', label="Scott 法则")
ax2.plot(xs, kde4(xs), 'b-', label="Silverman 法则")
ax2.plot(xs, stats.t.pdf(xs, 5), 'r--', label="真实 PDF")ax2.set_xlabel('$x$')
ax2.set_ylabel('密度')plt.show()

我们现在来看一个双峰分布，其中一个峰较宽，另一个峰较窄。我们预计这将是一个更难近似的密度，因为需要不同的带宽来准确解析每个特征。

>>> from functools import partial

>>> loc1, scale1, size1 = (-2, 1, 175)
>>> loc2, scale2, size2 = (2, 0.2, 50)
>>> x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1),
...                      np.random.normal(loc=loc2, scale=scale2, size=size2)])

>>> x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500)

>>> kde = stats.gaussian_kde(x2)
>>> kde2 = stats.gaussian_kde(x2, bw_method='silverman')
>>> kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2))
>>> kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5))

>>> pdf = stats.norm.pdf
>>> bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \
...               pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size

>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)

>>> ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12)
>>> ax.plot(x_eval, kde(x_eval), 'k-', label="Scott 法则")
>>> ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman 法则")
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2")
>>> ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5")
>>> ax.plot(x_eval, bimodal_pdf, 'r--', label="实际 PDF")

>>> ax.set_xlim([x_eval.min(), x_eval.max()])
>>> ax.legend(loc=2)
>>> ax.set_xlabel('$x$')
>>> ax.set_ylabel('密度')
>>> plt.show()

正如预期的那样，由于双峰分布的两个特征具有不同的特征尺寸，核密度估计（KDE）与真实概率密度函数（PDF）的接近程度并不如我们所期望的那样好。通过将默认带宽减半（Scott * 0.5），我们可以做得稍好一些，而使用比默认值小 5 倍的带宽则不够平滑。然而，在这种情况下，我们真正需要的是非均匀（自适应）带宽。

多变量估计

使用 scipy.stats.gaussian_kde，我们可以进行多变量估计，也可以进行单变量估计。我们以双变量情况为例。首先，我们生成一些随机数据，其中两个变量是相关的。

>>> def measure(n):
...     """测量模型，返回两个相关的测量值。"""
...     m1 = np.random.normal(size=n)
...     m2 = np.random.normal(scale=0.5, size=n)
...     return m1+m2, m1-m2

>>> m1, m2 = measure(2000)
>>> xmin = m1.min()
>>> xmax = m1.max()
>>> ymin = m2.min()
>>> ymax = m2.max()

然后我们将 KDE 应用于数据：

>>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
>>> positions = np.vstack([X.ravel(), Y.ravel()])
>>> values = np.vstack([m1, m2])
>>> kernel = stats.gaussian_kde(values)
>>> Z = np.reshape(kernel.evaluate(positions).T, X.shape)

最后，我们将估计的双变量分布以颜色图的形式绘制出来，并在上面绘制单个数据点。

>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)

>>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
...           extent=[xmin, xmax, ymin, ymax])
>>> ax.plot(m1, m2, 'k.', markersize=2)

>>> ax.set_xlim([xmin, xmax])
>>> ax.set_ylim([ymin, ymax])

风险提示与免责声明
本文内容基于公开信息研究整理，不构成任何形式的投资建议。历史表现不应作为未来收益保证，市场存在不可预见的波动风险。投资者需结合自身财务状况及风险承受能力独立决策，并自行承担交易结果。作者及发布方不对任何依据本文操作导致的损失承担法律责任。市场有风险，投资须谨慎。