Python 数学建模——高斯核密度估计

前言

  高斯核密度估计本是一种机器学习算法，在数学建模中也可以发挥作用。本文主要讨论用它来拟合变量的概率密度，获得概率密度函数 $f(x)$。

原理

  已知一个连续型随机变量 $X$ 的一系列观测值 $X_1,X_2,\cdots ,X_n$，我们可以用高斯核密度估计来拟合出 $X$ 的概率密度函数 $f(x)$。

代码实例

scipy 实现

  先给出代码，后面慢慢解释。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde# 为测试 gaussian_kde 效果，当场生成1000个beta分布的随机数
np.random.seed(114)
info = np.random.beta(a=2, b=5, size=1000)# 根据一系列观测值 info，拟合出概率密度
# 这个 gaussian_kde 有个神奇的参数 bw_method，说是计算估计器带宽的方法，可以调一下
kde = gaussian_kde(info)
x = np.linspace(min(info), max(info), 1000)
pdf = kde.evaluate(x)# 开始作图# 指定楷体以显示中文字体
plt.rcParams["font.sans-serif"] = ['KaiTi']
fig = plt.figure()
ax1 = fig.add_subplot(111, label="1")
ax2 = fig.add_subplot(111, label="2", frame_on=False)# 100 指的是 100 个直方条
ax1.hist(info,100, color='r', alpha=0.4)
ax1.set_xlabel("观测值")
ax1.set_ylabel("观测频数")ax2.plot(x, pdf, c = 'b')
ax2.set_xticks([])
ax2.set_ylabel("拟合概率密度")
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position('right')
plt.show()

  图画出来是这样的：

• 第 $6-7$ 行，按照 $\beta$ 分布生成了一个含有 $1000$ 个随机数的info。实际建模的时候，这个info应该是题目给我们的数据样本。
• 第 $11$ 行由info得到一个kde对象，kde.evaluate()实际上就是概率密度函数，传一个 $x$ 就返回一个 $f(x)$。
• 第 $12-13$ 行是取了横轴上的一些数据点获取概率密度用于后续画图，如果把x记作 $[x_1,\cdots ,x_n]$ 那么pdf就是 $[f(x_1),\cdots ,f(x_n)]$。
• 第 $18$ 行到最后都是在画图。

  如果你觉得拟合效果不佳，可以调整gaussian_kde的bw_method参数。这个东西是采样宽度，换句话说它越大 $f(x)$ 越粗糙，它越小 $f(x)$ 更容易过拟合。比如我在上面的kde = gaussian_kde(info)中加入参数bw_method=0.1，画出来这样的图：

  明显有些过拟合，但是很多时候我们需要这种过拟合。

seaborn 实现

  Seaborn 作为一个强大的 Python 可视化库，也内置了高斯核密度估计的功能。Seaborn 只需要一行代码即可画出核密度估计图，但是它无法返回pdf或者kde等对象，也就是说我们只能看到 $y=f(x)$ 的大致图像，却无法获取 $f(x_0)$ （$x_0$ 是某个具体值，比如 $x_0=5$）的值。如果你的目的只是为了可视化，展示随机变量 $X$ 的集中程度和均值情况，那么 seaborn 无疑是更方便的选择。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 为测试 gaussian_kde 效果，当场生成1000个beta分布的随机数
np.random.seed(114)
info = np.random.beta(a=2, b=5, size=1000)# seaborn 作图
sns.set_style("whitegrid") # 风格设置
sns.kdeplot(info,shade=True, color="g") # 概率密度函数底下填充绿色阴影
plt.show()

  画出来效果如下，该图象轮廓与上面的应该是一样的。如果有两个分布类似的样本，使用 seaborn 用不同颜色画出概率密度函数，将会对比鲜明，非常好看。

  除了kdeplot，seaborn 库里的distplot也可以进行高斯核密度估计，优点是它还能带上直方图。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 为测试 gaussian_kde 效果，当场生成1000个beta分布的随机数
np.random.seed(114)
info = np.random.beta(a=2, b=5, size=1000)# seaborn 作图
sns.set_style("whitegrid") # 风格设置
sns.distplot(info)
plt.show()

  你也可以调整distplot的参数，hist=False不画直方图，kde=False不画概率密度函数。