TP2 Expectation-Maximization¶
Le but de ce TP va être d'étudier le temps d'attente entre deux éruptions et la durée de l'éruption du geyser Old Faithful au Yellowstone National Park, Wyoming, USA.
Les données sont les suivantes
import numpy as np
import pandas as pd
import scipy.stats as sps
import matplotlib.pyplot as plt
df = pd.read_csv("Oldfaithful.csv")
df.head(10)
| eruptions | waiting | |
|---|---|---|
| 0 | 3.600 | 79 |
| 1 | 1.800 | 54 |
| 2 | 3.333 | 74 |
| 3 | 2.283 | 62 |
| 4 | 4.533 | 85 |
| 5 | 2.883 | 55 |
| 6 | 4.700 | 88 |
| 7 | 3.600 | 85 |
| 8 | 1.950 | 51 |
| 9 | 4.350 | 85 |
df.plot.scatter(x="eruptions", y="waiting", alpha=0.2)
<AxesSubplot:xlabel='eruptions', ylabel='waiting'>
Si on se concentre sur les temps d'éruption on a
plt.hist(df.eruptions, bins=40, alpha=0.3, density=True);
plt.xlabel("temps d'éruption");
On voit sur ces données deux groupes distincts. Chacun de ces groupes va être modélisés par une loi normale, ce qui nous amène à un modèle de mélange de Gaussienne pour modéliser nos données
$$ f(t) = \pi f_{\mu_1,\sigma_1}(t) + (1-\pi) f_{\mu_2,\sigma_2}(t) $$
où $f_{\mu,\sigma}$ représente la densité de la loi $\mathcal{N}(\mu,\sigma^2)$. On choisira $\mu_1 <\mu_2$ et $\pi \in]0,1[$ représentera la proportion des éruptions de courtes durées.
On va donc chercher à estimer:
- les durées moyennes des éruptions sur les deux groupes $\mu_1$ et $\mu_2$
- la variance de ces durées sur les deux groupes $\sigma_1$ et $\sigma_2$
- les proportions de chaque type d'eruption $\pi$
Pour faire cela, on va utiliser un algorithme appelé EM (pour expectation-maximization). Il s'agit d'un algorithme itératif d'optimisation (avec des applications diverses) qui va chercher les paramètres $\pi$, $\mu_1$, $\mu_1$ et $\sigma_1$ et $\sigma_2$ les mieux adaptés aux données. Cet algorithme va nous fournir les estimateurs du maximum de vraisemblance de ces paramètres.
Une description de la méthode est donnée sur wikipedia: https://fr.wikipedia.org/wiki/Algorithme_esp%C3%A9rance-maximisation
Principe de l'algorithme:¶
Dans ce cadre l'algorithme est le suivant pour passer de l'étape $N$ à $N+1$:
-(E-step) Pour chacune des données $x_i$ on calcule sa probabilité d'appartenir à l'un ou l'autre des groupes
$$ r_{1,i} = \frac{\pi^N f_{\mu^N_1, \sigma^N_1}(x_i)}{\pi^N f_{\mu^N_1, \sigma^N_1}(x_i) + (1-\pi^N) f_{\mu^N_2, \sigma^N_2}(x_i)}\qquad\text{et}\qquad r_{2,i} = \frac{(1-\pi^N) f_{\mu^N_2, \sigma^N_2}(x_i)}{\pi^N f_{\mu^N_1, \sigma^N_1}(x_i) + (1-\pi^N) f_{\mu^N_2, \sigma^N_2}(x_i)} = 1-r_{1,i}.$$
-(M-step) Pour chaque groupe on calcule le poid total de point alloué à un groupe
$$ n_1 = \sum_{i=1}^n r_{1,i}\qquad\text{and}\qquad n_2 = \sum_{i=1}^n r_{2,i}=n-n_1,$$
et on met à jour les différents paramètres du modèle de la manière suivante
$$ \pi^{N+1} = \frac{n_1}{n},$$
$$ \mu^{N+1}_1 = \frac{1}{n_1}\sum_{i=1}^n r_{1,i} x_i\qquad \text{et}\qquad \mu^{N+1}_2 = \frac{1}{n_2}\sum_{i=1}^n r_{2,i} x_i,$$
puis
$$ \sigma^{N+1}_1 = \sqrt{\frac{1}{n_1}\sum_{i=1}^n r_{1,i} (x_i- \mu^{N+1}_1)^2} \qquad \text{et}\qquad \sigma^{N+1}_2 = \sqrt{\frac{1}{n_2}\sum_{i=1}^n r_{2,i} (x_i- \mu^{N+1}_2)^2}.$$
On arrête l'algorithme quand la log-vraissemblance
$$ \mathcal{L}\mathcal{L}^N = \sum_{i=1}^n \log\Big(\pi^N f_{\mu^N_1, \sigma^N_1}(x_i) + (1-\pi^N) f_{\mu^N_2, \sigma^N_2}(x_i)\Big) $$
a convergé. En pratique on fixe $\varepsilon>0$ assez petit, et on arrête l'algorithme quand $$ |\mathcal{L}\mathcal{L}^{N+1} - \mathcal{L}\mathcal{L}^N |<\varepsilon.$$
D'autre critère d'arrêt sont possibles
Question 1) Implémenter l'algorithme précédent à l'aide d'une fonction (appelée EM_gm) qui prend en argument les données, les valeurs initiales pour nos paramètres, et une valeur de $\varepsilon$. Vous pouvez essayer d'utiliser un paramètre correspondant au nombre d'itération maximum que l'on autorise. Si ce nombre est atteint la fonction renvoie un message d'avertissement, disant que l'algorithme n'a pas forcement convergé.¶
def EM_gm(x, pi, mu1, mu2, sigma1, sigma2, err=1e-6, n_iter_max=100):
n = len(x)
ll_old = 0.
ll_new = 0.
# initialisation
r1 = np.zeros(n)
r2 = np.zeros(n)
pi_tmp, mu1_tmp, sigma1_tmp, mu2_tmp, sigma2_tmp = pi, mu1, sigma1, mu2, sigma2
for i in range(n_iter_max):
# E-step
# calcul des probas d'appratenance de chaque point
# à chaque Gaussiene
y = pi * sps.norm(mu1_tmp, sigma1_tmp).pdf(x)
r1 = y / ( y + (1 - pi) * sps.norm(mu2_tmp, sigma2_tmp).pdf(x) )
r2 = 1. - r1
# M-step
# mise à jour des paramètres
n1 = r1.sum()
n2 = n - n1
pi_tmp = n1 / n
mu1_tmp = np.dot(r1, x) / n1
mu2_tmp = np.dot(r2, x) / n2
sigma1_tmp = np.sqrt( np.dot(r1, (x - mu1_tmp)**2 ) / n1 )
sigma2_tmp = np.sqrt( np.dot(r2, (x - mu2_tmp)**2 ) / n2 )
# mise à jour de la log-vraisemblance
ll_new = np.sum( np.log( pi_tmp * sps.norm(mu1_tmp, sigma1_tmp).pdf(x) +
(1. - pi_tmp) * sps.norm(mu2_tmp, sigma2_tmp).pdf(x) ) )
if np.abs(ll_new - ll_old) < err:
break
ll_old = ll_new
else:
print("Attention, le nombre d'itération maximum a été atteint.")
return pi_tmp, mu1_tmp, sigma1_tmp, mu2_tmp, sigma2_tmp
pi, mu1, sigma1, mu2, sigma2 = EM_gm(df.eruptions, 0.5, 3, 3.5, 1., 1.)
pi, mu1, sigma1, mu2, sigma2
(0.3525114684863784, 2.0283755957492695, 0.25104025419862425, 4.2823267029208445, 0.42355791136439475)
xmin = df.eruptions.min()
xmax = df.eruptions.max()
x_plot = np.linspace(xmin, xmax, 100)
def gm_density(pi, mu1, sigma1, mu2, sigma2, x):
return pi * sps.norm(mu1, sigma1).pdf(x) + (1. - pi) * sps.norm(mu2, sigma2).pdf(x)
plt.hist(df.eruptions, bins=40, alpha=0.3, density=True, label="data")
plt.plot(x_plot, gm_density(pi, mu1, sigma1, mu2, sigma2, x_plot), label="estimation")
plt.xlabel("temps d'éruption")
plt.legend();
Question 2: A partir des sorties de la fonction EM_gm écrire une fonction qui agrège les différents points dans les deux catégories. On pourra utiliser les paramètres $r_1$ et $r_2$ pour faire la classification des points. On pourra illustrer les deux groupes graphiquement.¶
def agregation(x, pi, mu1, mu2, sigma1, sigma2):
n = len(x)
r = np.zeros(n)
# calcul des proba d'appartenance à chaque Gaussienne
y = pi * sps.norm(mu1, sigma1).pdf(x)
r = y / ( y + (1 - pi) * sps.norm(mu2, sigma2).pdf(x) )
classif = (r > 1/2)
return classif
classif = agregation(df.eruptions, pi, mu1, mu2, sigma1, sigma2)
plt.scatter(df.eruptions[classif], df.waiting[classif], alpha=0.3, label="short erupt")
plt.scatter(df.eruptions[np.invert(classif)], df.waiting[np.invert(classif)], alpha=0.3, label="long erupt")
plt.xlabel("eruptions")
plt.ylabel("waiting")
plt.legend();
Question 3: Calculer les distributions Bootstrap des paramètres $\pi$, $\mu_1$, $\sigma_1$, $\mu_2$ et $\sigma_2$ (on commencera avec $B=100$ pour faire les test, puis $B=2000$ pour les illustrations). Les illustrer avec des histogrammes. Tracer les densités de notre modèle de mélange correspondantes aux valeurs des échantillons Bootstrap et superposer y les données.¶
Rque: Pour simuler un nombre aléatoire de densité $f$
- on simule une variable $U$ de Bernoulli de paramètre $\pi$
- si $U=1$ on simule une loi normale de paramètre $(\mu_1,\sigma_1^2)$
- si $U=0$ on simule une loi normale de paramètre $(\mu_2,\sigma_2^2)$
def Estim_boot(x, pi, mu1, mu2, sigma1, sigma2, B):
n = len(x)
pi_boot = np.zeros(B)
mu1_boot = np.zeros(B)
sigma1_boot = np.zeros(B)
mu2_boot = np.zeros(B)
sigma2_boot = np.zeros(B)
U = sps.binom(1, pi).rvs((B, n))
m = U.sum(axis=1)
for b in range(B):
Y = np.concatenate([sps.norm(mu1, sigma1).rvs(m[b]), sps.norm(mu2, sigma2).rvs(n- m[b])])
pi_tmp, mu1_tmp, sigma1_tmp, mu2_tmp, sigma2_tmp = EM_gm(Y, pi, mu1, mu2, sigma1, sigma2)
pi_boot[b] = pi_tmp
mu1_boot[b] = mu1_tmp
sigma1_boot[b] = sigma1_tmp
mu2_boot[b] = mu2_tmp
sigma2_boot[b] = sigma2_tmp
return pi_boot, mu1_boot, sigma1_boot, mu2_boot, sigma2_boot
B = 2000
pi_boot, mu1_boot, sigma1_boot, mu2_boot, sigma2_boot = Estim_boot(df.eruptions, pi, mu1, mu2, sigma1, sigma2, B)
plt.hist(pi_boot, bins=50, density=True, alpha=0.3);
plt.hist(mu1_boot, bins=50, density=True, alpha=0.3);
plt.hist(mu2_boot, bins=50, density=True, alpha=0.3);
plt.hist(sigma1_boot, bins=50, density=True, alpha=0.3);
plt.hist(sigma2_boot, bins=50, density=True, alpha=0.3);
plt.hist(df.eruptions, bins=40, alpha=0.3, density=True, label="data")
for b in range(B):
plt.plot(x_plot,
gm_density(pi_boot[b], mu1_boot[b], sigma1_boot[b], mu2_boot[b], sigma2_boot[b], x_plot),
alpha=0.1, color="orange"
)
plt.plot(x_plot, gm_density(pi, mu1, sigma1, mu2, sigma2, x_plot), label="estimation")
plt.xlabel("temps d'éruption")
plt.legend();
On peut donc observer les variations d'estimation de la densité de nos données.
Question 4: On pourra modifier la fonction EM_gm pour enregistrer les différentes valeurs des paramètres $\pi$, $\mu_1$, $\sigma_1$, $\mu_2$ et $\sigma_2$ durant les différentes itérations. On illustrera graphiquement l'évolution le l'adequation du modèle avec les données¶
def EM_gm2(x, pi, mu1, mu2, sigma1, sigma2, err=1e-6, n_iter_max=100):
n = len(x)
ll_old = 0.
ll_new = 0.
# listes pour enregistrer les valeurs succéssives
# sert pour les illustrations
pi_rec = [pi]
mu1_rec = [mu1]
mu2_rec = [mu2]
sigma1_rec = [sigma1]
sigma2_rec = [sigma2]
ll_rec = []
# initialisation
r1 = np.zeros(n)
r2 = np.zeros(n)
pi_tmp, mu1_tmp, sigma1_tmp, mu2_tmp, sigma2_tmp = pi, mu1, sigma1, mu2, sigma2
for i in range(n_iter_max):
# E-step
# calcul des probas d'appratenance de chaque point
# à chaque Gaussiene
y = pi * sps.norm(mu1_tmp, sigma1_tmp).pdf(x)
r1 = y / ( y + (1 - pi) * sps.norm(mu2_tmp, sigma2_tmp).pdf(x) )
r2 = 1. - r1
# M-step
# mise à jour des paramètres
n1 = r1.sum()
n2 = n - n1
pi_tmp = n1 / n
mu1_tmp = np.dot(r1, x) / n1
mu2_tmp = np.dot(r2, x) / n2
sigma1_tmp = np.sqrt( np.dot(r1, (x - mu1_tmp)**2 ) / n1 )
sigma2_tmp = np.sqrt( np.dot(r2, (x - mu2_tmp)**2 ) / n2 )
# Enregistrement des mises à jours
pi_rec.append(pi_tmp)
mu1_rec.append(mu1_tmp)
mu2_rec.append(mu2_tmp)
sigma1_rec.append(sigma1_tmp)
sigma2_rec.append(sigma2_tmp)
# mise à jour de la log-vraisemblance
ll_new = np.sum( np.log( pi_tmp * sps.norm(mu1_tmp, sigma1_tmp).pdf(x) +
(1. - pi_tmp) * sps.norm(mu2_tmp, sigma2_tmp).pdf(x) ) )
# Enregistrement des mises à jours log vraiss
ll_rec.append(ll_new)
if np.abs(ll_new - ll_old) < err:
break
ll_old = ll_new
else:
print("Attention, le nombre d'itération maximum a été atteint.")
return (pi_tmp, mu1_tmp, sigma1_tmp, mu2_tmp, sigma2_tmp,
pi_rec, mu1_rec, sigma1_rec, mu2_rec, sigma2_rec, ll_rec)
(pi, mu1, sigma1, mu2, sigma2,
pi_rec, mu1_rec, sigma1_rec, mu2_rec, sigma2_rec, ll_rec) = EM_gm2(df.eruptions.to_numpy(), 0.5, 3, 3.5, 1., 1.)
plt.hist(df.eruptions, bins=40, alpha=0.3, density=True, label="data")
m = len(pi_rec)
ind = [0, 5, 10, 15]
for l in ind:
plt.plot(x_plot, gm_density(pi_rec[l], mu1_rec[l], sigma1_rec[l], mu2_rec[l], sigma2_rec[l], x_plot),
alpha=0.5, label="itération %s" % l)
plt.plot(x_plot, gm_density(pi, mu1, sigma1, mu2, sigma2, x_plot), alpha=0.5, label="estimation")
plt.xlabel("temps d'éruption")
plt.legend();
On voit ici que la convergence est très rapide, après 15 itérations, on avait quasiment convergé.