【Python】机器学习笔记11-核密度估计（Kernel Density Estimation）

本文的参考资料：《Python数据科学手册》；
本文的源代上传到了Gitee上；

本文用到的包：

%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.geoaxes import GeoAxesSubplot

from sklearn.base import BaseEstimator, ClassifierMixin
from sklearn.datasets import fetch_species_distributions, load_digits
from sklearn.model_selection import GridSearchCV, LeaveOneOut, train_test_split
from sklearn.neighbors import KernelDensity

sns.set()
plt.rc('font', family='SimHei')
plt.rc('axes', unicode_minus=False)

核密度估计（Kernel Density Estimation）

密度评估器是一种利用D维数据集生成D维概率分布估计的算法，GMM使用不同的高斯分布的加权汇总来表示概率分布估计。
核密度估计算法将高斯混合理念扩展到了逻辑极限，它通过对每一个点生成高斯分布的混合成分，获得本实质上是无参数的密度评估器（书本原话）。

核密度估计的自由参数是核类型和核带宽，前者指定每个点核密度分布的形状，后者指定每个点核的大小。

上边是书上讲的话，下面是我的理解：

核密度估计本质上是从有限的样本中尽可能地估算概率密度函数，类似于让直方图的区间趋向无穷小后得到的图形，当然由于样本数量是有限的不能直接实现这样的操作。

核密度估计示例

核密度估计在sklearn中由KernelDensity类实现，kernel参数指定需要用到的核函数，来适应不同的需求；
以下代码生成一组高斯分布的随机数据，分别绘制其直方图与KDE图，查看效果；

x_train = np.hstack((np.random.normal(2, 1.66, 200), np.random.normal(8, 2.33, 200)))  # 两个高斯分布组合到一起

model = KernelDensity(bandwidth=1.0, kernel='gaussian')
model.fit(x_train[:, np.newaxis])
x_range = np.linspace(x_train.min() - 1, x_train.max() + 1, 500)
x_log_prob = model.score_samples(x_range[:, np.newaxis])  # 这个东西返回概率的对数
x_prob = np.exp(x_log_prob)

plt.figure(figsize=(10, 10))
r = plt.hist(
    x=x_train,
    bins=50,
    density=True,
    histtype='stepfilled',
    color='red',
    alpha=0.5,
    label='直方图',
)
plt.fill_between(
    x=x_range,
    y1=x_prob,
    y2=0,
    color='green',
    alpha=0.5,
    label='KDE',
)
plt.plot(x_range, x_prob, color='gray')
plt.vlines(x=2, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.vlines(x=8, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.ylim(0, r[0].max() + 0.011)
plt.legend(loc='upper right')
plt.title('同一组数据的直方图与KDE图')

使用交叉检验确定带宽

在使用核密度估计时，如果带宽设置过小，会出现过拟合的现象，如果带宽设置过大，会出现欠拟合的现象，因此需要确定好最佳的带宽；
sklearn中的KernelDensity类支持使用GridSearchCV来寻找最佳参数，以下是一个示例：

x_train = np.hstack((np.random.normal(2, 1.66, 200), np.random.normal(8, 2.33, 200)))

grid = GridSearchCV(
    estimator=KernelDensity(kernel='gaussian'),
    param_grid={
   'bandwidth': 10 ** np.linspace(-1, 1, 100)},
    cv=LeaveOneOut(),
)
grid.fit(x_train[:, np.newaxis])
print(f'最佳带宽：{grid.best_params_["bandwidth"]}')

输出结果：

最佳带宽：0.7742636826811272

在球形空间中使用KDE

从sklearn中加载物种分布数据（fetch_species_distributions函数），返回的数据中包含了森林小稻鼠和褐喉树懒在南美中的分布数据；
我们将用KDE从这些分布数据中计算这两个物种的分布密度，并显示在地图上（还真是在球面上的数据）；

由于书本介绍的Basemap已经停止维护了，所以这里我是用了cartopy这个库来显示地图（费了不少力气）；
书本代码用到了一个在如今的sklearn中已经被废弃的函数construct_grids，我从网上找到了它的源码，凑合用一下；

def construct_grids(batch):  # 这个函数在现在版本的sklearn中被废弃了，从网上找到了以前的源码
    """Construct the map grid from the batch object

    Parameters
    ----------
    batch : Batch object
        The object returned by :func:`fetch_species_distributions`

    Returns
    -------
    (xgrid, ygrid) : 1-D arrays
        The grid corresponding to the values in batch.coverages
    """
    # x,y coordinates for corner cells
    xmin = batch.x_left_lower_corner + batch.grid_size
    xmax = xmin + (batch.Nx * batch.grid_size)
    ymin = batch.y_left_lower_corner + batch.grid_size
    ymax = ymin + (batch.Ny * batch.grid_size)

    # x coordinates of the grid cells
    xgrid = np.arange(xmin, xmax, batch.grid_size)
    # y coordinates of the grid cells
    ygrid = np.arange(ymin, ymax, batch.grid_size)

    return xgrid, ygrid

species_data = fetch_species_distributions()
species_names = ['褐喉树懒', '森林小稻鼠']

lat_and_lon = np.vstack((
    species_data.train['dd lat'],
    species_data.train['dd long'],
)).T

species_train = np.array