%
clear all
close all  force
  disp('Comparing Smirnov and Rosenblatt''s Tests')
  set(0, 'DefaultAxesFontSize', 16);
  lw =3;  fs = 16;  msize = 10;
s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);

m = 39; n= 35;

bmn = (m+n)/(2*m*n) * (min(m,n) - gcd(m, n))/(m + n + gcd(m,n)) ;
EWmn = 1/6 * (1 + 1/(m+n));
VarWmn = 1/45 * (1+1/(m+n) )*(1+ 1/(m+n) - ...
                 3/4 * (1/m + 1/n));

%
i = 1:m; j = 1:n;
pvals=[]; 
pvalr=[];
B=20000;
for b=1:B
x = sort(normrnd(-1 , 2, [1,  m]));
y = sort(trnd(10, [1, n]));

%SMIRNOV TEST
distances  =  [i/m - empiricalcdf(x, y)  j/n - empiricalcdf(y, x)];
Dmn =  max(distances) ;    
bmn = (m+n)/(2*m*n) * (min(m,n) - gcd(m, n))/(m + n + gcd(m,n));
zz =  sqrt(m * n/(m + n))*(Dmn + abs(m-n)/(6*m*n) +  bmn );
pvals  = [pvals  1 - kscdf(zz)];

%ROSENBLATT TEST
z = [x y]; 
rz = ranks(z);
Ri = rz(1:m);  Sj=rz(m+1:m+n);
Wmn = 1/(m+n)*(1/6 + 1/n * sum( (Ri - (1:m)).^2) + ...
      1/m * sum( (Sj - (1:n) ).^2 )) - 2 * m * n/(3*(m+n));
Wprimemn = 1/6 + (Wmn-EWmn)/sqrt(45*VarWmn);
pvalr = [pvalr 1-cvmcdf(Wprimemn)];

end
% 


sum(pvals < 0.05)/B
sum(pvalr < 0.05)/B
figure(1);
  hist(pvals, 100)
axis([0 0.1 0 15000 ])
 figure(2) 
hist(pvalr, 100)
axis([0 0.1 0 15000 ])


 
