python - Matplotlib 中的散点图轮廓

Question

我在 matplotlib 中生成了一个巨大的散点图（约 100,000 个点）。每个点在这个 x/y 空间中都有一个位置，我想生成包含点总数的某些百分位数的轮廓。

matplotlib 中是否有可以执行此操作的函数？我已经研究过轮廓（），但我必须编写自己的函数才能以这种方式工作。

谢谢！

Answer 1

基本上，您需要某种密度估计。有多种方法可以做到这一点：

1. 使用某种二维直方图（例如matplotlib.pyplot.hist2d或matplotlib.pyplot.hexbin）（您也可以将结果显示为等高线 - 只需使用numpy.histogram2d然后对结果数组进行等高线处理。）

2. 进行核密度估计 (KDE) 并对结果进行轮廓化。KDE 本质上是一个平滑的直方图。它不会将一个点落入特定的 bin 中，而是为周围的 bin 添加权重（通常呈高斯“钟形曲线”的形状）。

使用 2D 直方图简单易懂，但从根本上给出了“块状”结果。

“正确”地做第二个有一些皱纹（即没有一种正确的方法）。我不会在这里详细介绍，但是如果您想从统计角度解释结果，则需要阅读它（尤其是带宽选择）。

无论如何，这里有一个差异的例子。我将类似地绘制每一个，所以我不会使用等高线，但您可以使用等高线图轻松绘制 2D 直方图或高斯 KDE：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import kde

np.random.seed(1977)

# Generate 200 correlated x,y points
data = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 3]], 200)
x, y = data.T

nbins = 20

fig, axes = plt.subplots(ncols=2, nrows=2, sharex=True, sharey=True)

axes[0, 0].set_title('Scatterplot')
axes[0, 0].plot(x, y, 'ko')

axes[0, 1].set_title('Hexbin plot')
axes[0, 1].hexbin(x, y, gridsize=nbins)

axes[1, 0].set_title('2D Histogram')
axes[1, 0].hist2d(x, y, bins=nbins)

# Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
k = kde.gaussian_kde(data.T)
xi, yi = np.mgrid[x.min():x.max():nbins*1j, y.min():y.max():nbins*1j]
zi = k(np.vstack([xi.flatten(), yi.flatten()]))

axes[1, 1].set_title('Gaussian KDE')
axes[1, 1].pcolormesh(xi, yi, zi.reshape(xi.shape))

fig.tight_layout()
plt.show()

一个警告：点数非常多，scipy.stats.gaussian_kde会变得非常慢。通过近似来加速它是相当容易的——只需获取 2D 直方图并使用正确半径和协方差的高斯滤波器对其进行模糊处理。如果你愿意，我可以举个例子。

另一个警告：如果您在非笛卡尔坐标系中执行此操作，则这些方法都不适用！在球壳上获得密度估计要复杂一些。

Answer 2

我也有同样的问题。如果要绘制包含部分点的等高线，可以使用以下算法：

创建二维直方图

h2, xedges, yedges = np.histogram2d(X, Y, bibs = [30, 30])

h2 现在是包含整数的二维矩阵，整数是某个矩形中的点数

hravel = np.sort(np.ravel(h2))[-1] #all possible cases for rectangles 
hcumsum = np.sumsum(hravel)

丑陋的黑客，

让我们为 h2 2d 矩阵中的每个点给出矩形的累积点数，其中包含的点数等于或大于我们当前分析的点数。

hunique = np.unique(hravel)

hsum = np.sum(h2)

for h in hunique:
    h2[h2 == h] = hcumsum[np.argwhere(hravel == h)[-1]]/hsum

现在绘制h2的轮廓，它将是包含一些所有点的轮廓