使用Python进行二维散点图高斯求和

Question

我正在尝试建立一个人们可能宽泛称为自制KDE的东西。 我正在尝试评估相当大的数据点集的密度。 特别是，对于散点图中有许多数据点，我想使用颜色渐变来指示密度（请参见以下链接）。

为了举例说明，我提供了下面的一对随机（x，y）数据。 实际数据将分布在不同的比例上，因此X和Y网格点间距会有所不同。

import numpy as np
from matplotlib import pyplot as plt

def homemadeKDE(x, xgrid, y, ygrid, sigmaX = 1, sigmaY = 1):

    a = np.exp( -((xgrid[:,None]-x)/(2*sigmaX))**2 )
    b = np.exp( -((ygrid[:,None]-y)/(2*sigmaY))**2 ) 

    xweights = np.dot(a, x.T)/np.sum(a)
    yweights = np.dot(b, y.T)/np.sum(b)  

    return xweights, yweights

x = np.random.rand(10000)
x.sort()
y = np.random.rand(10000)

xGrid = np.linspace(0, 500, 501)
yGrid = np.linspace(0, 10, 11)

newX, newY = homemadeKDE(x, xGrid, y, yGrid)

我卡住的问题是如何将这些值投影回原始的x和y向量，以便我可以用它来绘制一个二维散点图(x,y)，其中z值表示密度，并根据给定的颜色映射进行着色，代码如下：

plt.scatter(x, y, c = z, cmap = "jet")

绘图和KDE方法实际上是受到这个很棒的回答的启发：answer

编辑1 为了消除一些困惑，想法是进行高斯KDE，这将在一个更粗糙的栅格上进行。 SigmaX和sigmaY分别反映了x和y方向核带宽。

- Fourier

---

权重仅是基于网格的点加权，遵循此定义：http://mccormickml.com/2014/02/26/kernel-regression/ - Fourier

@Gabriel 这类似于基于连续高斯函数的分箱对点进行着色向量。请查看我在链接中提供的基于scipy.stats的解决方案。 - Fourier

两个不同维度中的kde与您所瞄准的2D kde有任何关系的假设是从哪里来的？这难道不是一个数学问题吗？ - ImportanceOfBeingErnest

@Gabriel 没什么。我已经使用它，运行得非常完美。不过，我想写自己的代码来更清楚地理解它。 - Fourier

也许这会有所帮助，看一下如何生成多元KDE：http://sebastianraschka.com/Articles/2014_kernel_density_est.html#322-implementing-the-code-to-calculate-the-multivariate-gaussian-densities - Gabriel

Answer 1

我实际上 - 经过一点思考 - 能够自己解决问题。同时也感谢帮助和富有洞见的评论。

import numpy as np
from matplotlib import pyplot as plt

def gaussianSum1D(gridpoints, datapoints, sigma=1):

    a = np.exp( -((gridpoints[:,None]-datapoints)/sigma)**2 )

    return a

#some test data
x = np.random.rand(10000)
y = np.random.rand(10000)

#create grids
gridSize = 100
xedges = np.linspace(np.min(x), np.max(x), gridSize)
yedges = np.linspace(np.min(y), np.max(y), gridSize)

#calculate weights for both dimensions seperately
a = gaussianSum1D(xedges, x, sigma=2)
b = gaussianSum1D(yedges, y, sigma=0.1)

Z = np.dot(a, b.T).T

#plot original data
fig, ax = plt.subplots()
ax.scatter(x, y, s = 1)
#overlay data with contours 
ax.contour(xedges, yedges, Z, cmap = "jet")

- Fourier