
disp('Testing Parallelism for Several Lines');
disp('(See Section 9.5 in text (attributed to Sen, Adichie))');
disp('See the comments in ParallelSlopes.m and ParallelSlopes.txt');
disp('  (``Approach I``) for a discussion of the method');

disp(' ');
disp('EXAMPLE:');
disp('Test of cloud seeding of cyclones');
disp('  kk=2 samples: Seeded and not seeded (in that order)');
disp('  Y_ij = measure of precipitation');
disp('  X_ij = meridional circulation index for cyclones');
disp('  Sample sizes n1=11 (seeded) and n2=10 (unseeded)');
disp('See Table 9.5 page 436 in text (Problem 28)');

% Data for seeded cyclones
ww1 = [  24   0.180 ;
         28   0.175 ;
         30   0.178 ;
         37   0.021 ;
         43   0.260 ;
         47   0.715 ;
         52   0.441 ;
         57   0.205 ;
         71   0.417 ;
         87   0.498 ;
        115   0.603 ];

% Data for unseeded cyclones
ww2 = [ -7    0.138 ; 
        10    0.081 ; 
        17    0.072 ; 
        25    0.188 ; 
        44    0.075 ; 
        51    0.435 ; 
        53    0.423 ; 
        63    0.339 ; 
        75    0.519 ; 
        90    0.738 ];

kk=size(ww1,2);
nn1=size(ww1,1);  nn2=size(ww2,1);
nnn = [nn1 nn2];
nmax=max(nnn);

fprintf ('\nK=%d groups with sizes N1=%d and N2=%d\n', kk,nn1,nn2);

% Data as column vectors
xx1 = ww1(:,1);  yy1=ww1(:,2);
xx2 = ww2(:,1);  yy2=ww2(:,2);

disp(' ');
disp('Data:');
fprintf(' Seeded (n=%d):\n', nn1);
fprintf('   Y:');
for j=1:nn1;  fprintf(' %5.3f',yy1(j));  end;  fprintf('\n');
fprintf('   X:');
for j=1:nn1;  fprintf(' %3d  ',xx1(j));  end;  fprintf('\n');

fprintf(' Unseeded (n=%d):\n', nn2);
fprintf('   Y:');
for j=1:nn2;  fprintf(' %5.3f',yy2(j));  end;  fprintf('\n');
fprintf('   X:');
for j=1:nn2;  fprintf(' %3d  ',xx2(j));  end;  fprintf('\n');

disp(' ');
disp('The H_0 model is  Y_ij = beta X_ij + mu_i + error_ij.');
disp('  As in ParallelSlopes.m, minimizing over the mu_i for fixed');
disp('  leads to the system');
disp(' ');
disp('     Y_ij - Ybar_i = beta (X_ij - Xbar_i) + error_ij       (4)');
disp(' ');
disp('  Let beta be the overall least-squares estimate of beta in (4)');
disp('  for all (i,j) pairs (1 le i le kk, 1 le j le nn_i)');

% First, define 2xnmax matrices for X and Y
% Values beyond the end of ``short rows'' (nnn(i)<nmax) are zero
xx=zeros(2,nmax);  yy=zeros(2,nmax);
for j=1:nn1;  xx(1,j)=xx1(j);  yy(1,j)=yy1(j);  end;
for j=1:nn2;  xx(2,j)=xx2(j);  yy(2,j)=yy2(j);  end;

% Define the 2xnmax row-centered matrix  Xc_ij = X_ij - Xbar_i 
xxc=zeros(kk,nmax);   % Same shape as xx
for i=1:kk;   % For each row or group
  nn=nnn(i);  % Length of this X row
  % Find the i-th row mean
  xsum=0;  for j=1:nn;  xsum=xsum+xx(i,j);  end;
  xmean=xsum/nn;
  for j=1:nn;  xxc(i,j)=xx(i,j)-xmean;  end;
  end;

% Estimate of overall beta = beta_bar
bsum=0;  xsum=0;  % Numerator and denominator
for i=1:kk;
  nn=nnn(i);
  for j=1:nn;
    bsum=bsum+xxc(i,j)*yy(i,j);
    xsum=xsum+xxc(i,j)*xxc(i,j);
  end;  end;
betabar = bsum/xsum;
fprintf('\nOverall beta (beta_bar) = %.6f\n',betabar);
% For comparison, compute the within-group betas as well
% These are not used in the test procedure
% Note that values beyond the end of ``short rows'' are not used.
bbw=zeros(kk,1);
for i=1:kk;
  bsum=0;  xsum=0;  nn=nnn(i);
  for j=1:nn;
    bsum=bsum+xxc(i,j)*yy(i,j);
    xsum=xsum+xxc(i,j)*xxc(i,j);
    end;
  bbw(i)=bsum/xsum;  end;
fprintf('Beta values within groups:  %.6f  %.6f\n', bbw);

disp(' ');
disp('Aligned observations  Ystar_ij = Y_ij - betabar X_ij:');
yystar = yy - betabar*xx;
fprintf(' Seeded (n=%d):\n', nn1);
fprintf('   ');
for j=1:nn1;  fprintf(' %5.3f',yystar(1,j));  end;  fprintf('\n');
fprintf(' Unseeded (n=%d):\n', nn2);
fprintf('   ');
for j=1:nn2;  fprintf(' %5.3f',yystar(2,j));  end;  fprintf('\n');

% Find within-group ranks for yystar
yvals=zeros(nmax,1);  % Column-vector buffer for values from one row
rr=zeros(kk,nmax);    % For Friedman ranks: Same shape as xx and yy 
for i=1:kk;
  nn=nnn(i);
  % Copy a row (group) of yystar to column-vector yvals
  for j=1:nn;  yvals(j)=yystar(i,j);  end;
  % Find within-group midranks of column-vector yvals
  [midr tiegr] = getmidranks(yvals);
  % Copy midranks to the i-th row of rr
  for j=1:nn;  rr(i,j)=midr(j);  end;
  end;
% Within-group ranks are listed in the next table

% Under permutation, ranks are approximately uniformly distributed
%   with range (1,nn), mean (nn+1)/2, and variance approx nn^2/12
% Define two column vectors for row (within-group) values
tt=zeros(kk,1);  ci2=zeros(kk,1);
for i=1:kk;
    tsum=0;  csum=0;
    nn=nnn(i);
    for j=1:nn;
      tsum=tsum+xxc(i,j)*rr(i,j);
      csum=csum+xxc(i,j)*xxc(i,j);
      end;
    ti(i) = tsum/(nn+1);
    ci2(i) = csum;
    end;
disp(' ');
disp('Within-group ranks:');
fprintf(' Seeded:   ');
for j=1:nn1;  fprintf(' %4.1f', rr(1,j));  end;  fprintf('\n');
fprintf(' Unseeded: ');
for j=1:nn2;  fprintf(' %4.1f', rr(2,j));  end;  fprintf('\n');

disp(' ');
disp('Within-group Ti=Sum_j Xc(i,j)*Rank(i,j)/(n_i+1)  and');
disp('   Ci2=Sum(j=1,nn_i) (X_ij - Xbar_i)^2:');
fprintf(' Seeded:     T1=%6.2f   C1^2=%8.3f\n', ti(1), ci2(1));
fprintf(' Unseeded:   T2=%6.2f   C2^2=%8.3f\n', ti(2), ci2(2));

% Sen-Adichie`s test:
% Will have a chi-square distribution given H_0
vscore=0;
for i=1:kk;  vscore=vscore+ti(i)*ti(i)/ci2(i);  end;
vv = 12*vscore;
df=kk-1;
fprintf('\nSen-Adichie test:  X=%6.4f   P=%6.4f  (df=%d)\n',...
  vv, pvchisq(vv,df), df);

disp(' ')
disp('Compare with a permutation test of within-group ranks:');
disp('(Friedman-like permutations.)');
disp(' ');
disp('Initializing the random-number generator at');
Startseed=12345678;
fprintf('    %d\n',Startseed);
rand('state',Startseed);

nsims=10000;
fprintf('\nCarrying out ns=%d sets of within-group permutations:\n',nsims);

Ngr=0;
% Define two buffers
rstar=zeros(kk,nmax);  rvals=zeros(1,nmax);
% vscore as defined above does not simplify, so use it
%   as the observed statistic

fprintf('Observed Vscore = %g\n',vscore);
for ns=1:nsims;
    %% Randomly permute the two rows
    rstar(1,1:nn1)=rr(1,randperm(nn1));
    rstar(2,1:nn2)=rr(2,randperm(nn2));
    %% If running time is a consideration, the following
    %%   should be replaced by faster MATLAB code using sum()
    tt1=0;
    for j=1:nn1;  tt1=tt1+xxc(1,j)*rstar(1,j);  end;
    tt1=tt1/(nn1+1);
    tt2=0;
    for j=1:nn2;  tt2=tt2+xxc(2,j)*rstar(2,j);  end;
    tt2=tt2/(nn2+1);
    ttstar=tt1*tt1/ci2(1) + tt2*tt2/ci2(2);
    if ttstar>=vscore Ngr=Ngr+1;  end;
    end;

disp('Number of simulations with values >= Tscore and total number:');
disp(' ');
disp([Ngr nsims]);

phat=Ngr/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% CI for true P-value bracketing estimate of true pvalue:');
disp('(Since H and Fscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);

disp(' ');

