Comment calculer et tracer une estimation par noyau (densité) en 2D avec Python ?

L'estimation de densité par noyau (KDE) est une méthode non paramétrique permettant d'estimer les fonctions de densité de probabilité à partir de données d'échantillon. L'idée de base est de construire une fonction lisse qui peut être utilisée pour approximer la distribution de probabilité sous-jacente, en se basant sur un ensemble donné d'observations en deux dimensions.

Créer des données synthétiques

En utilisant Python, il est assez simple de calculer et de tracer un KDE 2D. La première étape consiste à importer les modules nécessaires, y compris numpy, scipy et matplotlib. Ensuite, préparez vos données pour le calcul - elles doivent être sous forme d'un tableau ou d'une liste de points bidimensionnels.

import numpy as np

mu = 5.0
sigma = 4.0

x = np.random.randn(10000) * sigma + mu

mu = 5.0
sigma = 10.0

y = np.random.randn(10000) * sigma + mu

Calculez une estimation de densité de noyau 2D.

Une fois que vos données sont prêtes, vous pouvez utiliser la fonction scipy.stats.gaussian_kde() pour calculer la KDE. Cette fonction renvoie un objet contenant toutes les informations nécessaires pour tracer le résultat.

from scipy import stats, mgrid, c_

data = np.stack((x,y), axis=1)

data.shape

xmin = data.min()
xmax = data.max()
ymin = data.min()
ymax = data.max()

kernel = stats.gaussian_kde(data.T)

Tracez le résultat en utilisant matplotlib.

Pour tracer un graphique, vous pouvez utiliser matplotlib ou seaborn, en fonction de la complexité du graphique souhaité. Avec matplotlib, il est assez simple de tracer un KDE 2D de base en utilisant les fonctions contourf() ou imshow(). Seaborn offre plus d'options de personnalisation et vous permet de créer des graphiques plus complexes avec la fonction kdeplot().

Avec imshow()

import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure(figsize=(10,4))

ax = fig.add_subplot(111)

X, Y = mgrid[xmin:xmax:100j, ymin:ymax:100j]

positions = c_[X.ravel(), Y.ravel()]

Z = np.reshape(kernel(positions.T).T, X.T.shape)

cax = ax.imshow(Z.T, interpolation='nearest',
            cmap=cm.jet, origin='lower', extent=[xmin, xmax, ymin, ymax])

plt.xlabel('x',fontsize=14)
plt.ylabel('y',fontsize=14)

plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)

ticks_at = [0, Z.max()/6.0, 2.0*Z.max()/6.0, 3.0*Z.max()/6.0, 4.0*Z.max()/6.0, 5.0*Z.max()/6.0, Z.max()]
cbar = fig.colorbar(cax,ticks=ticks_at,format='%1.2g')
cbar.set_label('Kernel density function')

plt.savefig("kernel2D.png", dpi=100, bbox_inches='tight')
plt.show()

Avec contourf()

fig = plt.figure(figsize=(10,4))

ax = fig.gca()

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

cfset = ax.contourf(X, Y, Z, cmap='Blues')

cset = ax.contour(X, Y, Z, colors='k')

ax.clabel(cset, inline=1, fontsize=10)
ax.set_xlabel('X')
ax.set_ylabel('Y')

plt.savefig("kernel2D_contourf.png", dpi=100, bbox_inches='tight')

plt.show()

Références