/* Anderson-Darling statistic with critical value by simulation */          
new; cls;
x=rndn(100,1);
times=10000;                 /* more than 200000 recommended(by normal version of GAUSS) */
pp=0.05;
call ADtest2(x,times,pp);

/* When meantimes=1, the same as above. As it is bigger(no limits), that will take more time. */
meantimes=10; 
call ADtest2a(x,times,meantimes,pp);



proc ADtest2(x,times,pp);
   local n,xbar,p,A2,stat,statm,i,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   x=sortc(x,1);
   p=cdfn((x-meanc(x))/stdc(x));
   A2=-n-sumc((2*seqa(1,1,n)-1)/n.*(ln(p)+ln(1-rev(p))));
   stat=A2*(1+0.75/n+2.25/n^2);
/* simulation */
   statm=zeros(times,1);
   i=1;
   do while i<=times;
      x=rndn(n,1);
      x=(x-meanc(x))/stdc(x);
      x=sortc(x,1);
      p=cdfn((x-meanc(x))/stdc(x));
      A2=-n-sumc((2*seqa(1,1,n)-1)/n.*(ln(p)+ln(1-rev(p))));
      statm[i]=A2*(1+0.75/n+2.25/n^2);
      i=i+1;
   endo;
   statm=sortc(statm,1);
   cval=statm[round(times*(1-pp))];
/* results */
   print "Anderson-Darling Test:";
   print/rz "n             =" n;
   print    "A*            =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%";
   retp(stat~cval);
endp;


/* Workable in Light version. (by using the means of each critical value many times) */

proc ADtest2a(x,times,meantimes,pp);
   local n,xbar,p,A2,stat,statm,i,j,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   x=sortc(x,1);
   p=cdfn((x-meanc(x))/stdc(x));
   A2=-n-sumc((2*seqa(1,1,n)-1)/n.*(ln(p)+ln(1-rev(p))));
   stat=A2*(1+0.75/n+2.25/n^2);
/* simulation */
   cval=0;
   j=1;
   do while j<=meantimes;
      statm=zeros(times,1);
      i=1;
      do while i<=times;
         x=rndn(n,1);
         x=(x-meanc(x))/stdc(x);
         x=sortc(x,1);
         p=cdfn((x-meanc(x))/stdc(x));
         A2=-n-sumc((2*seqa(1,1,n)-1)/n.*(ln(p)+ln(1-rev(p))));
         statm[i]=A2*(1+0.75/n+2.25/n^2);
         i=i+1;
      endo;
      statm=sortc(statm,1);
      cval=cval+statm[round(times*(1-pp))];
      j=j+1;
   endo;
   cval=cval/meantimes;
/* results */
   print "Anderson-Darling Test:";
   print/rz "n             =" n;
   print    "A*            =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%  (mean)";
   retp(stat~cval);
endp;