Python中的多元核密度估計

Question

我試圖使用SciPy的gaussian_kde函數來估計多元數據的密度。在我的代碼中，我採樣了一個3D多變量法線，並且適合內核密度，但我不知道如何評估我的適合度。

import numpy as np 
from scipy import stats 

mu = np.array([1, 10, 20]) 
sigma = np.matrix([[4, 10, 0], [10, 25, 0], [0, 0, 100]]) 
data = np.random.multivariate_normal(mu, sigma, 1000) 
values = data.T 
kernel = stats.gaussian_kde(values) 

我看到this，但不知道如何將其擴展到3D。

也不知道我怎麼開始評估擬合密度？我如何可視化這個？

Answer 1

有幾種方法可以在3D中可視化結果。

最簡單的方法是在您用來生成高斯KDE的點上評估高斯KDE，然後通過密度估計對點進行着色。

例如：

import numpy as np 
from scipy import stats 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 

mu=np.array([1,10,20]) 
sigma=np.matrix([[4,10,0],[10,25,0],[0,0,100]]) 
data=np.random.multivariate_normal(mu,sigma,1000) 
values = data.T 

kde = stats.gaussian_kde(values) 
density = kde(values) 

fig, ax = plt.subplots(subplot_kw=dict(projection='3d')) 
x, y, z = values 
ax.scatter(x, y, z, c=density) 
plt.show() 

如果你有一個更復雜的（即不是所有的趴在一個平面上）分佈的，那麼你可能要評估的規則的三維網格的KDE和可視化體積的等值面（3D輪廓）。這是最容易使用Mayavi的爲visualiztion：

import numpy as np 
from scipy import stats 
from mayavi import mlab 

mu=np.array([1,10,20]) 
# Let's change this so that the points won't all lie in a plane... 
sigma=np.matrix([[20,10,10], 
       [10,25,1], 
       [10,1,50]]) 

data=np.random.multivariate_normal(mu,sigma,1000) 
values = data.T 

kde = stats.gaussian_kde(values) 

# Create a regular 3D grid with 50 points in each dimension 
xmin, ymin, zmin = data.min(axis=0) 
xmax, ymax, zmax = data.max(axis=0) 
xi, yi, zi = np.mgrid[xmin:xmax:50j, ymin:ymax:50j, zmin:zmax:50j] 

# Evaluate the KDE on a regular grid... 
coords = np.vstack([item.ravel() for item in [xi, yi, zi]]) 
density = kde(coords).reshape(xi.shape) 

# Visualize the density estimate as isosurfaces 
mlab.contour3d(xi, yi, zi, density, opacity=0.5) 
mlab.axes() 
mlab.show()