%
%Illustration of Rosenblatt's Test
clear all
close all
 disp('Rosenblatt Test Simulated Example')
 lw = 2; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 16;
 msize = 8;
 
 %Note that both samples should be sorted
x =[...
   -5.75    -3.89    -3.69    -3.68    -3.54    -2.59 ...
   -2.53    -2.40    -2.39    -2.27    -2.02    -1.72 ...
   -1.64    -1.54    -1.27    -1.11    -1.08    -1.00 ...
   -0.88    -0.84    -0.47    -0.36    -0.29    -0.24 ...
   -0.20    -0.19     0.02     0.15     0.25     0.45 ...
    0.51     0.72     0.74     0.96     1.25     1.33 ...
    1.49     2.39     2.59 ];

y=[...
   -2.72    -2.18    -1.31    -1.17    -1.00    -0.94 ...
   -0.78    -0.65    -0.63    -0.52    -0.40    -0.37 ...
   -0.30    -0.21    -0.19    -0.11    -0.05     0.12 ...
    0.14     0.25     0.27     0.35     0.45     0.48 ...
    0.48     0.60     0.71     0.76     0.79     1.01 ...
    1.10     1.10     1.12     1.36     2.03 ];
 %%Note that both samples should be sorted, if not
 % x=sort(x); y=sort(y);
 
m = length(x); n=length(y);

m = length(x); n=length(y);
i = 1:m; j = 1:n;
z = [x y]; %combined
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)) % 0.6715
 pval1 = 1-cvmcdf(Wmn) %0.0149
 

% Rosenblatt's test with corrections
EWmn = 1/6 * (1 + 1/(m+n))  %0.1689
VarWmn = 1/45 * (1+1/(m+n) )*(1+ 1/(m+n) - ...
                 3/4 * (1/m + 1/n))   %0.0219

Wprimemn = 1/6 + (Wmn-EWmn)/sqrt(45*VarWmn)  %0.6728
 pval2 = 1-cvmcdf(Wprimemn)   %0.0148

% 
close all
xx = -7:0.01:5;
 
plot(xx, empiricalcdf(xx, x), 'r-', 'linewidth',2)
hold on
plot(xx, empiricalcdf(xx, y), 'k-','linewidth',2)
axis tight
legend('emp cdf from N(-1,2^2)','emp cdf from t_{10}','Location','NW')
print -depsc 'C:\STAT\Goodness\Goodnesseps\smirnovnt.eps'
