
disp('Two-Way Layout:');
disp('Data  X(i,j)  where');
disp('   1 le i le nn  Subjects or Blocks');
disp('   1 le j le kk  Treatments or Treatment Groups');
disp('Model:  X(i,j) = mu(i,j) plus error  or mu + theta_i + tau_j plus error');
disp('Test H_0:  mu(i,j) = theta_i  and there is no effect of Treatment');

disp(' ');
disp('Example 1:');
disp('Methods for rounding first base in baseball');
disp('k=3 methods with average times to second base for n=22 baseball players');
disp('   Method 1: Round out near first base');
disp('   Method 2: Narrow angle at first base');
disp('   Method 3: Wide angle at first base');
disp('Table 7.1 (p274) and Figure 7.1 (p275) in text');

xxr = [ 1   5.40   5.50   5.55 ; 
        2   5.85   5.70   5.75 ; 
        3   5.20   5.60   5.50 ; 
        4   5.55   5.50   5.40 ; 
        5   5.90   5.85   5.70 ; 
        6   5.45   5.55   5.60 ; 
        7   5.40   5.40   5.35 ; 
        8   5.45   5.50   5.35 ; 
        9   5.25   5.15   5.00 ; 
       10   5.85   5.80   5.70 ; 
       11   5.25   5.20   5.10 ; 
       12   5.65   5.55   5.45 ; 
       13   5.60   5.35   5.45 ; 
       14   5.05   5.00   4.95 ; 
       15   5.50   5.50   5.40 ; 
       16   5.45   5.55   5.50 ; 
       17   5.55   5.55   5.35 ; 
       18   5.45   5.50   5.55 ; 
       19   5.50   5.45   5.25 ; 
       20   5.65   5.60   5.40 ; 
       21   5.70   5.65   5.55 ; 
       22   6.30   6.30   6.25  ];

[nnr kkr] = size(xxr);

xx = xxr(:,2:kkr);       % Data without row numbers

[nn kk] = size(xx);

fprintf('\nn=%d blocks (subjects) and k=%d treatment groups\n',nn,kk);

disp(' ');
disp('Player, Round Out, Narrow Angle, Wide Angle');
%%disp(xx);
fprintf('   %2d.   %7.2f   %7.2f   %7.2f\n',xxr');  %'

fprintf('\nFRIEDMAN RANKS are WITHIN-BLOCK ranks, here 1 le R(i,j) le %d\n',kk);
fprintf('Times for n=%d players with within-subject ranks/midranks\n',nn);

%% Ranks and tie-groups for 22 blocks

rr=zeros(nn,kk);  tg=zeros(nn,kk);
sbuff=zeros(kk,1);  % Column-vector buffer
for i=1:nn;  % For all replications
  % Put values for that block into sbuff
  for j=1:kk;  sbuff(j)=xx(i,j);  end;
  [midr  tiegr] = getmidranks(sbuff);
  for j=1:kk;  rr(i,j)=midr(j);  tg(i,j)=tiegr(j);  end;
  end;

% Buffer for a table with row numbers then kk sets of data then ranks
xxd = zeros(nn,2*kk+1);

xxd(:,1) = [1:nn]';  %'
for j=1:kk;  xxd(:,2*j)=xx(:,j);  xxd(:,2*j+1)=rr(:,j);  end;

%%disp(xxd);
fprintf('   %2d.   %5.2f(%3.1f)   %5.2f(%3.1f)   %5.2f(%3.1f)\n',xxd');  %'

disp(' ');
disp('Rank sums and rank averages for each treatment group:');
rsum=sum(rr,1);   % 1xkk row vector of column sums of nnxkk matrix rr
ravg=rsum/nn;     % Also a 1xkk row vector

fprintf('   %2d. %7.2f   %6.2f\n',[ [1:kk]'  rsum'  ravg']');

%% Friedman statistic
%% Centered rank averages
ravgc = ravg - (kk+1)/2;

ss = (12*nn/(kk*(kk+1)))*sum(ravgc.*ravgc);

%Alternatively:
%ss=0;
%for j=1:kk;  ss=ss+ravgc(j)*ravgc(j);  end;
%ss = (12*nn/(kk*(kk+1)))*ss;

df=kk-1;
disp(' ');
disp('With no tie correction:');
Pval = pvchisq(ss,df);
fprintf('  S=%7.5f   P= %6.4f  (df=%d)\n', ss,Pval,df);

% Compute tie correction
tiesum=0;
for i=1:nn;  tgg=0;
  for j=1:kk;  tt=tg(i,j);
    if tt>1 tgg=tgg+(tt-1)*tt*(tt+1);  end;  end;
  tiesum=tiesum+tgg;  end;
tiecorr = tiesum/(nn*(kk-1)*kk*(kk+1));
ssp = ss/(1-tiecorr);
disp(' ');
disp('With tie correction:');
fprintf('  Tiesum=%g   Tiecorr=%g\n', tiesum,tiecorr);
Pval = pvchisq(ssp,df);
fprintf('  SP=%7.5f   P= %6.4f  (df=%d)\n', ssp,Pval,df);

disp(' ');
disp('Simulating P-values for the Friedman test:');
disp('Permuting method values independently for each subject');

disp(' ');
disp('Initializing the random-number generator at');
Startseed=12345678;
disp(Startseed);
rand('state',Startseed);

nsims=10000;

fprintf('Carrying out Nsims=%d Friedman permutations:\n',nsims);

% Note that S = A_n * FScore + C_n   for  A_N=12n/(k*(k+1)), C_N=-3n(k+1)
%   for  FScore = Sum(j=1,k) ranksum^2
% Thus simulating  S  is the same as simulating  FScore

%Fscore=0;
%for j=1:kk;  Fscore=Fscore+rsum(j)*rsum(j);  end;
Fscore = sum(rsum.*rsum);
disp('Fscore = Sum of Squares over Treatment Groups');
fprintf('  of Treatment-Group Rank Sums: %g\n',Fscore);
 
% Do nsims permutations for a test of Fscore
rrp=rr;   % Preserve original ranks
Ngr=0;
for ns=1:nsims;
  %% Permute rows of rrp independently
  for i=1:nn;  rrp(i,1:kk) = rrp(i,randperm(kk));  end;
  %% For a matrix rrp, sum(rrp,1) is the row vector of column sums
  %% Sum of squares of the sum of the columns of rrp
  Rscore=sum(sum(rrp,1).*sum(rrp,1));
  if Rscore >= Fscore  Ngr=Ngr+1;  end;
  end;

disp('Number of simulations with values >= Fscore 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]);




