基于Python的2D散点图自制高斯KDE密度可视化问题

首先，我们需要调整你的KDE函数来计算二维网格上的联合密度值，再通过插值或网格匹配的方式，把这些密度值映射回每个原始(x,y)数据点。下面是具体的步骤和代码：

1. 修正KDE函数，计算二维网格密度 #

你的原始函数计算的是x和y方向分开的权重，但我们需要的是二维联合密度。调整后的函数会生成一个与xGrid、yGrid对应的二维密度矩阵：

import numpy as np
from matplotlib import pyplot as plt
from scipy.interpolate import interpn

def homemadeKDE(x, xgrid, y, ygrid, sigmaX=1, sigmaY=1):
    # 计算每个x网格点对所有x数据的高斯核贡献
    x_kernels = np.exp(-((xgrid[:, None] - x) / (2 * sigmaX)) ** 2)
    # 计算每个y网格点对所有y数据的高斯核贡献
    y_kernels = np.exp(-((ygrid[:, None] - y) / (2 * sigmaY)) ** 2)
    
    # 二维联合密度：核的外积之和，再归一化（除以数据量和核的归一化因子）
    # 这里的归一化因子sigmaX*sigmaY*2π是高斯核的积分，确保密度积分接近1
    density_grid = np.dot(x_kernels, y_kernels.T) / (len(x) * sigmaX * sigmaY * 2 * np.pi)
    return density_grid

2. 生成数据并计算网格密度 #

注意：你原始代码中xGrid的范围是0-500，但x数据是0-1的随机数，这会导致大部分网格点没有数据。我们调整网格范围到数据的实际区间，让密度计算更合理：

# 生成原始数据
np.random.seed(42)  # 设置随机种子，结果可复现
x = np.random.rand(10000)
y = np.random.rand(10000)

# 生成匹配数据范围的网格
xGrid = np.linspace(x.min(), x.max(), 50)  # x方向50个网格点
yGrid = np.linspace(y.min(), y.max(), 50)  # y方向50个网格点

# 计算二维网格密度（根据数据尺度调整带宽）
density_grid = homemadeKDE(x, xGrid, y, yGrid, sigmaX=0.05, sigmaY=0.05)

3. 将网格密度映射到原始数据点 #

我们使用scipy.interpolate.interpn函数，通过插值得到每个原始(x,y)点对应的密度值。这种方法比直接匹配网格点更平滑：

# 准备网格坐标用于插值
grid_coords = (xGrid, yGrid)

# 对每个原始(x,y)点插值得到密度z值
z = interpn(grid_coords, density_grid, np.column_stack((x, y)), method='linear')

4. 绘制带密度颜色的散点图 #

现在可以用plt.scatter绘制散点图，用z值作为颜色：

plt.figure(figsize=(10, 8))
scatter = plt.scatter(x, y, c=z, cmap='jet', s=5, alpha=0.6)
plt.colorbar(scatter, label='Density')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scatter Plot with Homemade KDE Density Coloring')
plt.show()

关键说明 #

• 核带宽调整：sigmaX和sigmaY需要根据你的数据尺度调整。如果数据范围大，带宽可以适当增大；反之减小。比如上面的例子中x和y都是0-1，所以用0.05作为带宽。
• 效率优化：使用粗网格KDE的好处是降低计算复杂度（从O(n²)降到O(n*m)，m是网格点数量），非常适合大规模数据集。
• 插值方法：interpn支持多种插值方法，比如linear（线性插值）、nearest（最近邻），你可以根据需求选择。如果追求速度，nearest会更快。

内容的提问来源于stack exchange，提问作者Fourier