<jupyter_start><jupyter_text>Loan Sarazin & <NAME>  Mise en oeuvre de l'algorithme EM  Calcul de la valeur de la densité de probabilité d'un mélange de gaussienne en un point Les paramétres de la fonction gm_pdfx : le point où l'on calcule la valeur de la densitèmu : le vecteur des moyennes des gaussiennessigma : le vecteur des écart-types des gaussiennesp : le vecteur des probabilitès de la loi multinomiale associèe<jupyter_code>import numpy as np
import pandas as pd
from scipy.stats import norm
from scipy.stats import uniform
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
def gm_pdf(x, mu, sigma, p):
    #Initialisation de la variable de sortie
    resultat = 0.0
    #Contrôle de la cohérence des paramètres d'entrée
    #Le vecteur de moyenne doit avoir la même longueur que le vecteur p
    if len(mu) != len(p):
        print('Erreur de dimension sur la moyenne')
    # Le vecteur des écart-types doit avoir la même longueur que le vecteur p
    elif len(sigma) != len(p):
            print('Erreur de dimension sur les écarts-types')
    else:
    # Calcul de la valeur de la densité
        for i in range(0, len(p)):
            resultat = resultat + p[i] * norm.pdf(x, mu[i], sigma[i])
    return resultat<jupyter_output><empty_output><jupyter_text>Génération de nombre aléatoire suivant un mélange de gaussienneDans un mélange de gaussienne chaque densité de probabilité correspondant à une loi normale est pondérée par un coefficient plus petit que un.La somme de tous ces coefficients est égale à un. Une variable aléatoire qui suit une loi normale suit donc une des lois normales du mélange avec une probabilité égale au coefficient de pondération de cette même loi normale. Pour générer un échantillon suivant une loi normale, il faut donc procéder en deux étapes :Tirer aléatoirement un nombre entre 1 et N (le nombre de gaussiennes du mélange) suivant une loi mumltinomiale définie par les coefficients du mélange.Une fois ce nombre obtenu, on génère le nombre suivant la loi normale associée. Les paramètres de la fonction gm_rndmu : le vecteur des moyennes des gaussiennessigma : le vecteur des écart-types des gaussiennesp : le vecteur des probabilités de la loi multinomiale associée<jupyter_code>def gm_rnd(mu, sigma, p):
    # Initialisation de la variable de sortie
    resultat = 0.0
    #Contrôle de la cohérence des paramètres d'entrée
    #Le vecteur de moyenne doit avoir la même longueur que le vecteur p
    if len(mu) != len(p):
        print('Erreur de dimension sur la moyenne')
    # Le vecteur des écart-types doit avoir la même longueur que le vecteur p
    elif len(sigma) != len(p):
            print('Erreur de dimension sur sur les écarts-types')
    else:
    #Génération de l'échantillon
    # On échantillonne suivant une loi uniforme sur [0,1]
        u = uniform.rvs(loc = 0.0, scale = 1.0, size = 1)
    # % Chaque test suivant permet de définir un intervalle sur lequel la
    # probabilité d'appartenance de la variable uniforme est égale à l'une des
    # probabilités définie dans le vecteur p. Lorsque u appartient à l'un de
    # ces intervalles, c'est équivalent à avoir générer une variable aléatoire
    # suivant l'un des éléments de p. Par exemple, pour le premier test
    # ci-dessous, la probabilité que u appartienne à l'intervalle [0,p[0][ est
    # égale à p[0] puisque u suit une loi uniforme. Donc si u appartient à
    # [0,p[0][ cela est équivalent à avoir tirer suivant l'événement de probabilité p[0].
        if u < p[0]: # On test si on a généré un événement de probabilité p[0]
            resultat = sigma[0] * norm.rvs(loc = 0, scale = 1, size = 1) + mu[0]
            # Pour générer suivant une loi normale quelconque, il suffit de multiplier
            # une variable normale centrée réduite (moyenne nulle et écart-type égal à 1)
            # par l'écart-type désité et d'additionner la moyenne désirée au produit précédent.
        for i in range(1, len(p)):
            if (u > np.sum(p[0:i])) and (u <= np.sum(p[0:i+1])): # On test si on a généré
                # un événement de probabilité p[i]
                resultat = sigma[i] * norm.rvs(loc = 0.0, scale = 1.0, size = 1) + mu[i]
                # Pour générer suivant une loi normale quelconque, il suffit de multiplier
                # une variable normale centrée réduite (moyenne nulle et écart-type égal à 1)
                # par l'écart-type désité et d'additionner la moyenne désirée au produit précédent.
    return resultat<jupyter_output><empty_output><jupyter_text>Les paramètres de la densité de mélange de gaussienne<jupyter_code>p = np.array([0.2, 0.5, 0.3])
mu = np.array([-5.0, 0.0, 6.0])
sigma = np.array([1.0, 0.5, 1.0])<jupyter_output><empty_output><jupyter_text>Calcul de la densité de probabilité du mélange de gaussienne<jupyter_code>x = np.arange(-20,20, 0.001)
pointPdf = gm_pdf(x, mu, sigma, p)<jupyter_output><empty_output><jupyter_text>Tracé de la densité de probabilité du mélange de gaussienne<jupyter_code>plt.plot(x, pointPdf, 'r-')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Génération d'échantillons suivant  densité de probabilité du mélange de gaussienne<jupyter_code>nbEchantillon = 1000
donnees = np.array(nbEchantillon*[0],dtype = float)
for i in range(0, nbEchantillon, 1):
    donnees[i] = gm_rnd(mu, sigma, p)<jupyter_output><empty_output><jupyter_text>Tracé des échantillons<jupyter_code>plt.plot(donnees, 'b')
plt.title('Melange de gausiennes')
plt.xlabel('Index')
plt.ylabel('Donnees')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Tracé des échantillons<jupyter_code>plt.plot(donnees, 'g.')
plt.title('Melange de gausiennes')
plt.xlabel('Index')
plt.ylabel('Donnees')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Tracé de l'histogramme des échantillons<jupyter_code>plt.hist(donnees, bins = 30, density = False, edgecolor = "red")
plt.title('Melange de gausiennes')
plt.xlabel('Donnees')
plt.ylabel('Index')
plt.show()<jupyter_output><empty_output><jupyter_text>Tracé de l'histogramme des échantillons<jupyter_code>plt.hist(donnees, bins = 30, density = True, color = 'yellow', edgecolor = "red")
plt.title('Melange de gausiennes')
plt.xlabel('Donnees')
plt.ylabel('Index')
plt.show()<jupyter_output><empty_output><jupyter_text>L'algorithme EM<jupyter_code>nbMaxIterations = 40
mu_em = np.array([-0.0156, -4.9148, 5.9692])
sigma_em = np.array([1.3395, 1.3395, 1.3395])
alpha_em = np.array([0.4800, 0.2200, 0.3000])

def EM_algorithm(nbMaxIterations, mu_em, sigma_em, alpha_em, donnees):
    nbIteration = 1 #Initialisation de la variable d'arrêt
    nbComposante = len(alpha_em) #Nombre de composantes du mélange
    nbDonnees = len(donnees)  #Nombre de données
    p = np.zeros(shape=(nbComposante, nbDonnees))
    #Déclaration et initialisation de la matrice qui va contenir les probabilités
    #p(k|x,theta_courant)

    alpha_em_new = alpha_em
    sigma_em_carre_new = sigma_em
    mu_em_new = mu_em
    donneesP = np.zeros(shape=(nbEchantillon))

    while nbIteration < nbMaxIterations:
        #Création d'un array avec les lois à postériori
        for n in range(0, nbDonnees, 1):
            for k in range(0, nbComposante, 1):
                p[k, n] = alpha_em[k] * norm.pdf(x = donnees[n], loc = mu_em[k], scale = sigma_em[k])
            p[:, n] = p[:, n] / np.sum(p[:, n])
        for k in range(0, nbComposante, 1):
            alpha_em_new[k] = np.sum(p[k,:]) / nbDonnees #Moyenne de la loi de la gaussienne k
            for n in range(0, nbDonnees, 1):
                donneesP[n] = donnees[n] * p[k, n]
            mu_em_new[k]  = np.sum(donneesP) / np.sum(p[k, :])
            for n in range(nbDonnees):
                donneesP[n] = ((donnees[n] - mu_em_new[k]) ** 2) * p[k, n]
            sigma_em_carre_new[k] = np.sum(donneesP) / np.sum(p[k, :])
        mu_em = mu_em_new
        sigma_em = np.sqrt(sigma_em_carre_new)
        alpha_em = alpha_em_new
        nbIteration = nbIteration + 1
    return mu_em, sigma_em, alpha_em

mu_em, sigma_em, alpha_em = EM_algorithm(nbMaxIterations, mu_em, sigma_em, alpha_em, donnees)
print('Les paramètres estimés sont : ')
print('Moyennes des composantes du mélange', mu_em)
print('Ecrat type des composantes du mélange', sigma_em)
print('Probabilités des composantes du mélange', alpha_em)
print('La somme des probabilités des composantes du mélange vaut : ', np.sum(alpha_em))
plt.plot(x, pointPdf, 'r-', label = 'Originale')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')
pointPdfEstime = gm_pdf(x, mu_em, sigma_em, alpha_em)
plt.plot(x, pointPdfEstime, 'b-', label = 'Estimée')
plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()
#Comparaison en faisant varier les valeurs d'initialisation
plt.plot(x, pointPdf, 'r-', label = 'Originale')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')
pointPdfEstime = gm_pdf(x, mu_em, sigma_em, alpha_em)
plt.plot(x, pointPdfEstime, 'b-', label = 'Estimée init')
plt.legend(loc='upper left', shadow=True, fontsize='x-large')


mu_1 = np.array([0, -5, 6])
sigma_1 = np.array([1.3395, 1.3395, 1.3395])
alpha_1 = np.array([0.4800, 0.2200, 0.3000])
mu_1, sigma_1, alpha_1 = EM_algorithm(nbMaxIterations, mu_1, sigma_1, alpha_1, donnees)

pointPdfEstime = gm_pdf(x, mu_1, sigma_1, alpha_1)
plt.plot(x, pointPdfEstime, 'b-', label = 'Estimée 1')
plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()
#Comparaison en faisant varier les valeurs d'initialisation - nombre d'itérations
#Conservation des moyennes arrondies
mu_1 = np.array([0, -5, 6])
sigma_1 = np.array([1.3395, 1.3395, 1.3395])
alpha_1 = np.array([0.4800, 0.2200, 0.3000])

mu_2, sigma_2, alpha_2 = EM_algorithm(10, mu_1, sigma_1, alpha_1, donnees)

mu_3, sigma_3, alpha_3 = EM_algorithm(50, mu_1, sigma_1, alpha_1, donnees)



plt.figure(figsize = (10, 5))

plt.plot(x, pointPdf, 'r-', label = 'Originale')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')

pointPdfEstime = gm_pdf(x, mu_2, sigma_2, alpha_2)
plt.plot(x, pointPdfEstime, label = 'Itérations = 10')


pointPdfEstime = gm_pdf(x, mu_3, sigma_3, alpha_3)
plt.plot(x, pointPdfEstime, label = 'Itérations = 40')

plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()

#Comparaison en faisant varier les valeurs d'initialisation - nombre d'itérations
#Retour aux valeurs initiales des moyennes
mu_1 = np.array([-0.0156, -4.9148, 5.9692])
sigma_1 = np.array([1.3395, 1.3395, 1.3395])
alpha_1 = np.array([0.4800, 0.2200, 0.3000])

mu_2, sigma_2, alpha_2 = EM_algorithm(10, mu_1, sigma_1, alpha_1, donnees)

mu_3, sigma_3, alpha_3 = EM_algorithm(50, mu_1, sigma_1, alpha_1, donnees)
plt.figure(figsize = (10, 5))

plt.plot(x, pointPdf, 'r-', label = 'Originale')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')

pointPdfEstime = gm_pdf(x, mu_2, sigma_2, alpha_2)
plt.plot(x, pointPdfEstime, label = 'Itérations = 10')


pointPdfEstime = gm_pdf(x, mu_3, sigma_3, alpha_3)
plt.plot(x, pointPdfEstime, label = 'Itérations = 40')

plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()<jupyter_output><empty_output><jupyter_text>Conclusion On peut voir qu'en faisant varier très légèrement les paramètres d'initialisation, les estimations diffèrent beaucoup. Le nombre d'itérations a une influence négligeable lorsque les paramètres sont optimaux. En revanche, si les paramètres d'initialisation sont légèrement modifiés, et donc ne sont plus optimaux, le nombre d'itérations a une influence forte sur les résultats d'estimation.   Algorithme EM avec les données Galaxy<jupyter_code>data = pd.read_excel('2122_Galaxy.xlsx', header = None)
donnees = data.to_numpy()

plt.hist(donnees, bins = 30, density = False, edgecolor = "red")
plt.title('Données Galaxy')
plt.xlabel('Donnees')
plt.ylabel('Vitesse')
plt.show()
nb_vel = donnees.size
mu, sigma = np.mean(donnees), np.var(donnees)
print("La moyenne des vitesses vaut : {:.2f}".format(mu))
print("La variance des vitesses vaut : {:.2f}".format(sigma))

vel_min = donnees.min()
vel_max = donnees.max()
donnees.shape<jupyter_output><empty_output><jupyter_text>Tracé des vitesses<jupyter_code>plt.figure()
plt.plot(donnees, '.')<jupyter_output><empty_output><jupyter_text>Détermination du nombre de gaussiennesUtilisation de la fonction gaussian Mixture<jupyter_code>#3 gaussiennes
gm = GaussianMixture(n_components=3).fit(donnees)
# la méthode fit de GaussianMixture utilise l'algorithme EM.

mu = gm.means_.ravel()
alpha = gm.weights_
sigma = np.sqrt(gm.covariances_)
print("Les moyennes des Gaussiennes composantes sont les suivantes : \n")
print(mu)
print("\nLes poids alpha sont les suivants : \n")
print(alpha)
x = np.arange(9000,35000, 1000)
pointPdfEstime = gm_pdf(x, mu, sigma, alpha)

plt.figure()
plt.hist(donnees, bins = 30, density = True, edgecolor = "red")
plt.plot(x, pointPdfEstime.T, 'b-', label = 'Estimation')
plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()
#2 Gaussiennes
gm = GaussianMixture(n_components=4, covariance_type='full').fit(donnees.reshape(-1, 1))

nbMaxIterations = 40
mu_em_mod = np.round(gm.means_.ravel(), 4)
sigma_em_mod = np.round(np.sqrt(gm.covariances_.T.ravel()), 4)
alpha_em_mod = gm.weights_

print('Les paramètres estimés sont : ')
print('Moyennes des composantes du mélange', mu_em_mod)
print('Ecart type des composantes du mélange', sigma_em_mod)
print('Probabilités des composantes du mélange', alpha_em_mod)

pointPdf = gm_pdf(x, mu_em_mod, sigma_em_mod, alpha_em_mod)

plt.figure
plt.hist(donnees, bins = 40, density = True, color='skyblue', edgecolor = "red", label = 'Vitesses')
plt.plot(x, pointPdf, 'b-', label = 'Estimée')
plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.xlabel('Donnees')
plt.ylabel('pdf')
plt.title('Melange de gausiennes')
plt.show()<jupyter_output>Les paramètres estimés sont : 
Moyennes des composantes du mélange [19887.2015 33044.334  23203.2079  9710.1429]
Ecart type des composantes du mélange [1320.9724  921.7177 1606.6029  422.5107]
Probabilités des composantes du mélange [0.4764582  0.03658534 0.4015906  0.08536585]<jupyter_text>Comparaison des résultats en fonction du nombre de Gaussiennes<jupyter_code>x = np.arange(vel_min-100, vel_max+100, 5)

nb_gauss = [2, 3, 4, 5,6]
plt.figure(figsize = (15, 8))
plt.hist(donnees, bins = 40, density = True, color='skyblue', edgecolor = "red", label = 'Vitesses')

for i in nb_gauss:
    gm = GaussianMixture(n_components=i, covariance_type='full', n_init = 10).fit(donnees.reshape(-1, 1))

    nbMaxIterations = 40
    mu_em_mod = np.round(gm.means_.ravel(), 4)
    sigma_em_mod = np.round(np.sqrt(gm.covariances_.T.ravel()), 4)
    alpha_em_mod = gm.weights_

    pointPdf = gm_pdf(x, mu_em_mod, sigma_em_mod, alpha_em_mod)
    
    plt.plot(x, pointPdf, label = f'Estimée, {i} gaussiennes', linewidth = 3)
    
plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.xlabel('Donnees')
plt.ylabel('pdf')
plt.title('Melange de gausiennes')
plt.show()<jupyter_output><empty_output><jupyter_text>On peut voir sur la graphique ci-dessus que si on augmente trop le nombre de gaussiennes, on donne un poids très important à des données qui sont non significatives. Par exemple, avec 6 gaussiennes, on obtient un pic principal autour d'une valeur concentrant peu de données. Quand on utilise un nombre faible de gaussiennes, l'estimation est très proche de la réalité pour la majorité des données. Cependant, l'estimation ne prend pas en compte les données extrêmes ou les petits groupes de données. Il s'agit donc de trouver un compromis. D'après le graphique ci-dessus, on pourrait conclure que le nombre idéal de Gaussiennes à utiliser est de quatre. En effet, à partir de 5 gaussiennes, un pic se créé autour d'un petit groupe de données. Inversement, en dessous de 4 gaussiennes, l'estimation s'éloigne de la réalité pour les données 'centrales' (c'est-à-dire le plus gros groupe de données) en faisant disparaître les deux groupes de données.  Algorithme EM appliqué aux données Galaxy<jupyter_code>#Détermination des paramètres d'initialisation de manière empirique (lecture histogramme et exécution de l'algorithme EM)

moy_em_emp = np.array([9000, 20000, 23000, 33000])
sum_moy = np.sum(moy_em_emp)
ecart_type_em_emp = np.array([400, 1000, 1000, 400])
alpha_em_emp =moy_em_emp/sum_moy
#Algo EM du TP pour les paramètres définis ci-dessus
mu_emp, ecart_type_emp, alpha_emp = EM_algorithm(40, moy_em_emp, ecart_type_em_emp, alpha_em_emp, donnees)

#Gaussian mixture pour 4 gaussiennes
gm = GaussianMixture(n_components=4, covariance_type='full').fit(donnees.reshape(-1, 1))

nbMaxIterations = 40
mu_em_mod = np.round(gm.means_.ravel(), 4)
sigma_em_mod = np.round(np.sqrt(gm.covariances_.T.ravel()), 4)
alpha_em_mod = gm.weights_
pointPdf = gm_pdf(x, mu_em_mod, sigma_em_mod, alpha_em_mod)
plt.figure(figsize = (15, 8))

plt.hist(donnees, bins = 40, density = True, color='skyblue', edgecolor = "red", label = 'Vitesses')
plt.title('Densite de melange')
plt.xlabel('x')
plt.ylabel('pdf')

pointPdfEstime = gm_pdf(x, mu_emp, ecart_type_emp, alpha_emp)
plt.plot(x, pointPdfEstime, label = 'Itérations = 40')

plt.plot(x, pointPdf, 'b-', label = 'Estimation "optimale" Gaussian mixture')

plt.legend(loc='upper left', shadow=True, fontsize='x-large')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.117198256080771, 'fraction_numerical': 0.03051858650757228, 'mean_word_length': 4.704515706806283, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '42958', 'license': 'MIT', 'n_tokens_mistral': 7107, 'n_tokens_neox': 6658, 'n_words': 1809, 'path': 'TP2/MARIZY_Galaxy.ipynb', 'repo_name': 'LoanSarazin/sdia-estimation', 'size': '17437'}