Apprentissage de mélanges de lois Gaussiennes

by Joseph Razik, last modified on 2020-11-10
gaussiennes_1D

Apprentissage de GMM par l'algorithme EM

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

On définit une fonction pour calculer la probabilité qu'une valeur appartienne à une gaussienne à l'aide de sa fonction de densité ($\mu$ sa moyenne et $\sigma$ son écart-type)

$$\mathcal{N}(x, \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}\exp^{-\frac{1}{2\sigma^2}(x-\mu)^2}$$
In [2]:
def gauss(x, mu, sigma):
    """ fonction calculant la probabilité d'appartenance d'un point à une gaussienne (sorte de distance)
        sigma est l'écart-type, c'est-à-dire sqrt(variance). """
    return 1/(sqrt(2*pi)*sigma) * exp(-0.5*(x - mu)**2/sigma**2)

À quoi ressemble une gaussienne centrée-réduite typique $\mathcal{N}(0, 1)$ (de moyenne nulle et d'écart-type 1) ?

In [3]:
mu = 0.0
sigma = 1.0
G = [mu, sigma]
figure(figsize=(10,6))
plot(linspace(-5, 5, 1000), [gauss(x, mu, sigma) for x in linspace(-5, 5, 1000)])
Out[3]:
[<matplotlib.lines.Line2D at 0x7f6275558b20>]

La séquence de données à modéliser

In [4]:
data = [11, 9, 10, 18, 19, 18, 14, 15, 14]

L'ordre des données n'a pas d'importance pour une modélisation par gaussiennes

 Modélisation par une seule gaussienne

In [5]:
# on calcule simplement la moyenne et l'écart-type des données
mu_init = mean(data)
sigma_init = sqrt(var(data))  # écart-type = sqrt(variance)
print("moyenne: ", mu_init)
print("écart-type: ", sigma_init)

moyenne:  14.222222222222221
écart-type:  3.456966485800899
In [6]:
# on représente les données et le modèle gaussien associé
figure(figsize=(10,6))
plot(linspace(0, 30, 1000), [gauss(x, mu_init, sigma_init) for x in linspace(0, 30, 1000)])
plot(data, [0]*len(data), 'ro')
Out[6]:
[<matplotlib.lines.Line2D at 0x7f627344f9d0>]

Une estimation relative que le modèle correspond bien aux données peut être obtenue en calculant la log-vraisemblance. Toutefois, cette estimation en soit ne donne aucun indice sur la qualité absolue du modèle. Elle sert pour estimer l'amélioration du modèle d'une itération à l'autre.

La log-vraisemblance d'un système données-modéle est définie ainsi:

$$LL = \sum_{i=1}^N \log\big(\sum_{k=1}^M c_k\mathcal{N}(x_i, \mu_k, \sigma_k)\big)$$
In [7]:
def calc_LL(data, mu, sigma, poids):
    """ fonction qui retourne le Log-Likelihood (vraisemblance) d'un jeu de données par rapport à
        un modèle GMM donné. """
    ll = 0.0
    for i in range(len(data)):
        p = 0.0
        for k in range(len(mu)):
            p += poids[k] * gauss(data[i], mu[k], sigma[k])
        ll += log(p)
    return ll
In [8]:
calc_LL(data, [mu_init], [sigma_init], [1])
Out[8]:
-23.933969995560677

Modélisation par plusieurs gaussiennes

In [9]:
# on définit 3 gaussiennes avec des moyennes déterminées de manière empirique
# ou bien à partir des données (mu = mu_init ±sigma)
mu = [mu_init - sigma_init, mu_init, mu_init + sigma_init]
mu = [9, 10, 11]
# on garde les mêmes sigmas
sigma = [sigma_init] * 3

# on définit des pondérations identiques pour les 3 gaussiennes
poids = [1/3] * 3
In [10]:
mu, sigma, poids
Out[10]:
([9, 10, 11],
 [3.456966485800899, 3.456966485800899, 3.456966485800899],
 [0.3333333333333333, 0.3333333333333333, 0.3333333333333333])
In [11]:
# on représente les données et les 3 modèles gaussiens associés
figure(figsize=(10,6))
plot(linspace(0, 30, 1000), [poids[0]*gauss(x, mu[0], sigma[0]) for x in linspace(0, 30, 1000)])
plot(linspace(0, 30, 1000), [poids[1]*gauss(x, mu[1], sigma[1]) for x in linspace(0, 30, 1000)])
plot(linspace(0, 30, 1000), [poids[2]*gauss(x, mu[2], sigma[2]) for x in linspace(0, 30, 1000)])
plot(data, [0]*len(data), 'ro')
Out[11]:
[<matplotlib.lines.Line2D at 0x7f62733ce280>]

Il y a de gros recouvrements entre les différentes gaussiennes. il faut améliorer les modèles en appliquant l'algorithme EM

Quelle est sa log-vraisemblance ?

In [12]:
LL = []
LL.append(calc_LL(data, mu, sigma, poids))
print(LL[-1])

-30.29878671556341
In [13]:
def expectation(mu, sigma, poids):
    """ partie Expectation: on calcule les probabilités qu'une donnée appartienne à chacune des données. """
    gamma_i_k = empty((len(data), len(mu)))
    for i in range(len(data)):
        for k in range(len(mu)):
            gamma_i_k[i, k] = poids[k]*gauss(data[i], mu[k], sigma[k])
        # on normalise
        s = sum(gamma_i_k[i, :])
        gamma_i_k[i, :] /= s 
    return gamma_i_k
In [14]:
gamma_i_k = expectation(mu, sigma, poids)
In [15]:
gamma_i_k
Out[15]:
array([[0.3015765 , 0.34190751, 0.35651599],
       [0.35651599, 0.34190751, 0.3015765 ],
       [0.32865262, 0.34269477, 0.32865262],
       [0.14596018, 0.29725613, 0.55678369],
       [0.12946898, 0.28668364, 0.58384739],
       [0.14596018, 0.29725613, 0.55678369],
       [0.22673912, 0.33041509, 0.44284579],
       [0.2043591 , 0.32379351, 0.47184739],
       [0.22673912, 0.33041509, 0.44284579]])
In [16]:
def maximisation(mu, sigma, poids, gamma_i_k):
    """ partie Maximisation: on calcule les nouveaux paramètres des modèles. """
    for k in range(len(mu)):
        # les nouvelles moyennes
        mu[k] = sum([gamma_i_k[i, k]*data[i] for i in range(len(data))])
        mu[k] /= sum(gamma_i_k[:, k])
        # les nouveaux écart-types
        sigma[k] = sum([gamma_i_k[i, k]*(data[i] - mu[k])**2 for i in range(len(data))])
        sigma[k] /= sum(gamma_i_k[:, k])
        sigma[k] = sqrt(sigma[k])  # passage de la variance à l'écart-type
        # les nouveaux poids
        poids[k] = sum(gamma_i_k[:, k]) / len(data)
    return mu, sigma, poids
In [17]:
mu, sigma, poids = maximisation(mu, sigma, poids, gamma_i_k)
In [18]:
mu, sigma, poids
Out[18]:
([13.040386450592509, 14.010103151308401, 14.978131482907527],
 [3.343455828873287, 3.434170114597121, 3.3376721407093304],
 [0.22955242092500402, 0.32136993022516197, 0.449077648849834])
In [19]:
# on représente les données et les 3 modèles gaussiens associés
figure(figsize=(10,6))
plot(linspace(0, 30, 1000), [poids[0]*gauss(x, mu[0], sigma[0]) for x in linspace(0, 30, 1000)])
plot(linspace(0, 30, 1000), [poids[1]*gauss(x, mu[1], sigma[1]) for x in linspace(0, 30, 1000)])
plot(linspace(0, 30, 1000), [poids[2]*gauss(x, mu[2], sigma[2]) for x in linspace(0, 30, 1000)])
plot(data, [0]*len(data), 'ro')
Out[19]:
[<matplotlib.lines.Line2D at 0x7f627333e9d0>]
In [20]:
# calcul de la nouvelle vraisemblance
LL.append(calc_LL(data, mu, sigma, poids))
print(LL[-1])

-23.92936581332484
In [21]:
# on itère le processus 4 fois
figure(figsize=(12,12))

for i in range(4):
    gamma_i_k = expectation(mu, sigma, poids)
    mu, sigma, poids = maximisation(mu, sigma, poids, gamma_i_k)
    subplot(2, 2, i+1)
    plot(linspace(0, 30, 1000), [poids[0]*gauss(x, mu[0], sigma[0]) for x in linspace(0, 30, 1000)])
    plot(linspace(0, 30, 1000), [poids[1]*gauss(x, mu[1], sigma[1]) for x in linspace(0, 30, 1000)])
    plot(linspace(0, 30, 1000), [poids[2]*gauss(x, mu[2], sigma[2]) for x in linspace(0, 30, 1000)])
    plot(data, [0]*len(data), 'ro')
    LL.append(calc_LL(data, mu, sigma, poids))
    title("itération " + str(i) + ", LL= " + str(LL[-1]))
In [22]:
mu, sigma, poids
Out[22]:
([12.865548096200117, 13.97635602994447, 15.095563295109924],
 [3.2786545207900875, 3.4609698651934635, 3.2857341930157875],
 [0.2306689387650227, 0.3207143526469538, 0.44861670858802355])
In [23]:
# on itère le processus 30 fois
figure(figsize=(15,30))

for i in range(30):
    gamma_i_k = expectation(mu, sigma, poids)
    mu, sigma, poids = maximisation(mu, sigma, poids, gamma_i_k)
    subplot(10, 3, i+1)
    plot(linspace(0, 30, 1000), [poids[0]*gauss(x, mu[0], sigma[0]) for x in linspace(0, 30, 1000)])
    plot(linspace(0, 30, 1000), [poids[1]*gauss(x, mu[1], sigma[1]) for x in linspace(0, 30, 1000)])
    plot(linspace(0, 30, 1000), [poids[2]*gauss(x, mu[2], sigma[2]) for x in linspace(0, 30, 1000)])
    plot(data, [0]*len(data), 'ro')
    LL.append(calc_LL(data, mu, sigma, poids))
    title("itération " + str(i) + ", LL= " + str(LL[-1]))

Voyons l'évolution de la log-vraisemblance sur ces itérations successives du modèle

In [24]:
plot(LL)
Out[24]:
[<matplotlib.lines.Line2D at 0x7f62711b7a00>]

En utilisant une bibliothèque dédiée

In [25]:
from sklearn import mixture
In [26]:
g = mixture.GaussianMixture(1, covariance_type='spherical')
In [27]:
g.fit(array(data).reshape(-1,1))
Out[27]:
GaussianMixture(covariance_type='spherical')
In [28]:
g.means_
Out[28]:
array([[14.22222222]])
In [29]:
sqrt(g.covariances_)
Out[29]:
array([3.45696663])
In [30]:
g.weights_
Out[30]:
array([1.])
In [31]:
Z = [exp(g.score(x.reshape(-1,1))) for x in linspace(0, 30, 1000)]
In [32]:
plot(linspace(0, 30, 1000), Z)
plot(data, [0]*len(data), 'ro')
Out[32]:
[<matplotlib.lines.Line2D at 0x7f6272fedeb0>]
In [33]:
g3 = None
g3 = mixture.GaussianMixture(n_components=3,
                             covariance_type='spherical',
                             # n_init = 1,
                             # means_init = array([mu_init - sigma_init, mu_init, mu_init + sigma_init]).reshape(-1,1),
                             # means_init = array([10, 11, 9]).reshape(-1,1),
                             # precisions_init = array([sigma_init]*3),
                             # weights_init=[1/3, 1/3, 1/3],
                             # max_iter=1,
                             verbose=2,
                             verbose_interval=1,
                             # init_params = 'random',
                             # warm_start=True,
                            )
In [34]:
g3.fit(array(data).reshape(-1,1), 1)

Initialization 0
  Iteration 1	 time lapse 0.00291s	 ll change inf
  Iteration 2	 time lapse 0.00035s	 ll change 0.00000
Initialization converged: True	 time lapse 0.00330s	 ll -1.94861
Out[34]:
GaussianMixture(covariance_type='spherical', n_components=3, verbose=2,
                verbose_interval=1)
In [35]:
g3.means_
Out[35]:
array([[10.00001218],
       [14.33333434],
       [18.33333333]])
In [36]:
sqrt(g3.covariances_)
Out[36]:
array([0.81652579, 0.47140594, 0.47140558])
In [37]:
g3.weights_
Out[37]:
array([0.33333435, 0.33333232, 0.33333333])
In [38]:
Z = [exp(g3.score_samples(x.reshape(-1,1))) for x in linspace(0, 30, 1000)]
In [39]:
plot(linspace(0, 30, 1000), Z)
plot(data, [0]*len(data), 'ro')
Out[39]:
[<matplotlib.lines.Line2D at 0x7f6273032cd0>]
In [40]:
g3.n_iter_
Out[40]:
2
In [41]:
g3 = None
g3 = mixture.GaussianMixture(n_components=3,
                             covariance_type='spherical',
                             n_init = 1,
                             # means_init = array([mu_init - sigma_init, mu_init, mu_init + sigma_init]).reshape(-1,1),
                             means_init = array([10, 11, 9]).reshape(-1,1),
                             precisions_init = array([sigma_init]*3),
                             weights_init=[1/3, 1/3, 1/3],
                             # max_iter=1,
                             verbose=2,
                             verbose_interval=1,
                             init_params = 'random',
                             warm_start=True,
                            )
In [42]:
g3.fit(array(data).reshape(-1,1), 1)

Initialization 0
  Iteration 1	 time lapse 0.00065s	 ll change inf
  Iteration 2	 time lapse 0.00063s	 ll change 128.77986
  Iteration 3	 time lapse 0.00024s	 ll change 0.90619
  Iteration 4	 time lapse 0.00039s	 ll change 0.00009
Initialization converged: True	 time lapse 0.00195s	 ll -1.21883
Out[42]:
GaussianMixture(covariance_type='spherical', init_params='random',
                means_init=array([[10.],
       [11.],
       [ 9.]]),
                n_components=3,
                precisions_init=array([3.45696649, 3.45696649, 3.45696649]),
                verbose=2, verbose_interval=1, warm_start=True,
                weights_init=array([0.33333333, 0.33333333, 0.33333333]))
In [43]:
g3.means_
Out[43]:
array([[10.        ],
       [15.57107529],
       [ 9.        ]])