2

我正在使用 scipy.stats.kde.gaussian_kde() 进行 kde 分析，处理大量点需要时间（对于 250x250 网格的 100000 个点，需要 5 分钟）。

作为 gaussian_kde 的更快替代方法，我在这里找到了由Joe Kington编写的 fast_kde 函数。（加权 kde 也是选择 fast_kde 的一个因素）

而是绘制结果，我将其提取为格式 (xmin,xmax,ymin,ymax,value) 以供以后使用。我正在使用这种技术通过使用pcolormesh以原始形式提取结果。

这是问题陈述： 由网格（500,500）的fast_kde函数产生的结果不能由pcolormesh绘制，并且原始形式的输出也反映了相同的无效结果，但是imshow方法完美地绘制了这个结果。

生成一些随机的二维数据： 

from scipy import stats
def measure(n):
    "Measurement model, return two coupled measurements."
    m1 = np.random.normal(size=n)
    m2 = np.random.normal(scale=0.5, size=n)
    return m1+m2, m1-m2
m1, m2 = measure(2000)
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()

对数据执行核密度估计： 

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values)
Z = np.reshape(kernel(positions).T, X.shape)

将结果保存到文件：(x,y,value)

fid = open('output.csv','w')
Z1 = (kernel(positions).T, X.shape)
Z = kernel(positions).T
#for currentIndex,elem in enumerate(positions):
for currentIndex,elem in enumerate(Z):
  #if Z1[currentIneex]>0:
  s1 = '%f %f %f\n'%(positions[0][currentIndex], positions[1][currentIndex], Z[currentIndex] )
  fid.write(s1)
fid.close()

打印结果：(minx,maxx,miny,maxy,value) 

mshgrd = ax.pcolormesh(X,Y,Z) 
pths = mshgrd.get_paths() 
arr = mshgrd.get_array() 
for currentIndex,elem in enumerate(pths): 
 if arr[currentIndex]>0: bbox = elem.get_extents() 
 s2 = ",".join([str(i) for i in bbox.extents])+","+ str(arr[currentIndex]) +'\n' 
 print s2

绘制结果： 

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
          extent=[xmin, xmax, ymin, ymax])
ax.plot(m1, m2, 'k.', markersize=2)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
plt.show()

用于 fast_kde 的代码（问题区域） 

kernel = fast_kde(m1,m2,(500,500))
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
mshgrd = ax.pcolormesh(X,Y,kernel)
plt.show()

请建议我如何在此处添加图片（在哪里上传？）

4

0 回答 0