Estimation par noyau 2d (scipy + matplotlib)

Exemple d'estimation par noyau 2d utilisant scipy et matplotlib:

from scipy import stats, mgrid, c_, reshape, random, rot90

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

values = np.array([[ 0.04810362,  0.02606397],
       [ 0.21804558,  0.19724609],
       [ 0.02868612,  0.00172462],
       [ 0.40554094,  0.35376897],
       [ 0.43878371,  0.19474262],
       [ 0.37266064,  0.37235034],
       [ 0.05908098,  0.03568844],
       [ 0.45307499,  0.35184965],
       [ 0.11852442,  0.19123775],
       [ 0.14187719,  0.10859522],
       [ 0.49605238,  0.11632819],
       [ 0.04044018,  0.00787204],
       [ 0.19557305,  0.14191927],
       [ 0.11319108,  0.14575793],
       [ 0.20686109,  0.24453397],
       [ 0.12825906,  0.14751036],
       [ 0.06721044,  0.04319887],
       [ 0.12370241,  0.16036154],
       [ 0.22565724,  0.19159937],
       [ 0.37266064,  0.31799707],
       [ 0.1511976 ,  0.18545194],
       [ 0.0774111 ,  0.05279553],
       [ 0.06586416,  0.0569958 ],
       [ 0.06011658,  0.01246175],
       [ 0.03505506,  0.0020306 ],
       [ 0.14648561,  0.10317103],
       [ 0.08134638,  0.05424198],
       [ 0.04820718,  0.03045896],
       [ 0.49315271,  0.3565228 ],
       [ 0.03251784,  0.00161335],
       [ 0.03774762,  0.00734353],
       [ 0.12587717,  0.08194711],
       [ 0.41677719,  0.34954089],
       [ 0.15342413,  0.31607774],
       [ 0.35034347,  0.29660627],
       [ 0.05555994,  0.02828928],
       [ 0.11800662,  0.08386645],
       [ 0.33558616,  0.33916536],
       [ 0.33879653,  0.34063962],
       [ 0.02946282,  0.0016968 ],
       [ 0.35386452,  0.31045884],
       [ 0.03785118,  0.01157162],
       [ 0.18231738,  0.18567447],
       [ 0.31047288,  0.19805276],
       [ 0.10283507,  0.07913766],
       [ 0.06384474,  0.09015295],
       [ 0.30016866,  0.2060917 ],
       [ 0.08113926,  0.07668982],
       [ 0.07078326,  0.08063974],
       [ 0.02863434,  0.00225313],
       [ 0.38897136,  0.22230867],
       [ 0.03510684,  0.00553546],
       [ 0.42770278,  0.36172447],
       [ 0.08766354,  0.06161333],
       [ 0.22576079,  0.25660628],
       [ 0.0973464 ,  0.0717385 ],
       [ 0.59655738,  0.13218354],
       [ 0.17967659,  0.15104306],
       [ 0.23269932,  0.27379683],
       [ 0.06425898,  0.04066757],
       [ 0.13276392,  0.06013906],
       [ 0.32585153,  0.23766331],
       [ 0.04598064,  0.02447843],
       [ 0.05545638,  0.02965228],
       [ 0.0976053 ,  0.06956881],
       [ 0.27966377,  0.28720433],
       [ 0.16000019,  0.13368562],
       [ 0.42356038,  0.35401931],
       [ 0.19324295,  0.12570231],
       [ 0.09227195,  0.0608901 ],
       [ 0.09304865,  0.13465919],
       [ 0.04214892,  0.01329624],
       [ 0.15088691,  0.16662024],
       [ 0.0745632 ,  0.05068148],
       [ 0.14922996,  0.16720438],
       [ 0.2956638 ,  0.33148804],
       [ 0.10262796,  0.07232264],
       [ 0.18583842,  0.2508761 ],
       [ 0.10878978,  0.21972175],
       [ 0.18806495,  0.18892899],
       [ 0.23948249,  0.23866472],
       [ 0.15466686,  0.09268425],
       [ 0.27122363,  0.24216957],
       [ 0.34899718,  0.29974952],
       [ 0.24129479,  0.24970782],
       [ 0.39782572,  0.3470096 ],
       [ 0.12494513,  0.12892902],
       [ 0.0590292 ,  0.0327399 ],
       [ 0.37768331,  0.28645328],
       [ 0.63627261,  0.37452   ],
       [ 0.0862137 ,  0.1038108 ],
       [ 0.2467317 ,  0.31009722],
       [ 0.04810362,  0.01694019],
       [ 0.0323625 ,  0.0020306 ],
       [ 0.10195482,  0.11632819],
       [ 0.14239499,  0.109235  ],
       [ 0.19407143,  0.10489564],
       [ 0.03220716,  0.00319889],
       [ 0.19764425,  0.23963828],
       [ 0.02884146,  0.00837274]])

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

kernel = stats.kde.gaussian_kde(values.T)

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 = reshape(kernel(positions.T).T, X.T.shape)

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

ax.plot([xmin,xmax], [ymin,ymax], 'b--')

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")
plt.show()

Recherches associées