import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sps
%matplotlib inline
from TP_SEM_EM import *
$\renewcommand{\si}{\sigma} \renewcommand{\al}{\alpha} \renewcommand{\tta}{\theta} \renewcommand{\Tta}{\Theta} \renewcommand{\Si}{\Sigma} \renewcommand{\ld}{\ldots} \renewcommand{\cd}{\cdots} \renewcommand{\cN}{\mathcal{N}} \renewcommand{\R}{\mathbb{R}} \renewcommand{\p}{\mathbb{P}} \renewcommand{\f}{\frac} \renewcommand{\ff}{\frac{1}} \renewcommand{\ds}{\displaystyle} \renewcommand{\bE}{\mathbf{E}} \renewcommand{\bF}{\mathbf{F}} \renewcommand{\ii}{\mathrm{i}} \renewcommand{\me}{\mathrm{e}} \renewcommand{\hsi}{\hat{\sigma}} \renewcommand{\hmu}{\hat{\mu}} \renewcommand{\ste}{\, ;\, } \renewcommand{\op}{\operatorname} \renewcommand{\argmax}{\op{argmax}} \renewcommand{\lfl}{\lfloor} \newcommand{\ri}{\right} $
On rappelle qu'une mixture gaussienne a une densité de la forme $$x\mapsto \sum_{k=1}^K\alpha_kf_{\mu_k,\sigma_k}(x)=\sum_{k=1}^K\alpha_k\frac{1}{\sqrt{2\pi \sigma_k^2}}\mathrm{e}^{-\frac{(x-\mu_k)^2}{2\sigma_k^2}}.$$
On commence par définir une fonction qui permet d'évaluer la densité théorique d'un mélange de gaussiennes sur une grille x.
# Calcul de la densité théorique d'un GMM sur une grille x
def densite_theorique(K,mu,sigma,alpha,x):
y=np.zeros(len(x))
for j in range(K):
y+=alpha[j]*sps.norm.pdf(x,loc=mu[j],scale=sigma[j])
return y
Pour simuler une mixture gaussienne, on commence par fixer $K\geq 1$, et par choisir un paramètre $\theta:=(\alpha_k,\mu_k,\sigma_k)_{1\le k\le K}$, avec $\mu_1, \ldots, \mu_K\in \R$, $\sigma_1, \ldots, \sigma_K>0$ et $\alpha_1, \ldots, \al_K>0$ tels que $\al_1+\cd+\al_K=1$. On tire au hasard des données $(z_1,x_1), \ld,(z_n,x_n )$, copies indépendantes du vecteur aléatoire $(Z,X)$ de loi $\p_\theta$ donnée par :
Notons que le passage de $X$ à $Z$ se fait avec la formule de Bayes : pour tout $x\in \R$ et $k\in \{1, \ld, K\}$,
$$\label{eq:Oleron1}\p_\theta(Z=k | X=x)\;=\;\f{\al_kf_{\mu_k,\si_k}(x)}{\sum_{l=1}^K\al_lf_{\mu_l,\si_l}(x)}.$$
avec $$f_{\mu,\si}(x):=\ff{\sqrt{2\pi}\si}\me^{-\f{(x-\mu)^2}{2\si^2}}.$$
Les algorithmes SEM, $N$-SEM, EM ont pour but d'estimer les paramètres $\mu_1, \ld, \mu_K\in \R$, $\si_1, \ld, \si_K>0$ et $\al_1, \ld, \al_K>0$ ainsi que les valeurs des $z_i$ par la seule observation des $x_i$.
1. Ecrire une fonction echantillon(K,alpha,mu,sigma,n)
permettant de simuler $n$ réalisations indépendantes d'une telle loi. $K$ est le nombre de gaussiennes, alpha le vecteur des paramètres du mélange, mu les moyennes et sigma les écart-types (ce sont donc des listes de $K$ élements). La fonction doit renvoyer à la fois les Z_i et les X_i, par exemple sous forme de deux listes. Essayer de coder cette fonction SANS BOUCLE FOR.
2. Tester cette fonction avec K=3, alpha=(.4,.3,.3), mu=(-4,4,0), sigma=(1,1,1) et n = 1000 et tracer sur le même graphe l'histogramme des réalisations ainsi obtenues et le mélange correspondant.
# Simulation des Données
K=3
alpha_real=[.4,.3,.3]
mu_real=[-4,4,0]
sigma_real=[1,1,1]
n=1000
Z,X=echantillon(K,np.array(mu_real),np.array(sigma_real),alpha_real,n)
plt.hist(X,bins=round(3*n**(.3)),density=1,histtype="step",color = 'b') # remplacer density par normed si matplotlib <2.2
x=np.linspace(-8,8,num=1000)
y=densite_theorique(K,mu_real,sigma_real,alpha_real,x)
plt.plot(x,y,"r")
[<matplotlib.lines.Line2D at 0x137ca8c40>]
L'estimateur du Maximum de Vraisemblance (EMV), étant données les observations $(z_i,x_i)_{1\le i\le n}$, fournit les estimations suivantes pour le paramètre $\tta=(\al_k,\mu_k,\si_k)_{1\le k\le K}$: pour tout $k=1,\ld, K$,
$$\widehat{\al}_k=\f{|\{i\ste z_i=k\}|}{n},$$
$$\widehat{\mu}_k=\ff{n\widehat{\al}_k}\sum_{i\ste z_i=k} x_i,$$
$$\widehat{\si}^2_k=\ff{n\widehat{\al}_k}\sum_{i\ste z_i=k} (x_i-\widehat{\mu}_k)^2.$$
3. Ecrire une fonction estimation(K,Z,X) qui calcule l'EMV dans le cas où l'on connait les $(z_i,x_i)_{1\le i\le n}$. Cette fonction doit retourner les trois listes de taille K (alpha,mu,sigma). Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange.
(alpha_estime,mu_estime,sigma_estime)=estimation(K,Z,X)
print("\nVrai alpha et alpha estimé :")
print(alpha_real,alpha_estime)
print("\nVrai mu et mu estimé :")
print(mu_real,np.around(mu_estime,decimals=3))
print("\nVrai sigma et sigma estimé :")
print(sigma_real,np.around(sigma_estime,decimals=3))
Vrai alpha et alpha estimé : [0.4, 0.3, 0.3] [0.386, 0.282, 0.332] Vrai mu et mu estimé : [-4, 4, 0] [-3.981 3.997 0.042] Vrai sigma et sigma estimé : [1, 1, 1] [0.996 1.056 0.98 ]
Dans la suite de ce TP, on souhaite appliquer les algorithmes SEM, N-SEM et EM aux données $x_1, \ld, x_n$ pour estimer $\tta$. Estimer alors les $z_i$ : pour tout $i$, on choisit $$z_i:=\argmax_k \p_\theta(Z=k | X=x_i).$$ Comparer alors les valeurs de $\tta$ ainsi que des $z_i$ à leurs vraies valeurs.
4. Ecrire une fonction SEM(nbpas,K,x)
qui implémente l'algorithme précédent avec L étapes. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange.
L=int(100)
print("\nSEM :")
alpha_est,mu_est,sigma_est =0,0,0
(alpha_est,mu_est,sigma_est)=SEM(L,K,X)
print("\nVrai alpha et alpha estimé :")
print(alpha_real,alpha_est)
print("\nVrai mu et mu estimé :")
print(mu_real,np.around(mu_est,decimals=3))
print("\nVrai sigma et sigma estimé :")
print(sigma_real,np.around(sigma_est,decimals=3))
SEM : Vrai alpha et alpha estimé : [0.4, 0.3, 0.3] [1. 1. 1.] Vrai mu et mu estimé : [-4, 4, 0] [ 3.885 -3.938 0.059] Vrai sigma et sigma estimé : [1, 1, 1] [1.153 1.014 0.833]
5. Ecrire une fonction NSEM(N,L,K,x)
qui implémente l'algorithme précédent. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange.
L=100
N=10
print("\nN-SEM :")
(alpha_est,mu_est,sigma_est)=NSEM(N,L,K,X)
print("\nVrai alpha et alpha estimé :")
print(alpha_real,alpha_est)
print("\nVrai mu et mu estimé :")
print(mu_real,np.around(mu_est,decimals=3))
print("\nVrai sigma et sigma estimé :")
print(sigma_real,np.around(sigma_est,decimals=3))
N-SEM : Vrai alpha et alpha estimé : [0.4, 0.3, 0.3] [1. 1. 1.] Vrai mu et mu estimé : [-4, 4, 0] [ 3.94 -3.947 0.066] Vrai sigma et sigma estimé : [1, 1, 1] [1.104 1.014 0.869]
$$\mu_k^{t+1}=\ff{n\al_k^{t+1}}\sum_{i=1}^n \p_{\tta_{t}}(Z=k | X=x_i)x_i,$$ $$(\si_k^{t+1})^2=\ff{n\al_k^{t+1}}\sum_{i=1}^n \p_{\tta_{t}}(Z=k | X=x_i)(x_i-\mu_k^{t+1})^2.$$
6. Ecrire une fonction EM(L,K,x)
qui implémente l'algorithme précédent. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange.
print("\nEM :")
L=100
(alpha_est,mu_est,sigma_est)=EM(L,K,X)
EM :
print("\nVrai alpha et alpha estimé :")
print(alpha_real,alpha_est)
print("\nVrai mu et mu estimé :")
print(mu_real,np.around(mu_est,decimals=3))
print("\nVrai sigma et sigma estimé :")
print(sigma_real,np.around(sigma_est,decimals=3))
Vrai alpha et alpha estimé : [0.4, 0.3, 0.3] [0.33462064 0.31353918 0.35184018] Vrai mu et mu estimé : [-4, 4, 0] [-4.09 2.332 0.687] Vrai sigma et sigma estimé : [1, 1, 1] [0.934 2.286 2.523]
On se place cette fois en dimension deux. Chaque gaussienne du mélange est déterminée par sa moyenne (un vecteur de $\R^2$ et sa matrice de covariance (une matrice $2\times 2$).
7. Ecrire une nouvelle fonction échantillon2d
permettant d'échantillonner $n$ fois le mélange de gaussiennes de dimension 2 ayant les caractéristiques suivantes :
8. Appliquer la fonction précédentes pour réaliser une liste $X$ constituée de $n$ observations du mélange (donc $X$ est une liste de vecteurs de dimension 2).
9. Ecrire une nouvelle fonction EM2d
permettant d'estimer les paramètres du mélange à partir des observations $X$ et comparer le résultat obtenu avec les paramètres réels du mélange.
d=2
K=3
alpha=[.5,.3,.2]
mu=[[-5,-5],[5,5],[0,0]]
covariance=[np.array([[1,.5],[.5,1]])]*K
n=int(5e2)
L=100
# on échantillonne les données
(z0,x0)=echantillon2d(K,alpha,mu,covariance,n)
# on applique l'algorithme EM et on vérifie que les paramètres estimés sont proches des vrais paramètres
(alpha_em,mu_em,covariance_em)=EM2d(3,L,x0)
print("attention les classes estimées peuvent être dans un ordre différent des classes réelles !")
print("\nVrai alpha :",alpha)
print("\nEstimation de alpha :", alpha_em)
print("\nVrai mu :",mu)
print("\nEstimation de mu:",mu_em)
print("\nVraie covariance :",covariance)
print("\nEstimation de la covariance :", covariance_em)
attention les classes estimées peuvent être dans un ordre différent des classes réelles ! Vrai alpha : [0.5, 0.3, 0.2] Estimation de alpha : [0.23196642 0.2688371 0.49919648] Vrai mu : [[-5, -5], [5, 5], [0, 0]] Estimation de mu: [[-0.02430919 0.18438086] [ 4.9245234 4.91653035] [-5.09633677 -4.9550138 ]] Vraie covariance : [array([[1. , 0.5], [0.5, 1. ]]), array([[1. , 0.5], [0.5, 1. ]]), array([[1. , 0.5], [0.5, 1. ]])] Estimation de la covariance : [[[1.2292515 0.64568701] [0.64568701 1.28573749]] [[0.90058606 0.33126681] [0.33126681 0.82164365]] [[0.95682149 0.50140257] [0.50140257 1.13658659]]]