// Données

X = 1:4;
Y = [6, 5, 7, 10];

// régression horizontale
[a, b, sigma] = reglin(Y, X);

// régression verticale
[aa, bb, sigma] = reglin(X, Y);

// régression orthogonale
Ainit = [0, 1];

function [e] = resorth(A, X, Y)
    e = 1/(1 + A(2)^2)*(Y - A(1) - A(2)*X).^2
endfunction

[S, Aopt] = leastsq(list(resorth, X', Y'), Ainit);

// points projetés sur la droite

Xh = a*Y + b;
Yv = aa*X + bb;
normdo = sqrt(1 + Aopt(2)^2);
distcurv = (1/normdo)*(X + Aopt(2)*(Y - Aopt(1)));
Xorth = (1/normdo)*distcurv;
Yorth = Aopt(1) + (Aopt(2)/normdo)*distcurv;

// droite de régression verticale
xx1 = (y1 - bb)/aa;
xx2 = (y2 - bb)/aa;

// droite de régression orthogonale
xxx1 = (y1 - Aopt(1))/Aopt(2);
xxx2 = (y2 - Aopt(1))/Aopt(2);

// tracé
clf;
couleurs = [get(sdf(), 'color_map') ; 0.75, 0.75, 0.75];
xset('colormap', couleurs);

xsegs([X ; X], [Y ; Yv], 14) // segments verticaux
h2 = gce();
h2.thickness = 2;

xsegs([X ; Xorth], [Y ; Yorth], 14) // segments orthogonaux
h2 = gce();
h2.line_style = 2;

xpoly([xx1, xx2], [y1, y2]) // droite de régression verticale
h3 = gce();
h3.thickness = 2;
h3.foreground = 2;

xpoly([xxx1, xxx2], [y1, y2]) // droite de régression orthogonale
h5 = gce();
h5.line_style = 2;
h5.foreground = 2;

plot(X, Y, 'ok') // points
h5 = gce();
h5.children.mark_background = 5;
axe = gca();
axe.data_bounds = [0, 4 ; 5, 10.25];
axe.grid = [33, 33];
axe.tight_limits = 'on';
axe.isoview = 'on';
xtitle(' ', 'x', 'y')