/* Pearson chi-square test for normality */
/* This program is a GAUSS version of 'pearson.test' */
/* by Juergen Gross in library 'nortest' of R.       */
new; cls;
x=rndn(100,1);
call PearsonChisq(x);
print;
call PearsonChisqAdj(x);



proc PearsonChisq(x);
   local n,nc,num,c,prob,xp,h,p,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   nc=ceil(2*(n^(2/5)));
   num=floor(1+nc*cdfn((x-meanc(x))/stdc(x)));
   c=tabulate(num,nc);
   prob=1/nc*ones(nc,1);
   xp=n*prob;
   h=((c-xp)^2)./xp;
   p=sumc(h);
   pval=cdfchic(p,nc-1);
   print "stat=" p;
   print "pval=" pval;
   retp(p);
endp;

proc PearsonChisqAdj(x);
   local n,nc,num,c,prob,xp,h,p,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   nc=ceil(2*(n^(2/5)));
   num=floor(1+nc*cdfn((x-meanc(x))/stdc(x)));
   c=tabulate(num,nc);
   prob=1/nc*ones(nc,1);
   xp=n*prob;
   h=((c-xp)^2)./xp;
   p=sumc(h);
   pval=cdfchic(p,nc-1-2);
   print "stat=" p;
   print "pval=" pval "(adjusted)";
   retp(p);
endp;

proc tabulate(b,nb);
   local x,i;
   x=zeros(nb,1);
   i=1;
   do while i<=nb;
      x[i]=sumc(b.==i);
      i=i+1;
   endo;
   retp(x);
endp;