/* Frequency(monobit)Test[of NIST SP800-22] for 01 random sequence */
/* Based on Katsuichi Hirose's paper(2004) in Japanese. */
new; cls;
x=rndu(10000,1);
x=(x.>0.5);          /* x=rev(x); if you need the reverse order of it. */
call frequencytest(x);





proc frequencytest(x);
   local n,xi,Sn,stat,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   xi=2*x-1;                        /* -1 or 1 */
   Sn=sumc(xi);
   stat=Sn/sqrt(n);                 /* Not taking abs() here on purpose. */
   pval=erfc(abs(stat)/sqrt(2));    /* the same as 2*(1-cdfn(abs(stat))) */
/* results */
   print "Frequency Test:";
   print/rz "n   =" n;
   print    "stat=" stat "(with sign)";
   print    "pval=" pval;
   retp(stat);
endp;