donnees;

n_tests = 500;

% Estimation de la droite de regression par le maximum de vraisemblance :
[theta_Dorth_1,rho_Dorth_1] = estimation_3(x_donnees_bruitees,y_donnees_bruitees,n_tests);

% Affichage de la droite de regression estimee par le maximum de vraisemblance :
cos_theta_Dorth_1 = cos(theta_Dorth_1);
sin_theta_Dorth_1 = sin(theta_Dorth_1);
if abs(cos_theta_Dorth_1)<abs(sin_theta_Dorth_1)
	x_Dorth_1 = x_D;
	y_Dorth_1 = (rho_Dorth_1-x_Dorth_1*cos_theta_Dorth_1)/sin_theta_Dorth_1;
else
	y_Dorth_1 = y_D;
	x_Dorth_1 = (rho_Dorth_1-y_Dorth_1*sin_theta_Dorth_1)/cos_theta_Dorth_1;
end
plot(x_Dorth_1,y_Dorth_1,'b','LineWidth',3);
axis(bornes);
lg = legend('~Droite', ...
	'~Donnees bruitees', ...
	'~$D_\perp$ (maximum de vraisemblance)', ...
	'Location','Best');
set(lg,'Interpreter','Latex');

% Calcul et affichage de l'ecart angulaire :
EA_Dorth_1 = abs(theta_Dorth_1-theta_D);
fprintf('D_perp (maximum de vraisemblance) : ecart angulaire = %.2f degres\n',EA_Dorth_1/pi*180);

function [theta_Dorth_1,rho_Dorth_1] = estimation_3(x_donnees_bruitees,y_donnees_bruitees,n_tests)

    G = mean( [ x_donnees_bruitees.' y_donnees_bruitees.' ] );
    
    x_centre = x_donnees_bruitees - G(1);
    y_centre = y_donnees_bruitees - G(2);
        
    x = repmat(x_centre.', 1, n_tests);
    y = repmat(y_centre.', 1, n_tests);
    
    theta = rand(1, n_tests)*pi;
    
    dist = ( bsxfun(@times, x, cos(theta)) + bsxfun(@times, y, sin(theta)) ).^2;
    
    [min_val, min_index] = min( sum( dist ) );

    theta_Dorth_1 = theta(min_index);
    rho_Dorth_1 = G(1)*cos(theta_Dorth_1) + G(2)*sin(theta_Dorth_1);
    
    if rho_Dorth_1 < 0
       theta_Dorth_1 = theta_Dorth_1 - pi;
       rho_Dorth_1 =  G(1)*cos(theta_Dorth_1) + G(2)*sin(theta_Dorth_1);
    end
    
end