/* Test for the longest run of ones in a block[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 longestrun(x);



proc longestrun(x);
   local n,M,K,Nf,pii,NN,maxrun,y,c,count,i,j,vi,x2,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   if n<128;
      errorlog "ERROR: This test requires the sample of at least 128.";
   elseif n<6272;
      M=8; K=3; Nf=16; pii={0.2148,0.3672,0.2305,0.1875,0,0,0};
   elseif n<750000;
      M=128; K=5; Nf=49; pii={0.1174,0.2430,0.2493,0.1752,0.1027,0.1124,0};
   else;
      M=10000; K=6; Nf=75; pii={0.0882,0.2092,0.2483,0.1933,0.1208,0.0675,0.0727};
   endif;
   NN=floor(n/M);
   maxrun=zeros(NN,1);
   j=1;
   do while j<=NN;
      y=x[(j-1)*M+1:j*M];
      c=miss(0,0);
      count=1;
      i=1;
      do while i<=M-1;
         if y[i]/=y[i+1] and i/=(M-1);
            c=c|count;
            count=0;
         elseif y[i]/=y[i+1] and i==(M-1);
            c=c|count|1;
         elseif y[i]==y[i+1] and i==(M-1);
            count=count+1;
            c=c|count;
         endif;
         count=count+1;
         i=i+1;
      endo;
      c=c[2:rows(c)];
      maxrun[j]=maxc(c);
      j=j+1;
   endo;
   vi=zeros(7,1);
   if M==8;
      i=1;
      do while i<=3;
         vi[i]=sumc(maxrun.==i);
         i=i+1;
      endo;
      vi[4]=sumc(maxrun.>=4);
   elseif M==128;
      vi[1]=sumc(maxrun.<=4);
      i=2;
      do while i<=5;
         vi[i]=sumc(maxrun.==(i+3));
         i=i+1;
      endo;
      vi[6]=sumc(maxrun.>=9);
   else;
      vi[1]=sumc(maxrun.<=10);
      i=2;
      do while i<=6;
         vi[i]=sumc(maxrun.==(i+9));
         i=i+1;
      endo;
      vi[7]=sumc(maxrun.>=16);
   endif;
   x2=sumc((vi[1:K+1]-Nf*pii[1:K+1])^2/(Nf*pii[1:K+1]));
   pval=cdfchic(x2,K);
/* results */
   print "Test of Longest Run in a Block:";
   print/rz "n   =" n;
   print    "x2  =" x2;
   print    "pval=" pval;
   retp(x2);
endp;