
disp('One-Way Layouts:');
disp('Kruskal-Wallis test for equality of location parameters');
disp('  Model:  X(i,j) dist X + theta_i ---  Test H_0: Theta are identical');

disp('  ')
disp('Example 2:');
disp('Testing smoothness of paper in four batches');

disp(' ');

x1 = [ 387  415  438  445  455  460  477  580 ];
x2 = [ 392  393  397  414  418  429  433  458 ];
x3 = [ 340  350  390  400  430  430  440  450 ];
x4 = [ 340  348  348  354  372  378  412  428 ];

disp('Data in 4 treatment groups with sizes 8,8,8,8:');

n1=length(x1);  n2=length(x2);  n3=length(x3);  n4=length(x4);
kk=4;  nn=n1+n2+n3+n4;

disp(' ');
fprintf('  Group #1 (n=%d)',n1);
for j=1:n1;  fprintf (' %4d', x1(j));  end;  fprintf('\n');
fprintf('  Group #2 (n=%d)',n2);
for j=1:n2;  fprintf (' %4d', x2(j));  end;  fprintf('\n');
fprintf('  Group #3 (n=%d)',n3);
for j=1:n3;  fprintf (' %4d', x3(j));  end;  fprintf('\n');
fprintf('  Group #4 (n=%d)',n4);
for j=1:n4;  fprintf (' %4d', x4(j));  end;  fprintf('\n');

disp(' ');
disp('In one vertical table for analysis:');
disp('(Columns are Ordinal, Value, Treatment group, Midrank, Tiegroups');

dd = [ x1' ones(n1,1);  x2' 2*ones(n2,1);  x3' 3*ones(n3,1);  x4' 4*ones(n4,1)];  %

[midr tiegr] = getmidranks(dd);    % Midranks of 1st column of dd;

dd = [ [1:nn]'  dd  midr  tiegr];  %'  dd is now nnx5

fprintf('   %2d  %4d   %d  %5.1f  %2d\n', dd');  %'

disp(' ');
disp('Parallel display with midranks with treatment groups as columns:');
r1=midr(1:n1,1);  r2=midr(n1+1:n1+n2,1);
r3=midr(n1+n2+1:n1+n2+n3,1);  r4=midr(n1+n2+n3+1:nn,1);
ddx = [x1'  r1  x2'  r2  x3'  r3  x4'  r4];
fprintf(' %4d (%4.1f)   %4d (%4.1f)   %4d (%4.1f)   %4d (%4.1f)\n', ddx');  %'

disp(' ');
% Distinguished columns of the data dd
vals=dd(:,2);  group=dd(:,3);  midr=dd(:,4);

% Compute Kruskal-Wallis H statistic
% First, accumulate rank sums and group numbers by row
ranksum=zeros(kk,1);  num=zeros(kk,1);
for i=1:nn;  ranksum(group(i)) = ranksum(group(i)) + midr(i);
    num(group(i)) = num(group(i)) + 1;  end;
% ranksum and num are vectors, so we can divide componentwise
rankav=ranksum./num;

disp(' ');
disp('By Treatment Group: Number, Rank Sum, Rank Average:');
rdd = [ [1:kk]'  num  ranksum  rankav];  %'
fprintf('  Group #%d  n=%d   %5.1f  %6.2f\n', rdd');  %'


% Compute H statistic
hh=0;
ravc = rankav - (nn+1)/2;
for k=1:kk;  hh = hh + num(k)*ravc(k)*ravc(k);  end;
hh = (12*hh)/(nn*(nn+1));
% Alternatively:  hh = (12/(nn*(nn+1))*sum(num.*ravc.*ravc);

fprintf('\nKruskal-Wallis H = %g (with no tie correction)\n', hh);

% Adjust for tie correction
Tiesum=0;
for i=1:nn;  tt=tiegr(i);
  if tt>0  ntg=i;  Tiesum = Tiesum + (tt-1)*tt*(tt+1);
  else break;  end;  end;

disp(' ');
disp('Tie-group sizes:');
disp(tiegr(1:ntg)');  %'

Tiecorr = Tiesum/((nn-1)*nn*(nn+1));
disp('Tiesum and Tie correction:');
disp(' ');
disp([Tiesum  Tiecorr]);

if Tiesum>0
  hh = hh/sqrt(1-Tiecorr);  end;

fprintf('Kruskal-Wallis H = %g (with tie correction)\n', hh);

disp(' ');
disp('Large-sample chi-square approximation for H:');

df=kk-1;
Pval = pvchisq(hh,df);
fprintf('\n  P= %6.4f  (df=%d)\n', Pval,df);

disp(' ');
disp('Doing nsims simulations for the true P-value');
nsims=10000

% Initialize random-number generator from the system clock
disp('Initializing the RNG from the system clock at');

startseed=sum(100*clock)

rand('state',startseed);

% Note that H = A_n * HScore + C_n   for  A_N=12/(nn*(nn+1)), C_N=-3(nn+1)
%   for  Hscore = Sum ranksum^2/num
% Thus simulating  Hscore  is the same as simulating  H (hh)

% Observed Hscore:
Hscore=0;
for k=1:kk;  Hscore=Hscore + ranksum(k)*ranksum(k)/num(k);  end;
 
% We could also say:  Hscore=sum(ranksum.*ranksum./num);

fprintf('H=%g  and  Hscore=%g\n', hh,Hscore);

% Do nsims permutations for a test of Hscore
Ngr=0;  group=dd(:,3);  midr=dd(:,4);
for ns=1:nsims;
  %% Permuted copy of midranks:
  midrperm=midr(randperm(nn));
  %% Store rank sums for permuted midranks in rrsum
  rrsum=zeros(kk,1);
  for i=1:nn;  rrsum(group(i))=rrsum(group(i))+midrperm(i);  end;
  %% Calculate the hscore for the permuted midranks
  Hrand = sum(rrsum.*rrsum./num);
  if Hrand >= Hscore  Ngr=Ngr+1;  end;
  end;

disp(' ');
disp('Number of simulations with values >= Hscore 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 Hscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);
  

disp(' ');
disp('Test for an alternative of monotonic treatment-group medians:');
disp('The Jonckheere-Terpstra test statistic is very similar to');
disp('  the Kendall nonparametric correlation statistic between');
disp('  values and treatment-group numbers');

val=dd(:,2);  group=dd(:,3);
% Calculate J statistic
jj=0;
for i=1:nn;  for j=1:nn;
  g1=group(i);  g2=group(j);
  if g1<g2  val1=val(i);  val2=val(j);
    if val1<val2 jj=jj+1;
    elseif val1==val2 jj=jj+0.50;  end;  end;  % For nested ifs
  end;  end;  % For double loop

%% Display the statistic J
jj

disp('The Jonckheere-Terpstra test statistic mean and variance:');
disp('(WITHOUT tie correction.)');

jmean = (nn*nn - sum(num.*num))/4;
jvar = (nn*nn*(2*nn+3) - sum(num.*num.*(2*num+3)))/72;

zz = (jj-jmean)/sqrt(jvar);

Jmean_Jvar_Zval_NormPval = [jmean  jvar  zz  normpval(zz)/2]




