基于密度(Density-based)的聚类——核密度估计(KDE)

In density-based clustering, clusters are defined as areas of higher density than the remainder of the data set. Objects in sparse areas - that are required to separate clusters-are usually considered to be noise and border points. --wiki

在基于密度的聚类中，聚类定义为密度高于数据集其余部分的区域。稀疏区域中的对象（用于分隔cluster簇）通常被认为是噪声和边界点。

DBSCAN（Density-based spatial clustering of applications with noise带噪声的基于密度的空间聚类应用）与许多更新的方法相比，它具有定义明确的集群模型，称为”密度可达性“，类似于基于链接的聚类，它基于在特定距离阈值内的链接点。然而，它只连接满足密度标准的点，在原始变体中定义为该半径内其他对象的最小数量。聚类(cluster)由所有密度链接的对象加这些对象范围内所有的对象组成。DBSCAN的另一个有趣的特性是它的复杂度很低，它需要在数据库上进行线性数量的范围查询-并且每次运行都会发现基本相同的结果（它对于核心点和噪声点是确定性的，而对于边界点则不是），因此无需多次运行它，OPTICS（ordering points to identify the clustering structure用于标识聚类结构的排序点）是DBSCAN的一般化(generalization)，它消除了为范围参数

选择适当值的需要，并且产生与链接聚类相关的分层结果。

DBSCAN 和OPTICS的主要缺点是它们期望某种密度下降来检测聚类的边界。例如，在具有重叠高斯分布的数据集上，这些算法产生的聚类边界通常看来是任意的，因为聚类密度连续下降。在由高斯混合构成的数据集上，这些算法几乎总是比 诸如EM聚类这样的 能够对此类数据进行精确建模的方法 表现更好。

均值偏移（mean-shift）是一种聚类方法，其中，根据核密度估计，将每个对象移动到其活动中最密集的区域。最终，这些对象会聚到局部的密度最大值。类似于k-means聚类，这些“密度吸引子”可以用作数据集的代表，但均值-漂移可以检测类似于DBSCAN的任意形状的聚类。但是由于费时费力的迭代过程和密度估计，这导致聚类尾部过于碎片化。

在上一节中，我们介绍了高斯混合模型（Gaussian Mixture Model），它是由聚类estimator和密度estimator混合而成。回想一下密度估计器(density estimator)，如果要一个算法从D维数组中产生一个D维概率分布估计，高斯混合模型算法通过将密度表示为高斯分布的加权求和来实现它。从某种意义上说，核密度估计（kernel density estimation）是一种将高斯混合思想发挥到其逻辑极限的算法：它使用由每个点的一个高斯分量组成的混合，从而产生基本上非参数的密度估计量。

密度估计是一种旨在对生成数据集的概率分布进行建模。对于一维数据，我们熟悉的一种简单的density estimator有：直方图(histogram)。直方图将数据分成离散的桶，计算掉在每个桶中的点的个数，然后直观地可视化这个结果。

例如，我们从两个正态分布(normal distribution)中创建一些数据：

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np

def make_data(N, f=0.3, rseed=1):
    rand = np.random.RandomState(rseed)
    x = rand.randn(N)
    x[int(f * N):] += 5
    return x

x = make_data(1000)

可以使用plt.hist()函数创建基于标准计数的直方图。通过指定直方图的normed参数最终得到normalized直方图，其中桶的高度不反应计数，而是反应概率密度。

hist = plt.hist(x, bins=30, normed=True)

需要注意的是，对于均等分桶（equal binning），此种归一化仅改变y轴上的比例，而相对高度和使用计数的方法构建的直方图中的高度基本相同。选择此归一化，以便直方图的总面积等于1，我们可以通过下面的代码来查看直方图函数的输出来确认：

density, bins, patches = hist
widths = bins[1:] - bins[:-1]
(density * widths).sum()

使用直方图作为密度估计器的问题之一是，桶的大小和位置的选择可能会导致表示具有本质上不同的特征。例如，如果我们仅用20个点查看数据的样子，如何绘制这些桶，会导致我们对数据有完全不同的解析！例如下面这个例子：

x = make_data(20)
bins = np.linspace(-5, 10, 10)

fig, ax = plt.subplots(1, 2, figsize=(12, 4),
sharex=True, sharey=True,
subplot_kw={'xlim':(-4, 9),'ylim':(-0.02, 0.3)})
fig.subplots_adjust(wspace=0.05)
for i, offset in enumerate([0.0, 0.6]):
    ax[i].hist(x, bins=bins + offset, normed=True)
    ax[i].plot(x, np.full_like(x, -0.01), '|k',
    markeredgewidth=1)

在左侧，直方图清楚地表明这是双峰分布。在右侧，我们看到一个带有尾巴的几乎均匀分布的模型。如果没有看到前面的代码，你可能不会想到这两个直方图是使用相同的数据构建的。考虑到这一点，我们直观地看上图的直方图不是那么可信的，那么我们如何改进呢？

退一步来讲，我们可以将直方图视为一堆块，在那堆块里，在数据集的每个点上面堆一个块在桶中，我们可以直接到看到：

fig, ax = plt.subplots()
bins = np.arange(-3, 8)
ax.plot(x, np.full_like(x, -0.1), '|k',
markeredgewidth=1)
for count, edge in zip(*np.histogram(x, bins)):
    for i in range(count):
        ax.add_patch(plt.Rectangle((edge, i), 1, 1,alpha=0.5))
ax.set_xlim(-4, 8)
ax.set_ylim(-0.2, 8)

由上图我们看一看到许多块堆在不同的数值范围上，这些桶在两个方向上对其，它并不能够具体地反应出点，如果我们按照点来堆叠，而不是在一个范围类沿着对应的分桶来堆叠块，那么将会是下面这个样子：

x_d = np.linspace(-4, 8, 2000)
density = sum((abs(xi - x_d) < 0.5) for xi in x)
plt.fill_between(x_d, density, alpha=0.5)
plt.plot(x, np.full_like(x, -0.1), '|k', markeredgewidth=1)
plt.axis([-4, 8, -0.2, 8])

结果看起来虽然混乱，但是比标准的直方图更能直观地反应实际的数据特征。尽管如此，粗糙的边缘看起来还是不好看，也不能反应数据的任何真实属性。为了使它们平滑，我们可以使用光滑函数，例如高斯函数，让我们在每个点上使用标准的正态曲线而不是一个块。

from scipy.stats import norm
x_d = np.linspace(-4, 8, 1000)
density = sum(norm(xi).pdf(x_d) for xi in x)
plt.fill_between(x_d, density, alpha=0.5)
plt.plot(x, np.full_like(x, -0.1), '|k', markeredgewidth=1)
plt.axis([-4, 8, -0.2, 5])

经过平滑处理的图像在每个输入点的位置处都有高斯分布，可以更准确地了解数据分布的形状，以及更少的方差（对采样差异的响应变化要小得多）。

最后两幅图是一个维度上的核密度估计(kernel density estimation)的示例: 第一个使用所谓的“顶帽”核 （“tophat” kernel），第二个使用高斯核（Gaussian kernel）我们现在更详细地介绍核密度估计。

核密度估计（kernel density estimation）的自由参数(free parameters)是核(kernel)，核指定每个点的分布形状(shape of the distribution)，核带宽(kernel) 控制每个点的核大小(size of the kernel)。事实上，可以使用多核来进行内核密度估计：特别是，scikit_learn KDE 现实支持6个核（kernel)之一,可以在scikit-learn的density estimation 文档中进行了解。

尽管有多种版本的ptyhon 核密度估计(在scipy和statsmodels软件包中），但我更喜欢使用scikit-learn的版本，因为它更高效和灵活。它是在sklearn.neighbors.KernelDensity中实现的，使用六个核(kernel：gaussian, tophat, epanechnikov, exponential,linear, cosine)之一和几十个距离度量之一处理多维的KDE。由于KDE可能需要大量计算，因此scikit-learn estimation（估计）会在后台使用基于树的算法，并且使用atol（绝对公差）和rtol(相对公差)参数来权衡计算时间和准确性。我们将看到，可以使用scikit-learn的标准交叉验证工具来确定kernel bandwidth 内核带宽。

现在让我们使用scikit-learn 的KernelDensity estimator来重复上面的例子：

from sklearn.neighbors import KernelDensity
 # instantiate and fit the KDE model
kde = KernelDensity(bandwidth=1.0, kernel='gaussian')
kde.fit(x[:, None])
# score_samples returns the log of the probability density
logprob = kde.score_samples(x_d[:, None])
plt.fill_between(x_d, np.exp(logprob), alpha=0.5)
plt.plot(x, np.full_like(x, -0.01), '|k', markeredgewidth=1)
plt.ylim(-0.02, 0.22)

KernelDensity的用法：

sklearn.neighbors.KernelDensity(*, bandwidth=1.0, algorithm='auto', kernel='gaussian', metric='euclidean', atol=0, rtol=0, breadth_first=True, leaf_size=40, metric_params=None)

Bandwidth: float, default=1.0;核的带宽。

Algorithm: {‘kd_tree’, ’ball_tree’, ’auto’}, default=’auto’; 使用的树算法。

Kd_tree:

找到median value （中位数），在该点对它们进行二分类划分。

对上面的点进行划分，找到中位数，x的中位数为6，然后x=6把上面的点分成两部分：

对于上面分类后的左边部分y值为2,3,4,5,6，那么中位数为4，使用y=4这条线，将左边部分划分为两类，如下图：

然后对右边的进行分类，右边y值为1,5,6,7,9,中位数为6，使用y=6划分x=6右边。

然后整棵树的分叉看起来就像下面这个样子：

Ball tree

首先，设置整个数据点云的质心。距质心最大距离的点被选为第一个聚类和子节点的中心。距离第一个聚类中心最远的点被选为第二个聚类的中心点。然后将所有其他数据点分配给该节点，并将其分配到最接近的中心，即cluster 1和cluster2. 任何点只能是一个cluster中的成员，球面线可以相交，但是必须将点明确分配给一个集群。如果一个点恰好位于两个中心的点，并且到两边的距离相同，则必须将其分配给一个cluster。在每个cluster 中重复将数据点分为两个cluster/spheres，直到达到定义的深度。这会导致越来越多的cluster 包含越来越多的圆。

如上图所示，树的深度为2，centroid 1（质心1）是算法的开始。球体（2D）围绕所有数据点（灰色）放置从中心开始，选择cluster中最远的点，此处为3号或9号。这是cluster1的新中心，此处为紫色cluster的3号点。与3号点最远的距离是cluster2的中心。对于上图橙色的cluster, 这个中心点编号为9. 然后考虑紫色球提构成的所有数据点，以计算新的质心2。对橙色球体中的所有数据点进行同样的处理，得出质心3。然后，最远的点再次成为其中心新集群。数据点3是距离质心2以及新cluster中心的最远距离。

在centroid 2（质心点2）的另一侧，最大的距离在质心2与数据点1之间。它是第二个cluster 的中心。然后对橙色那一边也执行此步骤，从而再次产生两个cluster。但是橙色的那一边不平衡。

从下面可以看到生成的树（M是质心1的球体，并且包含所有数据点的起始球体）：

Metric: str, default=’euclidian’; 使用的距离度量方法。

Atol: float, default=0; 所需的结果绝对公差。较大公差会导致执行速度更快。

Rtol: float, default=0; 所需的相对公差，较大公差会导致执行速度更快。

Breadth_first: bool, default=True;如果True, 广度优先; 不然，深度优先。

Leaf_size: int, default=40; 指定基础树的叶大小。

Metric_params: dict, default=None；传递给树用于度量标准的其他参数。

通过交叉验证选择带宽

在KDE中选择带宽对于找到合适的density estimate（密度估计）是非常重要的。它用来控制和权衡密度估计（estimate of density）中的方差和偏差：太窄的带宽会导致一个高的方差估算(i.e., 过拟合），出现或者缺少单个点会对它造成很大的影响。太宽的带宽会导致它有较高的偏差估计（i.e., 欠拟合）,宽的kernel 会将数据结构清除掉。

在机器学习环境中，我们可以看到这种超参数的调试通常是通过交叉验证的方法来完成的。考虑到这一点，scikit-learn中的kernel density estimator(核密度估计）可以直接在scikit-learn’s标准网格搜索中使用。下面我们使用GridSearchCV 来优化之前数据集的带宽(bandwidth)。考虑到我们使用的数据集比较小，我们将会使用留一法 式交叉验证，这样可以最小化每个交叉验证实验的训练集大小的减少。

# from sklearn.grid_search import GridSearchCV
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import LeaveOneOut
bandwidths = 10 ** np.linspace(-1, 1, 100)
grid = GridSearchCV(KernelDensity(kernel='gaussian'),{'bandwidth': bandwidths},cv=LeaveOneOut())
grid.fit(x[:, None])

得到的最好的参数带宽为：1.123。

grid.best_params_

GridsearchCV的使用：

sklearn.model_selection.GridSearchCV(estimator, param_grid, *, scoring=None, n_jobs=None, refit=True, cv=None, verbose=0, pre_dispatch='2*n_jobs', error_score=nan, return_train_score=False)

详尽搜索估计器的指定参数值。通过参数网格进行交叉验证网格搜索来优化用于应用这些方法的估计器（estimator）参数。

GridSearchCV 使用一个“fit”和一个”score” 方法，它也可以使用一些在estimator中可以使用的功能：”score_samples”, ”predict”, ”predict_proba”, ”decision_function”, ”transform” 和 ”inverse_transform”。

Parameters:

Estimator: estimator 对象.

Scikit-learn estimator interface, scikit-learn estimator的接口。estimator要么提供score function，要么传递scoring。

Param_grid: 存参数的dict或list.

Scoring: 评估测试集上交叉验证模型性能的策略。（评分方法）

N_jobs: 并行运行的任务数。

Pre_dispatch: int,or str, default=n_jobs；控制在并行执行期间分派的jobs数量。

Cv: int, cross-validation generator or an iterable, default=None. 交叉验证拆分策略。

Attributes:

Cv_results_: dict of numpy(masked) ndarrays;

可以将key作为列标题和值作为列的字典，将其导入pandas的dataframe.

例如：

Best_estimator_: 通过搜索选择的估算器。

Best_score_: float; 在estimator 为Best_estimator时的平均交叉验证得分。

Best_params_: dic; 在保留数据上给出最佳结果的参数设置

LeaveOneOut的使用：

sklearn.model_selection.LeaveOneOut

Leave-One-Out cross-validator 留一法交叉验证器；

提供训练/测试数据集索引以将数据拆分为训练/测试集。每个样本用作测试集（单例）一次，而其余样本组成训练集。

leaveOneOut() 等同于 KFold(n_splits=n) 和LeavePOut(p=1) ，其中n 为样本数。

由于测试集数量众多（与样本数量相同），因此这种交叉验证方法成本可能比较大。对于大型数据集，应优先使用KFold, ShuffleSplit 或StratifiedKFold.

Methods:

Get_n_splits(X[,y,groups]) : 返回交叉验证中拆分迭代的次数。

Split(X[,y,groups]): 生成索引以将数据拆分成训练和测试数据集。

Example: KDE on a Sphere

下面举个使用KDE来对分布进行可视化的例子。该例子使用了一些可以从scikit learn中加载的数据：两种南美哺乳动物bradypus variegatus（棕喉树懒）和Microryzomys minutus（森林小稻鼠）的观测结果的地理分布。

from sklearn.datasets import fetch_species_distributions
data = fetch_species_distributions()
# Get matrices/arrays of species IDs and locations
latlon = np.vstack([data.train['dd lat'],data.train['dd long']]).T
species = np.array([d.decode('ascii').startswith('micro') for d in data.train['species']], dtype='int')

from mpl_toolkits.basemap import Basemap
from sklearn.datasets.species_distributions import construct_grids
xgrid, ygrid = construct_grids(data)
# plot coastlines with Basemap
m = Basemap(projection='cyl', resolution='c',llcrnrlat=ygrid.min(), urcrnrlat=ygrid.max(),llcrnrlon=xgrid.min(), urcrnrlon=xgrid.max())
m.drawmapboundary(fill_color='#DDEEFF')
m.fillcontinents(color='#FFEEDD')
m.drawcoastlines(color='gray', zorder=2)
m.drawcountries(color='gray', zorder=2)
# plot locations
m.scatter(latlon[:, 1], latlon[:, 0], zorder=3,c=species, cmap='rainbow', latlon=True)

上图并不能很好地说明物种的密度，因为物种范围地点可能会相互重叠。从这个图上面可能不会直观地意识到这一点，但是这个图上面显示了1600多个点。

下面让我们使用核密度估计以一种更可解析的方式来展示这个分布：作为底图上的密度平滑指数。由于此处的坐标系位于球面而非平面上，因此我们将使用Haversine距离测量方式，该测量方式将正确表示曲面上的距离。这里有一些样板代码：

X, Y = np.meshgrid(xgrid[::5], ygrid[::5][::-1])
land_reference = data.coverages[6][::5, ::5]
land_mask = (land_reference > -9999).ravel()
xy = np.vstack([Y.ravel(), X.ravel()]).T
xy = np.radians(xy[land_mask])
# Create two side-by-side plots
fig, ax = plt.subplots(1, 2)
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05)
species_names = ['Bradypus Variegatus', 'Microryzomys Minutus']
cmaps = ['Purples', 'Reds']
for i, axi in enumerate(ax):
     axi.set_title(species_names[i])
     # plot coastlines with Basemap
     m = Basemap(projection='cyl', llcrnrlat=Y.min(),
     urcrnrlat=Y.max(), llcrnrlon=X.min(),
     urcrnrlon=X.max(), resolution='c', ax=axi)
     m.drawmapboundary(fill_color='#DDEEFF')
     m.drawcoastlines()
     m.drawcountries()
     # construct a spherical kernel density estimate of the distribution
     kde = KernelDensity(bandwidth=0.03, metric='haversine')
     kde.fit(np.radians(latlon[species == i]))
     # evaluate only on the land: -9999 indicates ocean
     Z = np.full(land_mask.shape[0], -9999.0)
     Z[land_mask] = np.exp(kde.score_samples(xy))
     Z = Z.reshape(X.shape)
     # plot contours of the density
     levels = np.linspace(0, Z.max(), 25)
     axi.contourf(X, Y, Z, levels=levels, cmap=cmaps[i])

这个图就能够清晰地展示这两个物种的分布状态了。