/* Test for Frequency within a Block[of NIST SP800-22] for 01 random sequence */
/* Based on Katsuichi Hirose's paper(2004) in Japanese. */
/* Notice: I think the degree of freedom here is N-1 rather than N itself. */
/* However, it does not matter as N goes to a big number we usually use.   */
new; cls;
x=rndu(10000,1);
x=(x.>0.5);          /* x=rev(x); if you need the reverse order of it. */
M=100;
call blockfreq(x,M);





proc blockfreq(x,M);
   local n,NN,pii,i,x2,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   NN=floor(n/M);
   pii=zeros(NN,1);
   i=1;
   do while i<=NN;
      pii[i]=sumc(x[(i-1)*M+1:i*M])/M;
      i=i+1;
   endo;
   x2=4*M*sumc((pii-0.5)^2);
   pval=cdfchic(x2,NN-1);
/* results */
   print "Test for Frequency within a Block:";
   print/rz "n   =" n;
   print/rz "M   =" M;
   print    "x2  =" x2;; print/lz "  d.f.=" NN-1;
   print    "pval=" pval;
   retp(x2);
endp;