/* Length of Runs Test for U(0,1) sequence (The null of this test tends to be rejected.) */
/* 1 if up, 0 otherwise. Counts each length of runs for both up-runs and down-runs.      */
/* Attention! The formula by Loren Wahl on the web is wrong. His numerical example(N=40) */
/*            is based on the correct formula below.                                     */
/*            The note and formula in 'Ransuu no Chishiki' by Kazumasa Wakimoto(1970) on */
/*            the web is also wrong. Not either runs but both for this fromula.          */
/*            The correct formula is Fi=2/((i+3)!)*(n*(i^2+3*i+1)-(i^3+3*i^2-i-4)).      */
new; cls;
x=rndu(100,1);
call lengthofruns(x);


proc lengthofruns(x);
   local n,n1,d,y,c,count,i,Fi,x2,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   d=x[2:n]-x[1:n-1]; 
   y=(d.>0);                       /* 1 if x[i]>x[i-1], 0 otherwise. */
/* count each length of runs */
   n1=rows(y);
   c=miss(0,0);
   count=1;
   i=2;
   do while i<=n1;
      if y[i-1]/=y[i] and i/=n1;
         c=c|count;
         count=0;
      elseif y[i-1]/=y[i] and i==n1;
         c=c|count|1;
      elseif y[i-1]==y[i] and i==n1;
         count=count+1;
         c=c|count;
      endif;
      count=count+1;
      i=i+1;
   endo;
   c=c[2:rows(c)];
   count=zeros(n-2,1);
   i=1;
   do while i<=n-2;
      count[i]=sumc(c.==i);
      i=i+1;
   endo;
   Fi=zeros(n-2,1);
   i=1;
   do while i<=n-2;
      Fi[i]=2/((i+3)!)*(n*(i^2+3*i+1)-(i^3+3*i^2-i-4));
      i=i+1;
   endo;
   x2=((count-Fi)^2)./Fi;
   x2=delif(x2,x2.==miss(0,0));
   x2=sumc(x2);
   pval=cdfchic(x2,2);
/* results */
   print "Length of Runs Test:";
   print/rz "n   =" n;
   print    "x2  =" x2;
   print    "pval=" pval;
/* graphs */
   library pgraph;
   graphset;
   pqgwin auto;
   _plctrl=1;
   _plegctl={2 5 6 5.5};
   _plegstr="Theoretical\000Actual Counts";
   xtics(1,10,1,0);
   xlabel("length of runs");
   xy(seqa(1,1,10),Fi[1:10]~count[1:10]);
   retp(x2);
endp;