非参数估计-高斯核平滑Gaussian kernel smoothing-非参数密度估计

阅读原文时间：2021年04月20日阅读：1

1 非参数回归-核平滑

1.1 概念和计算

非参数化回归是指并不需要知道总的分布的情况下进行的一种统计推断回归方法。

核平滑是一种用来估计实值方程的统计方法，其实也就是一种非参数回归。核平滑来作为周围观察数据的加权平均值。权重由核确定，比如越近的数据权重越大。估计方程式平滑的，平滑程度由一个参数控制。

当预测变量的维度小时（p < 3），这个技术是最有效的，比如对于数据可视化。

计算过程如下：

参数意义如下：

令Y(X) 是一个关于X的连续函数。对于每一个X0，Nadaraya-Watson 加权平均值为 (smooth Y(X) estimation)

1.2 Nadaraya-Watson回归

Nadaraya-Watson回归（Nadaraya-Watson regression）是核回归（kernel regression）方法的一种，是使用核函数对数据进行**非参数估计（nonparametric estimation）**的回归模型。

Nadaraya-Watson回归具有非参数模型（nonparametric model）的一般性质。相比于使用参数估计的回归模型，Nadaraya-Watson回归的优势在于其复杂度由数据而不是预先给定的参数决定，且没有学习，即参数估计步骤；其不足在于模型的每次预测都需要处理所有样本，且在高维样本的情形下表现不佳。

1.3 高斯核

高斯核是最常用的核函数之一，表达式如下。

一维数据，通过高斯平滑后的结果如下：

2 高斯核平滑过程-Python实现

2.1 加载库和生成数据

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # Make numpy print 4 significant digits for prettiness
>>> np.set_printoptions(precision=4, suppress=True)
>>> np.random.seed(5) # To get predictable random numbers

>>> n_points = 40
>>> x_vals = np.arange(n_points)
>>> y_vals = np.random.normal(size=n_points)
>>> plt.bar(x_vals, y_vals)

标准的高斯分布如下：

>>> x = np.arange(-6, 6, 0.1) # x from -6 to 6 in steps of 0.1
>>> y = 1 / np.sqrt(2 * np.pi) * np.exp(-x ** 2 / 2.)
>>> plt.plot(x, y)

2.2 Full Width at Half Maximum (FWHM)

在统计上，通常我们通过标准差σ来定义高斯分布的形状。不过，当高斯分布用来做平滑时，我们通常使用**Full Width at Half Maximum (FWHM)**来描述高斯分布的宽度。

FWHM是核的宽度，高斯分布高度最大值的一半的两坐标距离。因此，对于上面的标准高斯分布，最大的高度约为0.4，宽度是0.2。这时，对应的x坐标分别是x = -1.175，和x = 1.175，那么FWHM = 2.35。

FWHM 和 σ 的互换公式如下：

>>> def sigma2fwhm(sigma):
...     return sigma * np.sqrt(8 * np.log(2))

>>> def fwhm2sigma(fwhm):
...     return fwhm / np.sqrt(8 * np.log(2))

我们的例子里：

>>> sigma2fwhm(1)
2.3548200450309493

2.3 分步进行平滑

做平滑的过程很简单。我们一个点一个点的来进行处理。对于每一个数据点，我们都会生成一个新的值。这个值的生成过程如下：

首先我们使用某个函数（在高斯平滑中，这个函数就是我们的高斯分布），根据原先该点的数值和它周围的点的数值，我们生成一些数值。。

比如，对于我们数据集里的第十四个点。我们使用FWHM 为4（x 轴）的高斯分布。为了生成第14个点的新值，我们将高斯分布的中心移动到x =13（第一个点为0）。（为了确保平滑后不进行值的整体缩放，我们将高斯曲线中的值除以曲线下的总面积，使这些值相加为一。）

>>> FWHM = 4
>>> sigma = fwhm2sigma(FWHM)
>>> x_position = 13 # 14th point
>>> kernel_at_pos = np.exp(-(x_vals - x_position) ** 2 / (2 * sigma ** 2))
>>> kernel_at_pos = kernel_at_pos / sum(kernel_at_pos)
>>> plt.bar(x_vals, kernel_at_pos)

高斯分布得到的第12 - 第16个点的值，分别是：

>>> kernel_at_pos[11:16]
array([ 0.1174,  0.1975,  0.2349,  0.1975,  0.1174])

真实数据的第12 - 第16个点的值，分别是：

>>> y_vals[11:16]
array([-0.2049, -0.3588,  0.6035, -1.6648, -0.7002])

然后，我们将高斯核的权重值，和我们原始数据相乘，对结果求和，就得到了我们平滑后的新值（对于第13个点）。

>>> y_by_weight = y_vals * kernel_at_pos # element-wise multiplication
>>> new_val = sum(y_by_weight)
>>> new_val
-0.34796859011845732

对于原始数据的每一个点，我们都重复上述的过程，就得到了高斯平滑后的数据。

>>> for x_position in x_vals:
...     kernel = np.exp(-(x_vals - x_position) ** 2 / (2 * sigma ** 2))
...     kernel = kernel / sum(kernel)
...     smoothed_vals[x_position] = sum(y_vals * kernel)
>>> plt.bar(x_vals, smoothed_vals)

平滑也可以被实现并理解为卷积（convolution）。

2.4 二维平滑

这时高斯分布不再是一个曲线，而是一个圆锥形。

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111, projection='3d')
>>> dx = 0.1
>>> dy = 0.1
>>> x = np.arange(-6, 6, dx)
>>> y = np.arange(-6, 6, dy)
>>> x2d, y2d = np.meshgrid(x, y)
>>> kernel_2d = np.exp(-(x2d ** 2 + y2d ** 2) / (2 * sigma ** 2))
>>> kernel_2d = kernel_2d / (2 * np.pi * sigma ** 2) # unit integral
>>> ax.plot_surface(x2d, y2d, kernel_2d)

和一维数据一样，对于每一个点，我们将kernel的值和原始数据值相乘，来得到这个点平滑之后的值。

2.5 为什么要进行平滑

主要原因是增加信噪比（signal to noise）。

比如下面的例子。

生成的一维信号如下所示：

>>> FWHM = 8
>>> sigma = fwhm2sigma(FWHM)
>>> x_position = 13 # 14th point
>>> sim_signal = np.exp(-(x_vals - x_position) ** 2 / (2 * sigma ** 2))
>>> plt.bar(x_vals, sim_signal)

然后我们向信号中增加噪音：

>>> noise = np.random.normal(size=n_points)
>>> sim_data = sim_signal + noise
>>> plt.bar(x_vals, sim_data)

下面使用高斯平滑

>>> smoothed_sim_data = np.zeros(y_vals.shape)
>>> for x_position in x_vals:
...     kernel = np.exp(-(x_vals - x_position) ** 2 / (2 * sigma ** 2))
...     kernel = kernel / sum(kernel)
...     smoothed_sim_data[x_position] = sum(sim_data * kernel)
>>> plt.bar(x_vals, smoothed_sim_data)