/* Geary statistic with critical value by simulation */
/* Notice: I almost verify the table in Y.Shibata(1981). */         
new; cls;
x=rndn(100,1);
times=10000;                 /* more than 200000 recommended(by normal version of GAUSS) */
pp=0.05;
call Geary2(x,times,pp);
call Geary2med(x,times,pp);

/* When meantimes=1, the same as above. As it is bigger(no limits), that will take more time. */
meantimes=10; 
print;
call Geary2a(x,times,meantimes,pp);
call Geary2meda(x,times,meantimes,pp);





proc Geary2(x,times,pp);
   local n,stat,statm,i,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   stat=sumc(abs(x-meanc(x)))/sqrt(n*sumc((x-meanc(x))^2));
/* simulation */
   statm=zeros(times,1);
   i=1;
   do while i<=times;
      x=rndn(n,1);
      x=(x-meanc(x))/stdc(x);
      statm[i]=sumc(abs(x-meanc(x)))/sqrt(n*sumc((x-meanc(x))^2));
      i=i+1;
   endo;
   statm=sortc(statm,1);
   cval=statm[round(times*(1-pp))];
/* results */
   print "Geary Test:";
   print/rz "n             =" n;
   print    "stat ( G1 )   =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%";
   retp(stat~cval);
endp;


proc Geary2med(x,times,pp);
   local n,stat,statm,i,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   stat=sumc(abs(x-median(x)))/sqrt(n*sumc((x-meanc(x))^2));
/* simulation */
   statm=zeros(times,1);
   i=1;
   do while i<=times;
      x=rndn(n,1);
      x=(x-meanc(x))/stdc(x);
      statm[i]=sumc(abs(x-median(x)))/sqrt(n*sumc((x-meanc(x))^2));
      i=i+1;
   endo;
   statm=sortc(statm,1);
   cval=statm[round(times*(1-pp))];
/* results */
   print "Geary Test(median version):";
   print/rz "n             =" n;
   print    "stat ( G2 )   =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%";
   retp(stat~cval);
endp;


/* Workable in Light version. (by using the means of each critical value many times) */

proc Geary2a(x,times,meantimes,pp);
   local n,stat,statm,i,j,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   stat=sumc(abs(x-meanc(x)))/sqrt(n*sumc((x-meanc(x))^2));
/* simulation */
   cval=0;
   j=1;
   do while j<=meantimes;
      statm=zeros(times,1);
      i=1;
      do while i<=times;
         x=rndn(n,1);
         x=(x-meanc(x))/stdc(x);
         statm[i]=sumc(abs(x-meanc(x)))/sqrt(n*sumc((x-meanc(x))^2));
         i=i+1;
      endo;
      statm=sortc(statm,1);
      cval=cval+statm[round(times*(1-pp))];
      j=j+1;
   endo;
   cval=cval/meantimes;
/* results */
   print "Geary Test:";
   print/rz "n             =" n;
   print    "stat ( G1 )   =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%  (mean)";
   retp(stat~cval);
endp;


proc Geary2meda(x,times,meantimes,pp);
   local n,stat,statm,i,j,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   stat=sumc(abs(x-median(x)))/sqrt(n*sumc((x-meanc(x))^2));
/* simulation */
   cval=0;
   j=1;
   do while j<=meantimes;
      statm=zeros(times,1);
      i=1;
      do while i<=times;
         x=rndn(n,1);
         x=(x-meanc(x))/stdc(x);
         statm[i]=sumc(abs(x-median(x)))/sqrt(n*sumc((x-meanc(x))^2));
         i=i+1;
      endo;
      statm=sortc(statm,1);
      cval=cval+statm[round(times*(1-pp))];
      j=j+1;
   endo;
   cval=cval/meantimes;
/* results */
   print "Geary Test(median version):";
   print/rz "n             =" n;
   print    "stat ( G2 )   =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%  (mean)";
   retp(stat~cval);
endp;