%stars = load('stars.txt'); 
  n = size(stars,1);
   X = [ones(n,1) stars(:,2)]; y = stars(:,3);
   bols  = X\y; [ignore,idx] = sort(stars(:,2));
   plot(stars(:,2),stars(:,3),'o',stars(idx,2),...
      X(idx,:)*bols,'-.' )
   legend('Data','OLS')
  %
   %   Least Absolute Deviation
   blad = medianregress(stars(:,2),stars(:,3));
   plot(stars(:,2),stars(:,3),'o',stars(idx,2),...
     X(idx,:)*bols,'-.',stars(idx,2),X(idx,:)*blad,...
     '-.') 
 legend('Data','OLS','LAD');
   %
   % Huber Estimation
   k = 1.345; % tuning parameters in Huber's weight function
   wgtfun  = @(e) (k*(abs(e)>k)-abs(e).*(abs(e)>k))./abs(e)+1;
   %  Huber's weight function
   wgt = rand(length(y),1);   % Initial Weights
  b0 = lscov(X,y,wgt);
   res = y - X*b0;            % Raw Residuals
  res = res/mad(res)/0.6745; % Standardized Residuals
   m = 30; 
   for i=1:m
       wgt = wgtfun(res);
       % Compute the weighted estimate using these weights
       bhuber = lscov(X,y,wgt);
       if all((bhuber-b0)<.01*b0)% Stop with convergence
         return;
       else
         res = y - X*bhuber;
         res = res/mad(res)/0.6745;
      end
    end
   plot(stars(:,2),stars(:,3),'o',stars(idx,2),X(idx,:)...
      *bols,'-.', stars(idx,2),X(idx,:)*blad,'-x',...
       stars(idx,2),X(idx,:)*bhuber,'-s')
      legend('Data','OLS','LAD','Huber');
  %
   % Least Trimmed Squares
   blts = lts(stars(:,2),y);
  plot(stars(:,2),stars(:,3),'o',stars(idx,2),X(idx,:)...
     *bols,'-.',stars(idx,2),X(idx,:)*blad,'-x',...
     stars(idx,2),X(idx,:)*bhuber,'-s', stars(idx,2),...
     X(idx,:)*blad,'-+')
     legend('Data','OLS','LAD','Huber','LTS');
   %
   % Least Median Squares
   [LMSout,blms,Rsq]=LMSreg(y,stars(:,2));
   plot(stars(:,2),stars(:,3),'o',stars(idx,2),X(idx,:)...
     *bols,'-.', stars(idx,2),X(idx,:)*blad,'-x',...
     stars(idx,2), X(idx,:)*bhuber,'-s',stars(idx,2),...
     X(idx,:)*blad,'-+', stars(idx,2),X(idx,:)*blms,'-d')
     legend('Data','OLS','LAD','Huber','LTS','LMS');