x=[4 7 8 9 12 18 19 21 25 30];
alpha = 0.05
 n=length(x);
%n=100
xx=sort(x)
y=normcdf(xx, 10, 5);  
%y=chi2cdf(x', 10);     %10 mean

    i1n = ( 2*(1:n) -1 )/(2 * n);
    u2 = sum( (y - i1n - mean(y) + 1/2).^2 )+ 1/(12 * n)    %0.1134
           a= (3*n - 8)/(60*(2*n - 3))
           b=(2*n - 3)/(84*n)
           df = 49 * n * (n-1)/(5*(2*n -3)^2)
    % 
     u2crit = a + b* chi2inv(1-alpha, df)    % 0.1875
     chi2 = (u2-a)/b                                       %3.7413
     pval = 1-chi2cdf(chi2, df)                    % 0.2152
    
nomega2 = 1/(12 * n) + sum( (normcdf(xx,10,5) - i1n).^2 ) 

u2
nomega2 - n*(mean(y) - 1/2)^2
%%
clear all
close all
alpha=0.05

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 ];

%
rs=ranks([x y]);
m=length(x); n=length(y);
r=rs(1:m); s=rs(m+1 : m+n);
Wmn = 1/(m*n*(m+n)) * ...
    (  m * sum( (sort(r)-(1:m)).^2 ) + n * sum( (sort(s)-(1:n)).^2 ) ) - ...
    (4 * m * n - 1)/(6 * (m+n))      %0.6715
pval = 1-cvmcdf(Wmn)                 %0.0149


%%
close all
clear all
pvals=[]; Ws=[];
B=10000;
m=9; n=15;
for i = 1:B
  x=sort(normrnd(0,1,[1 m]));
  y=sort(normrnd(0,1,[1 n]));
    
rs=ranks([x y]);
m=length(x); n=length(y);
r=rs(1:m); s=rs(m+1 : m+n);
Wmn = 1/(m*n*(m+n)) * ...
    (  m * sum( (r-(1:m)).^2 ) + n * sum( (s-(1:n)).^2 ) ) - ...
    (4 * m * n - 1)/(6 * (m+n));
Ws=[Ws Wmn];
pvals =[pvals 1-cvmcdf(Wmn)];
end
hist(pvals)

    
