/* D'Agostino's D with critical value by simulation */ 
/* Notice: I almost verify the table in D'Agostino(1972). */           
new; cls;
x=rndn(100,1);
times=10000;                 /* more than 200000 recommended(by normal version of GAUSS) */
pp=0.05;
call DAgostino2(x,times,pp);

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



proc DAgostino2(x,times,pp);
   local n,c,D,stat,statm,i,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   x=sortc(x,1);
   c=seqa(1,1,n)-1/2*(n+1);
   D=sumc(c.*x)/(n^(3/2)*sqrt(sumc((x-meanc(x))^2)));
   stat=sqrt(n)*(D-0.28209479)/0.02998598;
/* 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);
      c=seqa(1,1,n)-1/2*(n+1);
      D=sumc(c.*x)/(n^(3/2)*sqrt(sumc((x-meanc(x))^2)));
      statm[i]=sqrt(n)*(D-0.28209479)/0.02998598;
      i=i+1;
   endo;
   statm=sortc(statm,1);
   cval=statm[round(times*pp)]~statm[round(times*(1-pp))];
/* results */
   print "D'Agostino's D:";
   print/rz "n             =" n;
   print    "stat          =" stat;
   print    "                          Lower            Upper";
   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 DAgostino2a(x,times,meantimes,pp);
   local n,c,D,stat,statm,i,j,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   x=sortc(x,1);
   c=seqa(1,1,n)-1/2*(n+1);
   D=sumc(c.*x)/(n^(3/2)*sqrt(sumc((x-meanc(x))^2)));
   stat=sqrt(n)*(D-0.28209479)/0.02998598;
/* simulation */
   cval=0~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);
         c=seqa(1,1,n)-1/2*(n+1);
         D=sumc(c.*x)/(n^(3/2)*sqrt(sumc((x-meanc(x))^2)));
         statm[i]=sqrt(n)*(D-0.28209479)/0.02998598;
         i=i+1;
      endo;
      statm=sortc(statm,1);
      cval=cval+(statm[round(times*pp)]~statm[round(times*(1-pp))]);
      j=j+1;
   endo;
   cval=cval/meantimes;
/* results */
   print "D'Agostino's D:";
   print/rz "n             =" n;
   print    "stat          =" stat;
   print    "                          Lower            Upper";
   print    "critical value=" cval;; print/rz "         at" pp*100 "%  (mean)";
   retp(stat~cval);
endp;