def gauss(mu, sigma, x):
    """ fonction calculant la probabilité d'appartenance d'un point à une gaussienne (sorte de distance)
        attention, en dimension 2 cette fois
    """
    return 1/(sqrt(det(sigma)*2*pi**2))*exp(-0.5*(dot(dot((x-mu).T, inv(sigma)), (x-mu))))

# a quoi ressemble une gaussienne type
mu = array([0.0, 0.0])
sigma = array([[1.0, 0.0], [0.0, 1.0]])
G = [mu, sigma]
figure(figsize=(10,10))
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
Z = empty((X.shape[0], Y.shape[0]))
i = 0
for x in X:
    j = 0
    for y in Y:
        Z[i,j] = gauss(mu, sigma, array([x,y]))
        j+=1
    i+=1

X, Y = np.meshgrid(X, Y)
from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
contour(X, Y, Z) #, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)

#plot(linspace(-5, 5, 1000), [gauss(mu, sigma, t) for t in linspace(-5, 5, 1000)])

<matplotlib.contour.QuadContourSet instance at 0x304e050>

on génére dix valeurs centrées sur (0,0) et dix valeurs centrée sur (2,2) pour obtenir un modèle bimodale

V1 = randn(20,2)
V2 = randn(20,2)*2 + 2
print V1
print V2

[[-2.25739355 -0.23465985]
 [-0.06982311 -0.28916415]
 [ 0.6095009  -0.461931  ]
 [-2.27901933 -0.98775963]
 [ 0.10851471 -0.81672758]
 [ 1.25759542 -0.37583672]
 [ 0.18718211  1.51535751]
 [ 1.86532615  0.40349985]
 [-0.85207867  0.07319375]
 [ 0.42198088  0.04639263]
 [-0.17233774 -1.52433959]
 [ 0.41514289  0.38747269]
 [-0.85944505 -0.04495001]
 [-0.92029531  0.33579699]
 [-0.92511594 -1.56259708]
 [-0.59855765 -0.87654824]
 [ 1.92799762  0.20049376]
 [ 0.81682016 -0.97249486]
 [-1.36019736 -0.05955116]
 [-2.24242556  1.77607575]]
[[ 2.38637607 -0.34017295]
 [ 1.50268985 -0.18340479]
 [-1.12199187 -0.47052334]
 [ 3.44627906  3.35805622]
 [ 4.2159879   1.56895895]
 [ 2.874098    2.66313772]
 [ 2.26959182  3.38547026]
 [ 0.37737804 -1.14557744]
 [ 6.03421103  0.21785147]
 [-0.44180827  0.63267777]
 [-0.51219714  0.81684696]
 [ 2.64158735  3.23749432]
 [ 6.25407438  1.99205722]
 [ 1.9037868   3.87616099]
 [ 2.02412649  0.95644102]
 [-0.08818165  0.91164496]
 [ 0.15333521 -2.35243636]
 [ 2.78834746  5.43990552]
 [ 1.35539505  6.24108236]
 [ 5.38755681  2.09863897]]

On représente toutes les données par une seule gaussienne

mu = mean(V1.tolist()+V2.tolist())
print mu

0.849500585459

#help(cov)
V3 = array(V1.tolist()+V2.tolist()).T
sigma = cov(V3)
print sigma

[[ 4.50278168  1.80473736]
 [ 1.80473736  3.53013358]]

figure(figsize=(10,10))
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
Z = empty((X.shape[0], Y.shape[0]))
i = 0
for x in X:
    j = 0
    for y in Y:
        Z[i,j] = gauss(mu, sigma, array([x,y]))
        j+=1
    i+=1

X, Y = np.meshgrid(X, Y)
from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
contour(X,Y,Z)
#ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.scatter(V3[0,:], V3[1,:], 0, c='r')

<mpl_toolkits.mplot3d.art3d.Patch3DCollection at 0x3203090>

Maintenant on va essayer d'apprendre un modèle à 2 gaussiennes. Pour cela on part de la mono-gaussienne précédente et on décale leur centroide de pas grand chose

mu_1 = mu - array([sigma[0,0], sigma[1,1]])
sigma_1 = sigma
mu_2 = mu + array([sigma[0,0], sigma[1,1]])
sigma_2 = sigma
G1 = [mu_1, sigma_1]
G2 = [mu_2, sigma_2]

# on affiche cette modélisation initiale
figure(figsize=(10,10))
X = arange(-10, 10, 0.25)
Y = arange(-10, 10, 0.25)
Z1 = empty((X.shape[0], Y.shape[0]))
Z2 = empty((X.shape[0], Y.shape[0]))
i = 0
for x in X:
    j = 0
    for y in Y:
        Z1[i,j] = gauss(mu_1, sigma_1, array([x,y]))
        Z2[i,j] = gauss(mu_2, sigma_2, array([x,y]))
        j+=1
    i+=1

X, Y = np.meshgrid(X, Y)
from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
contour(X,Y,Z1)
contour(X,Y,Z2)
#ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.scatter(V3[0,:], V3[1,:], 0, c='r')
ylim((-10,10))
xlim((-10,10))

(-10, 10)

pour chaque valeur échantillon on calcul la distance avec chacune des gaussiennes

d_1_2 = {(x[0],x[1]):(gauss(mu_1, sigma_1, x), gauss(mu_2, sigma_2, x)) for x in V3.T.tolist()}

on affecte les points à la gaussienne la plus proche.

X_1 = [x for x in V3.T.tolist() if d_1_2[(x[0], x[1])][0] < d_1_2[(x[0], x[1])][1]]
X_2 = [x for x in V3.T.tolist() if d_1_2[(x[0], x[1])][0] >= d_1_2[(x[0], x[1])][1]]

print str(len(X_1)) + " " + str(len(X_2))

15 25

On recalcul les moyennes et ecart-types

mu_1 = mean(X_1)
sigma_1 = cov(array(X_1).T)
mu_2 = mean(X_2)
sigma_2 = cov(array(X_2).T)

on affiche les nouvelles gaussiennes

figure(figsize=(10,10))
X = arange(-10, 10, 0.25)
Y = arange(-10, 10, 0.25)
Z1 = empty((X.shape[0], Y.shape[0]))
Z2 = empty((X.shape[0], Y.shape[0]))
i = 0
for x in X:
    j = 0
    for y in Y:
        Z1[i,j] = gauss(mu_1, sigma_1, array([x,y]))
        Z2[i,j] = gauss(mu_2, sigma_2, array([x,y]))
        j+=1
    i+=1

X, Y = np.meshgrid(X, Y)
from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
contour(X,Y,Z1)
contour(X,Y,Z2)
#ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.scatter(V3[0,:], V3[1,:], 0, c='r')
ylim((-10,10))
xlim((-10,10))

(-10, 10)

Et on recommence ... au moins 4 fois pour voir

# l'initialisation est un peu différente de la façon précédente, mais pour que la convergence soit moins rapide.
mu_1 = mu - 0.0003*array([sigma[1,1], sigma[0,0]])
sigma_1 = sigma
mu_2 = mu + 0.0003*array([sigma[1,1], sigma[0,0]])
sigma_2 = sigma
G1 = [mu_1, sigma_1]
G2 = [mu_2, sigma_2]

for k in xrange(4):
    figure(figsize=(10,10))

    d_1_2 = {(x[0],x[1]):(gauss(mu_1, sigma_1, x), gauss(mu_2, sigma_2, x)) for x in V3.T.tolist()}
    X_1 = [x for x in V3.T.tolist() if d_1_2[(x[0], x[1])][0] < d_1_2[(x[0], x[1])][1]]
    X_2 = [x for x in V3.T.tolist() if d_1_2[(x[0], x[1])][0] >= d_1_2[(x[0], x[1])][1]]
    mu_1 = mean(X_1)
    sigma_1 = cov(array(X_1).T)
    mu_2 = mean(X_2)
    sigma_2 = cov(array(X_2).T)

    X = arange(-10, 10, 0.25)
    Y = arange(-10, 10, 0.25)
    Z1 = empty((X.shape[0], Y.shape[0]))
    Z2 = empty((X.shape[0], Y.shape[0]))
    i = 0
    for x in X:
        j = 0
        for y in Y:
            Z1[i,j] = gauss(mu_1, sigma_1, array([x,y]))
            Z2[i,j] = gauss(mu_2, sigma_2, array([x,y]))
            j+=1
        i+=1

    X, Y = np.meshgrid(X, Y)
    from mpl_toolkits.mplot3d import Axes3D
    ax = gca(projection='3d')
    contour(X,Y,Z1)
    contour(X,Y,Z2)
    #ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
    ax.scatter(V3[0,:], V3[1,:], 0, c='r')
    ylim((-10,10))
    xlim((-10,10))

Joseph Razik - UTLN/LIS

Gaussiennes - Apprentissage naïf en 2D