找出散点图中最密集区域的中心

Question

我有一组 x 和 y 坐标。我想绘制它们并找出密集绘制的区域。我想找到那个区域中心的坐标。

我已经编写了一个用于在 x-y 平面上绘制点的代码。

import numpy as np
import matplotlib.pyplot as plt
x = np.array([111,152,153,14,155,156,153,154,153,152,154,151,200,201,200])
y = np.array([112,151,153,12,153,154,155,153,152,153,152,153,202,202,204])
plt.scatter(x,y)
plt.show()

但是，我无法找到密集绘制区域的中心。有人可以帮帮我吗？

Answer 1

首先，您需要区分（并确定）密度和聚类。您的问题没有清楚地定义您想要什么，所以我假设您想要评估找到 最大密度 的点之间的位置。老实说，我正准备完成一些令人耳目一新的任务，因此在这里选择了 density 选项;)

单独的术语密度未得到充分说明，因为有多种方法可以评估/估计/计算点云密度。我建议使用 KDE（核密度估计器），因为它很容易适应一个人的需求，并且允许交换核（平方、线性、高斯、余弦以及任何其他适用的）。

注意：完整代码如下；我现在将使用较小的 sn-ps 逐步完成它。

预赛

从您提供的点开始，以下 sn-p 会在您的笛卡尔坐标中计算一个适当的均匀间隔网格：

RESOLUTION = 50
LOCALITY = 2.0

dx = max(pts_x) - min(pts_x)
dy = max(pts_y) - min(pts_y)

delta = min(dx, dy) / RESOLUTION
nx = int(dx / delta)
ny = int(dy / delta)
radius = (1 / LOCALITY) * min(dx, dy)

grid_x = np.linspace(min(pts_x), max(pts_x), num=nx)
grid_y = np.linspace(min(pts_y), max(pts_y), num=ny)

x, y = np.meshgrid(grid_x, grid_y)

如果您愿意，可以通过绘制plt.scatter(grid_x, grid_y) 轻松检查它。所有这些初步计算只是确保手头有所有必需的值：我希望能够为 KDE 内核传播指定某种 resolution 和 locality 设置；此外，我们需要x 和y 方向上的最大点距离来生成网格，计算步长delta，分别计算x 和y 方向上的网格单元数nx 和ny，以及根据我们的位置设置计算内核radius。

核密度估计

为了使用内核估计点密度，我们确实需要两个函数，如下一个 sn-p 中所述

def gauss(x1, x2, y1, y2):
    """
    Apply a Gaussian kernel estimation (2-sigma) to distance between points.

    Effectively, this applies a Gaussian kernel with a fixed radius to one
    of the points and evaluates it at the value of the euclidean distance
    between the two points (x1, y1) and (x2, y2).
    The Gaussian is transformed to roughly (!) yield 1.0 for distance 0 and
    have the 2-sigma located at radius distance.
    """
    return (
        (1.0 / (2.0 * math.pi))
        * math.exp(
            -1 * (3.0 * math.sqrt((x1 - x2)**2 + (y1 - y2)**2) / radius))**2
        / 0.4)


def _kde(x, y):
    """
    Estimate the kernel density at a given position.

    Simply sums up all the Gaussian kernel values towards all points
    (pts_x, pts_y) from position (x, y).
    """
    return sum([
        gauss(x, px, y, py)
        # math.sqrt((x - px)**2 + (y - py)**2)
        for px, py in zip(pts_x, pts_y)
    ])

第一个，gauss 同时执行两项任务：它采用x1, x2, y1, y2 定义的两个点，计算它们的欧几里得距离并使用该距离来评估高斯核的函数值。当距离为0（或非常小）时，高斯核被转换为近似产生1.0，并将其2-sigma 固定在先前计算的radius。此特性由上述locality 设置控制。

确定最大值

幸运的是，numpy 提供了一些简洁的辅助函数来将任意 Python 函数应用于向量和矩阵，因此计算很简单：

kde = np.vectorize(_kde)  # Let numpy care for applying our kde to a vector
z = kde(x, y)

xi, yi = np.where(z == np.amax(z))
max_x = grid_x[xi][0]
max_y = grid_y[yi][0]
print(f"{max_x:.4f}, {max_y:.4f}")

在你的情况下（给定高斯核设置和我的网格假设），最大密度是

155.2041, 154.0800

绘图

您的点云（蓝色十字）与已识别的最大值（红色十字）显示在第一张图片中。第二张图片显示了使用第一个代码 sn-p 中的设置通过高斯 KDE 计算的估计密度。

完整代码

import math
import matplotlib.pyplot as plt
import numpy as np

pts_x = np.array([
    111, 152, 153, 14, 155, 156, 153, 154, 153, 152, 154, 151, 200, 201, 200])
pts_y = np.array([
    112, 151, 153, 12, 153, 154, 155, 153, 152, 153, 152, 153, 202, 202, 204])


RESOLUTION = 50
LOCALITY = 2.0

dx = max(pts_x) - min(pts_x)
dy = max(pts_y) - min(pts_y)

delta = min(dx, dy) / RESOLUTION
nx = int(dx / delta)
ny = int(dy / delta)
radius = (1 / LOCALITY) * min(dx, dy)

grid_x = np.linspace(min(pts_x), max(pts_x), num=nx)
grid_y = np.linspace(min(pts_y), max(pts_y), num=ny)

x, y = np.meshgrid(grid_x, grid_y)


def gauss(x1, x2, y1, y2):
    """
    Apply a Gaussian kernel estimation (2-sigma) to distance between points.

    Effectively, this applies a Gaussian kernel with a fixed radius to one
    of the points and evaluates it at the value of the euclidean distance
    between the two points (x1, y1) and (x2, y2).
    The Gaussian is transformed to roughly (!) yield 1.0 for distance 0 and
    have the 2-sigma located at radius distance.
    """
    return (
        (1.0 / (2.0 * math.pi))
        * math.exp(
            -1 * (3.0 * math.sqrt((x1 - x2)**2 + (y1 - y2)**2) / radius))**2
        / 0.4)


def _kde(x, y):
    """
    Estimate the kernel density at a given position.

    Simply sums up all the Gaussian kernel values towards all points
    (pts_x, pts_y) from position (x, y).
    """
    return sum([
        gauss(x, px, y, py)
        # math.sqrt((x - px)**2 + (y - py)**2)
        for px, py in zip(pts_x, pts_y)
    ])


kde = np.vectorize(_kde)  # Let numpy care for applying our kde to a vector
z = kde(x, y)

xi, yi = np.where(z == np.amax(z))
max_x = grid_x[xi][0]
max_y = grid_y[yi][0]
print(f"{max_x:.4f}, {max_y:.4f}")


fig, ax = plt.subplots()
ax.pcolormesh(x, y, z, cmap='inferno', vmin=np.min(z), vmax=np.max(z))
fig.set_size_inches(4, 4)
fig.savefig('density.png', bbox_inches='tight')

fig, ax = plt.subplots()
ax.scatter(pts_x, pts_y, marker='+', color='blue')
ax.scatter(grid_x[xi], grid_y[yi], marker='+', color='red', s=200)
fig.set_size_inches(4, 4)
fig.savefig('marked.png', bbox_inches='tight')