如何在python中绘制核密度估计(KDE)和三维数据的交叉点

Question

我有3D数据集(X,Y,Z)。我想执行KDE，绘制数据及其估计。然后，获取零交叉点并使用KDE绘制它。我的尝试如下。我有以下问题:

1. 行X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]和positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()])在这里(kde文档)，他们会对可视化原始数据的真实估计有任何影响吗?我真的不明白为什么我必须使用我的min和max来执行KDE，然后使用ravel() ?

2. 为什么我必须转换f = np.reshape(kernel(positions).T, X.shape)中的数据

4. 我无法绘制KDE估计的原始数据和KDE估计/原始数据的零交叉:

5. 过零应该是矢量吗?在下面的代码中，它是元组

   df = pd.read_csv(file, delimiter = ',')
   Convert series from data-frame into arrays
   X = np.array(df['x']) 
   Y = np.array(df['y']) 
   Z = np.array(df['z'])
   data = np.vstack([X, Y, Z])
   # perform KDE
   kernel = scipy.stats.kde.gaussian_kde(data)
   density = kernel(data)
   fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
   x, y, z = data
   scatter = ax.scatter(x, y, z, c=density)
   xmin = values[0].min()
   xmax = values[0].max()
   ymin = values[1].min()
   ymax = values[1].max()
   zmin = values[2].min()
   zmax = values[2].max()
   X,Y, Z =      np.mgrid[xmin:xmax:100j,ymin:ymax:100j,zmin:zmax:100j]
   positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()])
   
   f = np.reshape(kernel(positions).T, X.shape)
   derivative = np.gradient(f)
   dz, dy, dx = derivative
   xdiff = np.sign(dx)   # along X-axis 
   ydiff = np.sign(dy)   # along Y-axis 
   zdiff = np.sign(dz)   # along Z-axis
   xcross = np.where(xdiff[:-1] != xdiff[1:])
   ycross = np.where([ydiff[:-1] != ydiff[1:]])
   zcross = np.where([zdiff[:-1] != zdiff[1:]])
   Zerocross =  xcross + ycross + zcross
   

Answer 1

行X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]和positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()])，如这里(kde文档)，它们对可视化原始数据的真实估计有任何影响吗?我真的不明白为什么我必须使用我的min和max来执行KDE，然后使用ravel() ?

这两行设置了一个由x、y、z位置组成的网格，KDE将在这些位置被求值。在上面的代码中，它们只被用来估计核密度函数的导数。由于它们目前没有被用于任何与绘图相关的事情，因此它们不会影响可视化。

xmin, xmax等用于确保网格覆盖数据中x, y, z值的全部范围。语法xmin:xmax:100j相当于np.linspace(xmin, xmax, 100)，即np.mgrid在xmin和xmax之间返回100个均匀间隔的点。

np.mgrid返回的X, Y和Z数组的形状都是(100, 100, 100)，而kernel(positions)的positions参数必须是(n_dimensions, n_points)。np.vstack([X.ravel(),Y.ravel(),Z.ravel()])行只是将np.mgrid的输出重新塑造成这种形式。.ravel()将每个(100, 100, 100)数组平展成一个(1000000,)向量，np.vstack将它们在第一维上连接成一个(3, 1000000)点数组。

为什么我必须转换f = np.reshape(kernel(positions).T, X.shape)中的数据

你不:-)。kernel(positions)的输出是一个1D向量，所以对它进行转置不会产生任何影响。

我无法绘制KDE估计的原始数据和KDE估计/零交叉的原始数据:

你试了什么?上面的代码似乎估计了核密度函数梯度的过零点，但没有包含任何代码来绘制它们。你想拍什么样的情节?

过零应该是向量吗?在下面的代码中，它是元组

当你调用np.where(x)，其中x是一个多维数组，你得到一个元组，其中x是非零的索引。由于xdiff[:-1] != xdiff[1:]是一个3D数组，您将返回一个包含三个一维索引数组的元组，每个维度一个索引数组。

您可能不希望在np.where([ydiff[:-1] != ydiff[1:]])中使用额外的方括号，因为在这种情况下，[ydiff[:-1] != ydiff[1:]]将被视为(1, 100, 100, 100)数组而不是(100, 100, 100)，因此您将获得包含4个索引数组而不是3个数组的元组(第一个数组将全为零，因为第一个维度的大小为1)。