#!/usr/bin/env python
# coding: utf-8
#
#
# # TP Algorithmes Stochastiques, SEM, N-SEM, EM
# In[1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sps
get_ipython().run_line_magic('matplotlib', 'inline')
# In[2]:
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}
# $
#
#
# ## Rappels de cours
# ### Mélange de gaussiennes
# 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.
# In[3]:
# 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
# ### Simulation des données
#
# 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 :
# - $Z\in \{1, \ld, K\}$ de loi $(\al_1, \ld, \al_K)$,
# - pour tout $k\in \{1, \ld, K\}$, conditionellement à $Z=k$, $X\sim \cN(\mu_k, \si_k^2)$.
#
# 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.
# In[4]:
# 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")
# ### Estimation des paramètres dans le cas où l'on connaît les $z_i$
# 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.
# In[5]:
(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))
#
# ## Algorithmes
#
# 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.
#
# ### Algorithme SEM
# - initialisation arbitraire des classes $z_i\in \{1,\ldots,K\}$
# - répéter un grand nombre de fois:
# - a) calculer $\widehat{\theta}$
# - b) pour tout $i=1,\ldots,n$, tirer la valeur de $z_i$ selon la loi $$\mathbb{P}_{\widehat\theta}(Z_i=\bullet | X_i=x_i)$$
#
# 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.
# In[6]:
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))
# ### Algorithme SEM à $N$ tirages
# - dupliquer $N$ fois le jeu d'observations $(x_1, \ld, x_n)$, qui devient donc $(x_{i}^{(j)})_{1\le i\le n,\, 1\le j\le N}$
# - appliquer SEM à ce jeu de données étendu
#
# 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.
# In[7]:
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))
# ### Algorithme EM
# - initialisation arbitraire du paramètre $\theta_0$,
# - étant donné le paramètre $\theta_t$, répéter pour $t=0,1,2,\dots$ et jusqu'à convergence
# - Calcul de la matrice
# $$\left[\p_{\tta_{t}}(Z=k | X=x_i)\right]_{1\le i\le n, \, 1\le k\le K}=\left[\f{\al_k^tf_{\mu_k^t,\si^t_k}(x_i)}{\sum_{l=1}^K\al_l^tf_{\mu_l^t,\si_l^t}(x_i)}\right]_{1\le i\le n, \, 1\le k\le K}$$
# - Calcul de $\tta_{t+1}$ : pour tout $k=1,\ld,K$,
# $$\al_k^{t+1}=\ff n\sum_{i=1}^n \p_{\tta_{t}}(Z=k | X=x_i),$$
# $$\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.
# In[8]:
print("\nEM :")
L=100
(alpha_est,mu_est,sigma_est)=EM(L,K,X)
# In[9]:
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))
# ## EM en dimension supérieure
# 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 :
# - K=3
# - alpha=(.5,.3,.2)
# - les 3 gaussiennes ont les moyennes $\mu_1=(-5,-5)$, $\mu_2 = (5,5)$, $\mu_3=(0,0)$. On pourra rassembler ces moyennes dans une liste $mu=[[-5,-5],[5,5],[0,0]]$
# - les 3 gaussiennes ont la même covariance $$\begin{pmatrix} 1 & 0.5 \\ 0.5 & 1\end{pmatrix}.$$
#
# 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.
# In[10]:
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)