這篇文章將向您展示如何：

使用高斯核估計2個分布的PDF(PDF即Probability Density Function，概率密度函數)
使用Matplotlib表示PDF，其中密度圖周圍帶有標記的輪廓線
如何提取輪廓線
如何在3D(3維空間)中繪製上述高斯核
如何使用2D直方圖繪製相同的PDF

首先，生成由3個blob組成的輸入數據集：(關於blob：scikit中的make_blobs方法常被用來生成聚類算法的測試數據，直觀地說，make_blobs會根據用戶指定的特征數量、中心點數量、範圍等來生成幾類數據，這些數據可用於測試聚類算法的效果。）

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
from sklearn.datasets.samples_generator import make_blobsn_components = 3
X, truth = make_blobs(n_samples=300, centers=n_components, 
                      cluster_std = [2, 1.5, 1], 
                      random_state=42)plt.scatter(X[:, 0], X[:, 1], s=50, c = truth)
plt.title(f"Example of a mixture of {n_components} distributions")
plt.xlabel("x")
plt.ylabel("y");

為了擬合高斯核，我們指定了一個網格，該網格將在每個軸上使用100點插值(例如mgrid(xmin：xmax：100j))：

# Extract x and y
x = X[:, 0]
y = X[:, 1]# Define the borders
deltaX = (max(x) - min(x))/10
deltaY = (max(y) - min(y))/10xmin = min(x) - deltaX
xmax = max(x) + deltaXymin = min(y) - deltaY
ymax = max(y) + deltaYprint(xmin, xmax, ymin, ymax)# Create meshgrid
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]

我們將使用scipy的gaussian_kde方法擬合高斯核：

positions = np.vstack([xx.ravel(), yy.ravel()])
values = np.vstack([x, y])
kernel = st.gaussian_kde(values)
f = np.reshape(kernel(positions).T, xx.shape)

用帶注釋的輪廓繪製核

fig = plt.figure(figsize=(8,8))
ax = fig.gca()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
cfset = ax.contourf(xx, yy, f, cmap='coolwarm')
ax.imshow(np.rot90(f), cmap='coolwarm', extent=[xmin, xmax, ymin, ymax])
cset = ax.contour(xx, yy, f, colors='k')
ax.clabel(cset, inline=1, fontsize=10)
ax.set_xlabel('X')
ax.set_ylabel('Y')
plt.title('2D Gaussian Kernel density estimation')

這個matplotlib對象稱為QuadContour set(設置在代碼中)。我們可以通過迭代訪問allsegs對象以編程方式訪問等高線。可從labelTexts訪問計算出的標簽。

plt.figure(figsize=(8,8))for j in range(len(cset.allsegs)):
    for ii, seg in enumerate(cset.allsegs[j]):
        plt.plot(seg[:,0], seg[:,1], '.-', label=f'Cluster{j}, level{ii}')plt.legend()

3D KDE圖

我們將使用matplotlib的axis3d(來自mplot3d)。可以將密度繪製為曲麵：

fig = plt.figure(figsize=(13, 7))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(xx, yy, f, rstride=1, cstride=1, cmap='coolwarm', edgecolor='none')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('PDF')
ax.set_title('Surface plot of Gaussian 2D KDE')
fig.colorbar(surf, shrink=0.5, aspect=5) # add color bar indicating the PDF
ax.view_init(60, 35)

或作為線框：

fig = plt.figure(figsize=(13, 7))
ax = plt.axes(projection='3d')
w = ax.plot_wireframe(xx, yy, f)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('PDF')
ax.set_title('Wireframe plot of Gaussian 2D KDE');

使用2D直方圖表示

呈現相同信息的另一種方法是使用2D直方圖。設定參數normed = False將返回實際頻率，設為True將返回PDF(概率密度函數)。

h =plt.hist2d(x, y)
plt.colorbar(h[3])

本文的完整的代碼可在Github找到。

python二維密度圖的例程

用帶注釋的輪廓繪製核

3D KDE圖

參考資料