/* Combination test(Greater than 0.5 or not) for U(0,1) sequence */
new; cls;
x=rndu(1000,1);
length=10;
call combtest(x,length);



proc combtest(x,length);
   local n,ng,y,Fi,k,yn,count,i,x2,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   if not(x>=0 and x<=1);
      errorlog "Error: This is not U(0,1).";
      retp(".");
   endif;
   n=rows(x);
   ng=floor(n/length);
   y=(x.>0.5);
   Fi=zeros(length+1,1);
   k=0;
   do while k<=length;
      Fi[k+1]=ng*nCr(length,k)*(1/2)^length;
      k=k+1;
   endo;
   yn=zeros(ng,1);
   i=1;
   do while i<=ng;
      yn[i]=sumc(y[length*(i-1)+1:length*i]);
      i=i+1;
   endo;
   count=zeros(length+1,1);
   i=0;
   do while i<=length;
      count[i+1]=sumc(yn.==i);
      i=i+1;
   endo;
   x2=sumc(((count-Fi)^2)./Fi);
   pval=cdfchic(x2,length);
/* results */
   print "Combination Test(x>0 or not):";
   print/rz "n     =" n;
   print/rz "length=" length;
   print    "x2    =" x2;; print/lz "  d.f.=" length;
   print    "pval  =" pval;
/* graphs */
   library pgraph;
   graphset;
   pqgwin auto;
   _plctrl=1;
   _plegctl={2 5 6 5.5};
   _plegstr="Theoretical\000Actual Counts";
   xlabel("# of 1's");
   ylabel("groups");
   xy(seqa(0,1,length+1),Fi~count);
   retp(x2);
endp;

proc nCr(n,r);
   if r==0 or n==r;
      retp(1);
   endif;
   retp( exp(sumc(ln(seqa(n,-1,n)))-sumc(ln(seqa(r,-1,r)))-sumc(ln(seqa(n-r,-1,n-r)))) );
endp;