donnees_n;

% Tracé des données d'apprentissage (croix bleues) :
figure('Name','Estimation par le maximum de vraisemblance','Position',[0.33*L,0,0.33*L,0.5*H]);
plot(D_app(1,:),D_app(2,:),'+b','MarkerSize',10,'LineWidth',2);

xlabel('$x$','Interpreter','Latex');
ylabel('$y$','Interpreter','Latex');
axis([-taille taille -taille taille]);
xlim([-taille taille]);
ylim([-taille taille]);
axis equal;
hold on;

% Tirages aléatoires de paramètres pour la paire d'ellipses :
nb_tirages = 20000;
parametres_test = zeros(nb_tirages,n,5);
parametres_test(:,:,1) = 2*taille/5*(rand(nb_tirages,n)+1);		    % Demi-grand axe
parametres_test(:,:,2) = rand(nb_tirages,n);				        % Excentricité
parametres_test(:,:,3) = (3*taille/5)*(4*rand(nb_tirages,n)-2);		% Abscisse du centre
parametres_test(:,:,4) = (3*taille/5)*(4*rand(nb_tirages,n)-2);		% Ordonnée du centre
parametres_test(:,:,5) = 2*pi*rand(nb_tirages,n);			        % Angle du grand axe

% Estimation d'une paire d'ellipses par le maximum de vraisemblance :
parametres_estim = max_vraisemblance_n(D_app,parametres_test,sigma,n);
parametres_estim = reshape(parametres_estim,n,5);

% Tracé des ellipses estimées par le maximum de vraisemblance (trait rouge) :
for i = 1:n
    [x,y] = points_ellipse(parametres_estim(i,:),theta_affichage);
    plot([x x(1)],[y y(1)],'r-','LineWidth',3);
end
drawnow;

pause(0.5);

% Algorithme EM :
difference_score = 1;
seuil = 0.00001;
score = calcul_score_n(parametres_VT,parametres_estim,n);
k = 0;
while abs(difference_score)>seuil

	% Calcul des probabilités d'appartenance aux deux classes :
	probas = probabilites_EM_n(D_app,parametres_estim,n,sigma);

	% Partition des données :
    [~,classes] = max(probas);
    
    hold off;
    for i = 1:n
        D_app_i = D_app(:,classes == i);
	    plot(D_app_i(1,:),D_app_i(2,:),'+','MarkerSize',10,'LineWidth',2,'Color',colors(i,:));
        hold on;

        % Estimation en moindres carrés pondérés :
        X = moindres_carres_ponderes(D_app,probas(i,:));
        parametres_estim(i,:) = conversion(X);
        [x,y] = points_ellipse(parametres_estim(i,:),theta_affichage);
        plot([x x(1)],[y y(1)],'-','LineWidth',3,'Color',colors(i,:));
    end


    % Calcul du nouveau score :
	score_nouv = calcul_score_n(parametres_VT,parametres_estim,n);
	difference_score = score_nouv-score;
	score = score_nouv;

    axis([-taille taille -taille taille]);
    xlim([-taille taille]);
    ylim([-taille taille]);

    F = getframe(gcf);
    [X, Map] = frame2im(F);
    [fimage, map] = rgb2ind(X, 256);
    if k == 0
        imwrite(fimage, map, "saves/exo5.gif", 'gif', 'Loopcount', inf);
    else
        imwrite(fimage, map, "saves/exo5.gif", 'gif', 'WriteMode', 'append');
    end
    
    k = k + 1;
end

% Affichage du score final :
fprintf('Score de l''estimation par EM : %.3f\n',score);