/* D'Agostino's D */
/* Based on C.Minotani(1992) pp.38-39 and checked by 'gofnorm' in library 'fBasics' of R. */   
new; cls;
x=rndn(100,1);
call DAgostinoD(x);



proc DAgostinoD(x);
   local n,c,D,Y;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   x=sortc(x,1);
   n=rows(x);
   c=seqa(1,1,n)-1/2*(n+1);
   D=sumc(c.*x)/(n^(3/2)*sqrt(sumc((x-meanc(x))^2)));
   Y=sqrt(n)*(D-0.28209479)/0.02998598;
/* results */
   print "D'Agostino's D:";
   print/rz "n   =" n;
   print    "stat=" Y;
   retp(Y);
endp;