
disp('Rank Regression for a Multiple Regression');
disp('  Y_i = mu + Sum(j=1,p) X_{ij} beta_j + error_i');
disp('Specialized to p=1, so   Y_i = beta*X_i + mu');
disp('See Sections 9.6-9.7 in text');

disp(' ');
disp('Data: Amount of pre-packed engine lubricant and engine survival times:');
disp('(Source of data unknown.)');
disp(' ');
disp('Lubricant  Survival          Lubricant Survival');

ww2 = [ 1.   29     23     11.   84    429;
        2.   21     36     12.   37    489;
        3.   30     67     13.   94    504;
        4.   45    104     14.   56    925;
        5.   48    147     15.   67    980;
        6.   92    164     16.   79   1254;
        7.   74    247     17.   70   2961;
        8.   61    304     18.   93   5351;
        9.   99    355     19.   76   9915;
       10.   24    408     20.   72  11301  ];

fprintf (' %3d.  %2d  %5d          %3d.   %2d  %5d\n', ww2');  %'

disp(' ');
disp('Let X=Lubricant and Y=Survival (in days)');

% Rewrite data as one nnx3 array
ww = [ ww2(:,1) ww2(:,2) ww2(:,3);  ww2(:,4) ww2(:,5) ww2(:,6)];
nn=size(ww,1);      % Number of rows

%% Sorting data into one nnx3 matrix
ww=sortrows(ww,1);  % Sort on ordinal number

% Fetch Xs and Ys as nnx1 column vectors
xx=ww(:,2);  yy=ww(:,3);

% Are there any X- or Y-ties?
[midrx  xties] = getmidranks(xx);
xg=0;
for i=1:nn;  
  if xties(i)>0 xg=i;  else break;  end;  end;
[midry  yties] = getmidranks(yy);
yg=0;
for i=1:nn;  
  if yties(i)>0 yg=i;  else break;  end;  end;
fprintf('\nWithin-sample tie groups:  X: %d   Y: %d\n', xg, yg);
if xg==0 && yg==0  disp('NO TIES in Xs or Ys');  end;

%% Find slope and intercept of simple regression
%% (After we find the slope, the intercept will be the median
%%   of the Walsh averages of  Y_i - Sum(j=1,p) X_{ij} beta_j,
%%   so that finding the slope is the hard part.)

% Collect slopes and absolute-X-differences in an mmx2 matrix
%   for mm=nn*(nn-1)/2
mm=nn*(nn-1)/2;
ww=zeros(mm,2);  nw=0;
for i=1:nn;  for j=i+1:nn;
  xdiff=xx(j)-xx(i);
  if xdiff~=0  nw=nw+1;
    ww(nw,1)=(yy(j)-yy(i))/xdiff;
    ww(nw,2)=abs(xdiff);  end;
  end;  end;

ww = ww(1:nw,:);    %% Trim off any unused rows
% Sort ww by slope then X-difference
ww = sortrows(ww, [1 2]);

disp(' ');
disp('Finding the minimum value of D(beta):');
fprintf('  nw=%d nodes ranging from  %g  to  %g\n', nw,ww(1,1),ww(nw,1)); 

% Compute Qsum = Sum(i=1,nn) (Rank(X_i) - (nn+1)/2)*X_i
%   =  Sum(i=1,nn) (Rank(X_i) - (nn+1)/2)*(X_i - Xbar)
% Note that midranks of X_i are already in the array midrx
midrxc = midrx - (nn+1)/2;
qsum=0;
for i=1:nn;  qsum=qsum+midrxc(i)*xx(i);  end;
fprintf('  Slope for beta<%g: %g   for beta>%g: %g\n',...
  ww(1,1),-qsum, ww(nw,1),qsum);

disp('Finding the lowest point on the curve Y = D(beta):');
% lslope is initially the slope just before ww(1)
% The beta value at ww(i) is ww(i,1)
lslope=-qsum;
for i=1:nw;
  rslope=lslope+ww(i,2);  % Add |X_{j_2} - X_{j_1}|
  %% Define beta and end loop if the slope changes sign
  %%   or has a flat spot
  if abs(lslope)<1e-12
     beta=(ww(i,1)+ww(i+1,1))/2;  break;
  elseif lslope<0 && rslope>0
     beta=ww(i,1);  break;
  end;
  %% The options above set beta and exit the loop
  %% Otherwise update lslope and go to next value
  lslope=rslope;
  end;
fprintf('  At minimum value, lslope=%g  rslope=%g     Qsum=%g\n', ...
  lslope, rslope, qsum);

%% Also find D(beta)
resid=zeros(nn,1);
for i=1:nn;  resid(i)=yy(i)-beta*xx(i);  end;
[midr tiegr] = getmidranks(resid);
midrc = midr - (nn+1)/2;
dbeta = sum(midrc.*resid);
fprintf('  Minimum attained at  beta=%g   D(beta)=%g\n', beta, dbeta);

%% Set mu = median of Walsh averages of Y-beta*X values
wav=zeros(nn*nn,1);  nw=0;
for i=1:nn;  for j=i:nn;
  nw=nw+1;  wav(nw) = (resid(i)+resid(j))/2;  end;  end;
wav=wav(1:nw);    %% Trim off unused values
mu=median(wav);
fprintf ('  Intercept:  mu=%g  (median of Y_i-beta X_i Walsh averages)\n', mu);

%% Find average-absolute and RMS errors
resid=yy-beta*xx-mu;
err_avabs = sum(abs(resid))/nn;
err_rms = sqrt(sum(resid.*resid)/nn);

disp(' ');
disp('Regression line and errors:');
disp('  E1 is mean of absolute values of residuals');
disp('  E2 is root-mean-square of residuals');
fprintf('  Rank-Regression:  Y = %7.4f X + %8.3f   E1=%7.2f   E2=%7.2f\n', ...
   beta, mu, err_avabs, err_rms);

disp(' ');
disp('See TheilRegression.txt for the corresponding results');
disp('  for Theil`s and least-squares regressions.');
disp(' ');

