Tạo một ô mật độ cho các vùng có mật độ cao, điểm cho các vùng thưa thớt

Question

Tôi cần tạo một ô có chức năng giống như ô mật độ cho các vùng có mật độ cao trên ô, nhưng dưới một số ngưỡng sử dụng các điểm riêng lẻ. Tôi không thể tìm thấy bất kỳ mã hiện có nào trông giống như những gì tôi cần trong thư viện hình thu nhỏ matplotlib hoặc từ tìm kiếm trên google. Tôi có một mã số làm việc tôi đã viết bản thân mình, nhưng nó là hơi phức tạp và (quan trọng hơn) mất một thời gian dài không thể chấp nhận được khi số lượng các điểm/thùng là lớn. Đây là mã:

import numpy as np 
import math 
import matplotlib as mpl 
import matplotlib.pyplot as plt 
import pylab 
import numpy.random 

#Create the colormap: 
halfpurples = {'blue': [(0.0,1.0,1.0),(0.000001, 0.78431373834609985, 0.78431373834609985), 
(0.25, 0.729411780834198, 0.729411780834198), (0.5, 
0.63921570777893066, 0.63921570777893066), (0.75, 
0.56078433990478516, 0.56078433990478516), (1.0, 0.49019607901573181, 
0.49019607901573181)], 

    'green': [(0.0,1.0,1.0),(0.000001, 
    0.60392159223556519, 0.60392159223556519), (0.25, 
    0.49019607901573181, 0.49019607901573181), (0.5, 
    0.31764706969261169, 0.31764706969261169), (0.75, 
    0.15294118225574493, 0.15294118225574493), (1.0, 0.0, 0.0)], 

    'red': [(0.0,1.0,1.0),(0.000001, 
    0.61960786581039429, 0.61960786581039429), (0.25, 
    0.50196081399917603, 0.50196081399917603), (0.5, 
    0.41568627953529358, 0.41568627953529358), (0.75, 
    0.32941177487373352, 0.32941177487373352), (1.0, 
    0.24705882370471954, 0.24705882370471954)]} 

halfpurplecmap = mpl.colors.LinearSegmentedColormap('halfpurples',halfpurples,256) 

#Create x,y arrays of normally distributed points 
npts = 1000 
x = numpy.random.standard_normal(npts) 
y = numpy.random.standard_normal(npts) 

#Set bin numbers in both axes 
nxbins = 25 
nybins = 25 

#Set the cutoff for resolving the individual points 
minperbin = 1 

#Make the density histrogram 
H, yedges, xedges = np.histogram2d(y,x,bins=(nybins,nxbins)) 
#Reorient the axes 
H = H[::-1] 

extent = [xedges[0],xedges[-1],yedges[0],yedges[-1]] 

#Compute all bins where the density plot value is below (or equal to) the threshold 
lowxleftedges = [[xedges[i] for j in range(len(H[:,i])) if H[j,i] <= minperbin] for i in range(len(H[0,:]))] 
lowxrightedges = [[xedges[i+1] for j in range(len(H[:,i])) if H[j,i] <= minperbin] for i in range(len(H[0,:]))] 
lowyleftedges = [[yedges[-(j+2)] for j in range(len(H[:,i])) if H[j,i] <= minperbin] for i in range(len(H[0,:]))] 
lowyrightedges = [[yedges[-(j+1)] for j in range(len(H[:,i])) if H[j,i] <= minperbin] for i in range(len(H[0,:]))] 

#Flatten and convert to numpy array 
lowxleftedges = np.asarray([item for sublist in lowxleftedges for item in sublist]) 
lowxrightedges = np.asarray([item for sublist in lowxrightedges for item in sublist]) 
lowyleftedges = np.asarray([item for sublist in lowyleftedges for item in sublist]) 
lowyrightedges = np.asarray([item for sublist in lowyrightedges for item in sublist]) 

#Find all points that lie in these regions 
lowdatax = [[x[i] for j in range(len(lowxleftedges)) if lowxleftedges[j] <= x[i] and x[i] <= lowxrightedges[j] and lowyleftedges[j] <= y[i] and y[i] <= lowyrightedges[j]] for i in range(len(x))] 
lowdatay = [[y[i] for j in range(len(lowyleftedges)) if lowxleftedges[j] <= x[i] and x[i] <= lowxrightedges[j] and lowyleftedges[j] <= y[i] and y[i] <= lowyrightedges[j]] for i in range(len(y))] 

#Flatten and convert into numpy array 
lowdatax = np.asarray([item for sublist in lowdatax for item in sublist]) 
lowdatay = np.asarray([item for sublist in lowdatay for item in sublist]) 

#Plot 
fig1 = plt.figure() 
ax1 = fig1.add_subplot(111) 
ax1.plot(lowdatax,lowdatay,linestyle='.',marker='o',mfc='k',mec='k') 
cp1 = ax1.imshow(H,interpolation='nearest',extent=extent,cmap=halfpurplecmap,vmin=minperbin) 
fig1.colorbar(cp1) 

fig1.savefig('contourtest.eps') 

Mã này tạo ra một hình ảnh trông như thế này:

Tuy nhiên, khi sử dụng trên dữ liệu lớn hơn thiết lập chương trình phải mất vài giây đến vài phút. Bất kỳ suy nghĩ về làm thế nào để tăng tốc độ này lên? Cảm ơn!

Answer 1

này nên làm điều đó:

import matplotlib.pyplot as plt, numpy as np, numpy.random, scipy 

#histogram definition 
xyrange = [[-5,5],[-5,5]] # data range 
bins = [100,100] # number of bins 
thresh = 3 #density threshold 

#data definition 
N = 1e5; 
xdat, ydat = np.random.normal(size=N), np.random.normal(1, 0.6, size=N) 

# histogram the data 
hh, locx, locy = scipy.histogram2d(xdat, ydat, range=xyrange, bins=bins) 
posx = np.digitize(xdat, locx) 
posy = np.digitize(ydat, 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 = xdat[ind][hhsub < thresh] # low density points 
ydat1 = ydat[ind][hhsub < thresh] 
hh[hh < thresh] = np.nan # fill the areas with low density by NaNs 

plt.imshow(np.flipud(hh.T),cmap='jet',extent=np.array(xyrange).flatten(), interpolation='none', origin='upper') 
plt.colorbar() 
plt.plot(xdat1, ydat1, '.',color='darkblue') 
plt.show() 

Answer 2

Vấn đề của bạn là bậc hai - đối với npts = 1000, bạn có kích thước mảng đạt 10^6 điểm và hơn bạn lặp qua các danh sách này với tính năng hiểu danh sách.
Bây giờ, đây là vấn đề về hương vị của khóa học, nhưng tôi thấy rằng việc hiểu danh sách có thể mang lại một mã hoàn toàn khó thực hiện, và đôi khi chúng hơi nhanh hơn một chút ... nhưng đó không phải là quan điểm của tôi.
Quan điểm của tôi là cho các hoạt động mảng lớn bạn có chức năng NumPy như:

np.where, np.choose etc. 

Thấy rằng bạn có thể đạt được điều đó chức năng của comprehensions danh sách với NumPy, và mã của bạn sẽ chạy nhanh hơn.

Tôi có hiểu chính xác, nhận xét của bạn không?

#Find all points that lie in these regions 

bạn đang thử nghiệm một điểm trong một đa giác? nếu có, hãy xem xét point in polygon bên trong matplotlib.

Answer 3

Sau một đêm để ngủ và đọc qua đề xuất của Oz123, tôi đã tìm ra. Bí quyết là tính toán bin nào mỗi x, y rơi vào (xi, yi), sau đó kiểm tra nếu H [xi, yi] (thực ra, trong trường hợp của tôi H [yi, xi]) nằm dưới ngưỡng. Đoạn mã dưới, và chạy rất nhanh cho một số lượng lớn các điểm và có nhiều bụi:

import numpy as np 
import math 
import matplotlib as mpl 
import matplotlib.pyplot as plt 
import pylab 
import numpy.random 

#Create the colormap: 
halfpurples = {'blue': [(0.0,1.0,1.0),(0.000001, 0.78431373834609985, 0.78431373834609985), 
0.25, 0.729411780834198, 0.729411780834198), (0.5, 
0.63921570777893066, 0.63921570777893066), (0.75, 
0.56078433990478516, 0.56078433990478516), (1.0, 0.49019607901573181, 
0.49019607901573181)], 

    'green': [(0.0,1.0,1.0),(0.000001, 
    0.60392159223556519, 0.60392159223556519), (0.25, 
    0.49019607901573181, 0.49019607901573181), (0.5, 
    0.31764706969261169, 0.31764706969261169), (0.75, 
    0.15294118225574493, 0.15294118225574493), (1.0, 0.0, 0.0)], 

    'red': [(0.0,1.0,1.0),(0.000001, 
    0.61960786581039429, 0.61960786581039429), (0.25, 
    0.50196081399917603, 0.50196081399917603), (0.5, 
    0.41568627953529358, 0.41568627953529358), (0.75, 
    0.32941177487373352, 0.32941177487373352), (1.0, 
    0.24705882370471954, 0.24705882370471954)]} 

halfpurplecmap = mpl.colors.LinearSegmentedColormap('halfpurples',halfpurples,256) 

#Create x,y arrays of normally distributed points 
npts = 100000 
x = numpy.random.standard_normal(npts) 
y = numpy.random.standard_normal(npts) 

#Set bin numbers in both axes 
nxbins = 100 
nybins = 100 

#Set the cutoff for resolving the individual points 
minperbin = 1 

#Make the density histrogram 
H, yedges, xedges = np.histogram2d(y,x,bins=(nybins,nxbins)) 
#Reorient the axes 
H = H[::-1] 

extent = [xedges[0],xedges[-1],yedges[0],yedges[-1]] 

#Figure out which bin each x,y point is in 
xbinsize = xedges[1]-xedges[0] 
ybinsize = yedges[1]-yedges[0] 
xi = ((x-xedges[0])/xbinsize).astype(np.integer) 
yi = nybins-1-((y-yedges[0])/ybinsize).astype(np.integer) 

#Subtract one from any points exactly on the right and upper edges of the region 
xim1 = xi-1 
yim1 = yi-1 
xi = np.where(xi < nxbins,xi,xim1) 
yi = np.where(yi < nybins,yi,yim1) 

#Get all points with density below the threshold 
lowdensityx = x[H[yi,xi] <= minperbin] 
lowdensityy = y[H[yi,xi] <= minperbin] 

#Plot 
fig1 = plt.figure() 
ax1 = fig1.add_subplot(111) 
ax1.plot(lowdensityx,lowdensityy,linestyle='.',marker='o',mfc='k',mec='k',ms=3) 
cp1 = ax1.imshow(H,interpolation='nearest',extent=extent,cmap=halfpurplecmap,vmin=minperbin) 
fig1.colorbar(cp1) 

fig1.savefig('contourtest.eps') 

Answer 4

Đối với hồ sơ, đây là kết quả của một nỗ lực mới sử dụng scipy.stats.gaussian_kde thay vì một biểu đồ 2D. Người ta có thể hình dung các kết hợp khác nhau của chia lưới màu và đường viền tùy theo mục đích.

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

# parameters 
npts = 5000   # number of sample points 
bins = 100   # number of bins in density maps 
threshold = 0.01 # density threshold for scatter plot 

# initialize figure 
fig, ax = plt.subplots() 

# create a random dataset 
x1, y1 = np.random.multivariate_normal([0, 0], [[1, 0], [0, 1]], npts/2).T 
x2, y2 = np.random.multivariate_normal([4, 4], [[4, 0], [0, 1]], npts/2).T 
x = np.hstack((x1, x2)) 
y = np.hstack((y1, y2)) 
points = np.vstack([x, y]) 

# perform kernel density estimate 
kde = gaussian_kde(points) 
z = kde(points) 

# mask points above density threshold 
x = np.ma.masked_where(z > threshold, x) 
y = np.ma.masked_where(z > threshold, y) 

# plot unmasked points 
ax.scatter(x, y, c='black', marker='.') 

# get bounds from axes 
xmin, xmax = ax.get_xlim() 
ymin, ymax = ax.get_ylim() 

# prepare grid for density map 
xedges = np.linspace(xmin, xmax, bins) 
yedges = np.linspace(ymin, ymax, bins) 
xx, yy = np.meshgrid(xedges, yedges) 
gridpoints = np.array([xx.ravel(), yy.ravel()]) 

# compute density map 
zz = np.reshape(kde(gridpoints), xx.shape) 

# plot density map 
im = ax.imshow(zz, cmap='CMRmap_r', interpolation='nearest', 
       origin='lower', extent=[xmin, xmax, ymin, ymax]) 

# plot threshold contour 
cs = ax.contour(xx, yy, zz, levels=[threshold], colors='black') 

# show 
fig.colorbar(im) 
plt.show()