我不久前写了一些代码,使用高斯kde来绘制简单的密度散点图。然而,对于大于100000点的数据集,它只是“永远”运行(几天后我就把它杀死了)。一个朋友在R中给了我一些代码,可以在几秒钟内创建这样一个密度图(plot_fun.R),看起来matplotlib应该也能做同样的事情。

我认为正确的地方是2d直方图,但我正在努力使密度“正确”。我修改了我在this question找到的代码来实现这一点,但是密度没有显示出来,看起来只有密度最大的可能点得到了任何颜色。

下面是我使用的代码:

# initial data
x = -np.log10(np.random.random_sample(10000))
y = -np.log10(np.random.random_sample(10000))
#histogram definition
bins = [1000, 1000] # number of bins
thresh = 3 #density threshold
#data definition
mn = min(x.min(), y.min())
mx = max(x.max(), y.max())
mn = mn-(mn*.1)
mx = mx+(mx*.1)
xyrange = [[mn, mx], [mn, mx]]
# histogram the data
hh, locx, locy = np.histogram2d(x, y, range=xyrange, bins=bins)
posx = np.digitize(x, locx)
posy = np.digitize(y, locy)
#select points within the histogram
ind = (posx > 0) & (posx <= bins[0]) & (posy > 0) & (posy <= bins[1])
hhsub = hh[posx[ind] - 1, posy[ind] - 1] # values of the histogram where the points are
xdat1 = x[ind][hhsub < thresh] # low density points
ydat1 = y[ind][hhsub < thresh]
hh[hh < thresh] = np.nan # fill the areas with low density by NaNs
f, a = plt.subplots(figsize=(12,12))
c = a.imshow(
np.flipud(hh.T), cmap='jet',
extent=np.array(xyrange).flatten(), interpolation='none',
origin='upper'
)
f.colorbar(c, ax=ax, orientation='vertical', shrink=0.75, pad=0.05)
s = a.scatter(
xdat1, ydat1, color='darkblue', edgecolor='', label=None,
picker=True, zorder=2
)

会产生这样的情节:

KDE代码在这里:

^{pr2}$

会产生这样的情节:

当然,问题是,这段代码在大型数据集上不可用。

我的问题是:如何使用2d直方图来生成这样的散点图?ax.hist2d没有产生有用的输出,因为它会给整个绘图着色,而我所有的努力都失败了,我要让上面的2d直方图数据正确地为图的密集区域着色,结果总是要么没有着色,要么只有一小部分最密集的点被着色。显然,我只是不太理解代码。