%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   M. Welling %%%%%%% March 2000 %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Linear Regression %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
colordef black;

N_points_1 = input('Number of points you may want to select for the cluster (max 1000)= ');
xs_1=zeros(1,N_points_1);
ys_1=zeros(1,N_points_1);

figure(1);clf;
axis([-100 100 -100 100]);
hold on; 

%%%%%%%% acquire points %%%%%%%%%%%%%%%%%%%
for i=1:N_points_1,

    fprintf(1,'%d - Select a point with the mouse: ',i);
    [x,y] = ginput2(1);
    fprintf(1,'(%3.3f,%3.3f)\n',x,y);
    plot(x,y,'r+');
    xs(i)=x;
    ys(i)=y;

end;

X = [xs;ones(1,N_points_1)]';
Y = ys';

th =  inv(X'*X)*X'*Y;

sg2 = (Y-X*th)'*(Y-X*th)/N_points_1;

t = -100:0.1:100;
f = th(2) + th(1)*t;

plot(t,f,'g','linewidth',1);
plot(t,f+2*sqrt(sg2),'m','linewidth',1);
plot(t,f-2*sqrt(sg2),'m','linewidth',1);