/* Watson U^2-statistic test modified for normality */
/* This program is a GAUSS version of Watson U^2-statistic in 'gofnorm' */
/* (Fortran code by Paul Johnson) located in library 'fBasics' of R.    */
new; cls;
x=rndn(100,1);
call WatsonModified(x);



proc WatsonModified(x);
   local n,fn1,fn2,fx,a,b,sum3,sum4,sum5,zbar,cvm,W;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   x=sortc(x,1); 
   x=(x-meanc(x))/stdc(x);
   fn1=seqa(1,1,n)/n;
   fn2=(2*seqa(1,1,n)-1)/(2*n);
   fx=cdfn(x);
   a=(2*seqa(1,1,n)-1).*ln(fx);
   b=(2*seqa(1,1,n)-1).*ln(1-rev(fx));
   sum3=sumc(a+b);
   sum4=sumc((-fn2+fx)^2);
   fn1=abs(fn1-fx);
   cvm=1/(12*n)+sum4;
   sum5=sumc(fx);
   zbar=sum5/n;
   W=cvm-n*(zbar-0.5)*(zbar-0.5);
   W=W*(1+0.5/n);
   print "U^2=" W;
   retp(W);
endp;