python - 如何在 python 中绘制 3D 数据的核密度估计 (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()])为 here (kde 文档)它们对可视化原始数据的真实估计有什么影响吗？我真的不明白为什么我必须使用最小值和最大值来执行 KDE，然后使用 ravel()？

2. 为什么我必须转置 f = np.reshape(kernel(positions).T, X.shape)

中的数据

3. 代码正确吗？

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

line X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] and positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()]) as here (kde documentation) will they have any effects in visualising the real estimation for the original data?. I don't really understand why I have to use my min and max to perform KDE and then use 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 的输出 reshape 为:形式。 .ravel()将每个 (100, 100, 100) 数组展平为 (1000000,) 向量，并且 np.vstack将它们连接到第一个维度以形成一个 (3, 1000000) 点数组。

why I have to transpose the data in f = np.reshape(kernel(positions).T, X.shape)

你不:-)。 kernel(positions) 的输出是一个一维向量，因此转置它不会产生任何效果。

I failed to plot the original data with KDE estimation and KDE estimation/ original data with zero crossing:

你尝试了什么？上面的代码似乎估计了核密度函数梯度的零交叉，但不包含任何绘制它们的代码。你想创作什么样的情节？

Should zero crossings be vector ?. In the code below it's tuple

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

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