2012-05-03 20 views
10

मुझे साजिश पर उच्च घनत्व वाले क्षेत्रों के लिए एक घनत्व साजिश की तरह काम करने की आवश्यकता है, लेकिन कुछ सीमाओं के नीचे व्यक्तिगत बिंदुओं का उपयोग करता है। मुझे कोई मौजूदा कोड नहीं मिला जो मुझे matplotlib थंबनेल गैलरी या Google खोजों से जो चाहिए उसे समान दिखता है। मेरे पास एक कामकाजी कोड है जिसे मैंने स्वयं लिखा है, लेकिन यह कुछ हद तक मुश्किल है और (अधिक महत्वपूर्ण बात यह है कि अंक/डिब्बे बड़ी संख्या में एक अस्वीकार्य रूप से लंबा समय लगता है। यहाँ कोड है:उच्च घनत्व वाले क्षेत्रों के लिए कुशलता से घनत्व साजिश बनाएं, स्पैस क्षेत्रों के लिए अंक

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') 

इस कोड को एक छवि है कि इस तरह दिखता है पैदा करता है:

countour test

हालांकि, जब बड़े डेटा पर इस्तेमाल किया सेट कार्यक्रम मिनट के लिए कई सेकंड लेता है। इसे कैसे गति दें इस पर कोई विचार? धन्यवाद!

+0

कुछ दिन पहले मेरी प्रेमिका मुझे सुंदर भूखंडों वह आर के [ 'smoothScatter'] (http://rfunction.com/archives/595) समारोह है, जो फ़ायदेमंद एक को जोड़ती है के साथ बनाया गया है पता चला है स्कैटर प्लॉट और घनत्व मानचित्र। मैं तुरंत निराश हो गया कि matplotlib में कोई समकक्ष नहीं था, इसलिए मुझे इस पुराने प्रश्न को इसके बारे में यहां खुशी हुई। – Julien

उत्तर

13

यह करना चाहिए:

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() 

image

+0

अच्छा, यह मेरे अंतिम समाधान के समान विचार है लेकिन कोड की कम पंक्तियों में व्यक्त किया गया है। धन्यवाद! – Singularity

+0

क्या वही काम करने का कोई तरीका है, लेकिन गतिशील साजिश फिर से स्केल के साथ? उदाहरण के लिए जहां मानक विचलन बहुत अलग हैं? – chiffa

+0

'np.histogram2d' भी काम करता है, 'scipy' आयात करने की कोई आवश्यकता नहीं है – Mathias711

2

आपकी समस्या वर्गबद्ध है - npts = 1000 के लिए, आपके पास सरणी आकार 10^6 अंक तक पहुंच गया है, और आप इन सूचियों पर सूची समझ के साथ पुन: प्रयास करते हैं।
अब, यह निश्चित रूप से स्वाद का विषय है, लेकिन मुझे लगता है कि सूची समझ पूरी तरह से कोड उत्पन्न कर सकती है जो कि पालन करना मुश्किल है, और वे कभी-कभी थोड़ा तेज होते हैं ... लेकिन यह मेरा मुद्दा नहीं है।
मेरे बिंदु आप की तरह numpy कार्यों है कि बड़े सरणी के संचालन के लिए है:

np.where, np.choose etc. 

देखें कि आप NumPy साथ सूची comprehensions की है कि कार्यक्षमता प्राप्त कर सकते हैं, और अपने कोड तेजी से चलाना चाहिए।

क्या मैं सही ढंग से आपकी टिप्पणी समझता हूं?

#Find all points that lie in these regions 

क्या आप बहुभुज के अंदर एक बिंदु के लिए परीक्षण कर रहे हैं? यदि हां, तो matplotlib के अंदर point in polygon पर विचार करें।

1

रात को सोने के लिए और ओज़ 123 के सुझावों के माध्यम से पढ़ने के बाद, मैंने इसे समझ लिया। यह चाल गणना करने के लिए है कि कौन सा बिन प्रत्येक एक्स, वाई बिंदु (xi, yi) में आता है, फिर जांच करें कि एच [xi, yi] (वास्तव में, मेरे मामले में एच [यी, xi]) दहलीज के नीचे है। नीचे दिए गए कोड है, और अंक की बड़ी संख्या के लिए बहुत तेजी से चलता है और अधिक स्वच्छ है:

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') 
+0

मैंने आपको अपने सुझाव को लागू करने के लिए एक अपवर्तित दिया :-) हमेशा नमस्ते बिल्टिन के साथ काम करने का प्रयास करें, यह सूची समझों की तुलना में तेज़ है – Oz123

4

रिकॉर्ड के लिए, यहाँ है 2 डी हिस्टोग्राम की बजाय scipy.stats.gaussian_kde का उपयोग करके एक नए प्रयास का नतीजा। कोई उद्देश्य के आधार पर रंग जाल और समोच्च के विभिन्न संयोजनों को कल्पना कर सकता है।

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() 

Smooth scatter plot

संबंधित मुद्दे