2015-07-20 175 views
1

我有3D數據集(X,Y,Z)。我想執行KDE,繪製數據和估計。然後,獲得過零點並用KDE繪製它。我的嘗試如下。我有以下問題:如何在python中繪製3D數據的核心密度估計(KDE)和過零點?

  1. X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()])here(KDE文檔)將他們在可視化的原始數據真實估計任何影響?我真的不明白爲什麼我必須使用我的min和max來執行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 
    
+0

三維數據加密度總共是四維的。以有意義的方式可視化這些數據非常困難。 2-D KDE非常容易操作(如果使用'seaborn',只需一行)。也許考慮通過PCA降維來將3D轉換爲2D而不會丟失太多信息。 –

+1

因爲'scipy.stats.gaussian_kde'確切地計算了KDE,所以永遠不會出現零交叉(除了在無窮遠處)。你想解決什麼問題? –

+2

@JoeKington它看起來像她想估計內核密度函數的導數的零交叉 –

回答

2

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]positions = np.vstack([X.ravel(),Y.ravel(),Z.ravel()])here(KDE文檔),他們將在可視化的原始數據真實估計任何影響?我不明白爲什麼我必須使用我的minmax來執行KDE,然後使用ravel()

這兩行建立了一個x,y,z位置的網格,KDE將被評估。在上面的代碼中,它們僅用於估計內核密度函數的導數。由於它們目前沒有用於與繪圖有關的任何事情,因此它們不會影響可視化。

xmin,xmax等用於確保網格覆蓋數據中x,y,z值的全部範圍。語法xmin:xmax:100j等效於np.linspace(xmin, xmax, 100),即np.mgridxminxmax之間返回100個均勻間隔的點。

通過np.mgrid將每個返回的XYZ陣列具有形狀(100, 100, 100),而positions參數kernel(positions)需要是(n_dimensions, n_points)。行np.vstack([X.ravel(),Y.ravel(),Z.ravel()])只是將np.mgrid的輸出重新整形爲這種形式。 .ravel()將每個(100, 100, 100)數組變爲(1000000,)矢量,並將它們連接到第一維以創建(3, 1000000)點數組。

爲什麼我要轉的數據f = np.reshape(kernel(positions).T, X.shape)

你不:-)。 kernel(positions)的輸出是一維矢量,因此移調它將不起作用。

我失敗過零密謀與KDE估計和KDE估計/原始數據的原始數據:

那你試試?上面的代碼似乎估計了內核密度函數梯度的零交叉,但不包括任何代碼來繪製它們。你想製作什麼樣的情節?

過零應該是矢量嗎?在它下面的元組

當你調用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)。

+0

非常感謝您的詳細解答。但正如你所提到的那樣,在導數中使用了網格,1)這是否會影響爲原始數據的導數得到正確的過零點?我將KDE結果繪製爲散點圖(x,y,z,c =密度)。但是,2)我不知道展示流程的最佳方式;使用KDE繪製原始數據,顯示一階導數和零交點(峯值)。您可以給我關於如何顯示從原始數據發生的變化,直到過零點顯示峯值的指導。 3)如果我想獲得過零點的數量,是否應該只有len(Zerocross)? – Yasmin

+0

另外,4)100j或50j如何影響一階導數?當我爲他們兩人分散情節。這些點位於X和Y軸上的相同位置,但它們的密度值稍有變化。 – Yasmin