/* Some useful statistical procedues. For the GRBL2 package

  This should be included in the program file -- put
  #INCLUDE GRBL2_STATS.INC ;
 at the END of your program file.

*/


@ Compute basic statistics

G_STATS will use the following global variables. 
These are typically set by processing a GRBL2_BATCH command file.

The following globals are required:
   FILENAME : Name of gauss data file. 

   NX    : Number of variables to compute statistics for.
   XNAME : variables to compute statistics for. 
   XINDEX : pointers, into FILENAME's columns, of the XNAME variables

   SLCTLIST   :  Observation selection critieria. 0 to use all obs

   WEIGHT     :  Weight variable.
   WINDEX     :  Weight index (0 if unweighted)

The following are optional. If ID is specified, some stats on panel sizes
are produced.
   IDNAME     : If specified, variable name of  "panel" identifier. 
		 or, a "#"
   IDINDEX    : If specified, index into FILENAME.
		Of, if IDNAME="#", the number of rows in "balanced" panels.

  BYVAR      : If specified, variable name of the "by" identifier
  BYINDEX    : Index for above

  STATLIST      : Space delimited string of statistics to be generated.
	       If not specified, a basic set of values are computed.
	       The following statistics are supported:
		  MEAN MEDIAN MIN MAX N NMISS SD SUM 1 5 10 25 50 75 90 95 99 
               Or
		   *  
               Or, if blank (not specified)
		  MEAN SD MIN MAX N 

	       to display all the above statistics.
               Where the 1... 99 refer to quintiles. 
               I.e.; 5 means the value at which 5% of the observations are less then.
	       Data is displayed in order entered in STATS.
               Only 5 items are displayed per table (hence, if you specify 
	       everything, 4 tables are created).
	       Or, include a " $" to force a "table break".
	       For example:
		  STATS MEAN SD SUM $ MEDIAN 1 5 95 99 ;
	       Would produced 2 tables, the first with three columns of statistics, the
	       2nd with five.

   TITLEU   : One or more line title to display.

@
g_stats:     

clear newmat,fmt,nname,ddat,inthis,stat_vals,tabledo,ax,axp,axps,statsdo;
clear lstatlist;

open f1=^filename ;

lstatlist=g_strip(upper(statlist)) ;
if lstatlist$=="*" ;
 lstatlist="MEAN MEDIAN MIN MAX N NMISS SD SUM 1 5 10 25 50 75 90 95 99 "; 
endif;
if lstatlist$=="" ; lstatlist="MEAN MIN MAX SD N "; endif;


vlist="MEAN"|"MEDIAN"|"MIN"|"MAX"|"N"|"NMISS"|"SD"|"SUM"|"1"|"5"|"10"|"25"|"50"|"75"|"90"|"95"|"99"; 

if f1 < 0 ;
  call grbl2_error("ERROR:  no FILENAME specified ",1);
  stop;
endif;

dat=readr(f1,rowsf(f1));
nr_orig=rowsf(f1);  

call close(f1);

useme=g_doslct(slctlist,filename);
if scalmiss(useme);
    call grbl2_error(" ERROR: not able to open file " $+ filename,1);
      stop;
endif;

doby=0;
if missrv(byindex,0) ne 0 ;
   Id2index=byindex;
   id2name=byvar ;
   doby=1 ;
endif;
used_obs=g_mkdata("X"|"ID"|"ID2"|"WEIGHT",miss(1,1),useme,0);

nr_sel=rows(x);
nr_nomiss=rows(packr(x)) ;


if doby==1 ;
   bigx=x; bigid=id; bigw=weight;
endif;


tabledo=zeros(20,6);      @ each row is a table, col is stat to use @
statsdo=zeros(20,1);     @ flag =1 if "do this stat" (corresponds to vlist@

doweight=0;
if missrv(wtindex,0)>0 ;
    doweight=1 ; 
endif;


ir=1; ic=1;
nstats=0;
do until lstatlist$=="";
  if ic>5 ;			@ max 5 values per row @
     ir=ir+1; ic=1; 
  endif;
  {aword,lstatlist}=token(lstatlist);
   aword=g_strip(aword);    
   if aword$=="$" ;
        if ic>1 ; 
          ir=ir+1;
   	  ic=1;
        endif;
  else;
     ii=indcv(aword,vlist);
     if scalmiss(ii); format /ld 1,8;
       call grbl2_error("ERROR: unknown statistic: "  $+ aword,1);
	stop;
     endif; 
     tabledo[ir,ic]=ii ;
     statsdo[ii]=1;
     tabledo[ir,6]=ic ;
     ic=ic+1;
     nstats=nstats+1;
   endif;
endo;

if ic==1 ; 
   ir=ir-1; 
endif;
if ir==0 ;
  call grbl2_error("ERROR: no stats specified ",1);
  stop;
endif;

tabledo=tabledo[1:ir,.];

stat_Vals=zeros(rows(xindex),20);   @ stats stored in here @

   
@  1     2    3   4  5   6    7   8  9  10 11 12 13  14  15 16 17  
 MEAN MEDIAN MIN MAX N NMISS SD  SUM 1  5  10 25 50  75  90 95 99 
@
dosort=sumc(statsdo[2|9|10|11|12|13|14|15|16|17]) ;

hlines={};
for j (1,rows(tabledo),1);
  inthis=tabledo[j,6];
  hline=" Variable "; 
  for jj (1,inthis,1);
     ii=tabledo[j,jj];
     aa=vlist[ii];
     if doweight==1 ;
       if ii==1 or ii==7 or ii==8 ;
	  aa="Wtd. "$+aa;
       endif;
     endif;
     aa=g_padwrd(aa,13,1);
     hline=hline$+aa;
  endfor;
  hlines=hlines$|hline;
endfor;


OUTPUT ON;
"***********************************************************";
 foo=g_parse_title(" " $+ titleu,"   ");
 $foo;  
"***********************************************************" ;
  format /ld 1,8;   ":  " ;;$datestr(date);; " " ;;$timestr(time);
  ": Dataset= " $filename;
   "   # of rows in dataset: " ;; format 10,0 ; nr_orig ;
   format 8,8;

  if not(missrv(slctlist,0) eq 0);  @ display other selection criteria @
       "Selection criteria: ";
         oo=g_vuslct(slctlist,FILEname,0);
       "   Rows removed using selection criteria: ";;
        format /ld 6,0 ; (nr_orig-nr_sel) ; 
   endif;
   if nr_sel>nr_nomiss;
       "   Rows  with NO missing values: ";;
        format /ld 6,0 ; (nr_nomiss) ; 
   endif;
 ?;
OUTPUT OFF;


if doby==0;
  output on ;
  goto stats_make_tables ;
endif;

bylist=unique(id2,1);
iby=1;

OUTPUT ON;
" BY variable ";; format /ld 1,8; $byvar;; " has " ;; format /ld 6,0;
rows(bylist);; " unique values.";
output off ;

ugly0:;		@ gosub in gosub seems iffy.. so use ugly gotos instead @
    bydx=g_mkindx(id2.==bylist[iby]);
    x=bigx[bydx,.];
    if doweight==1 ; weight=bigw[bydx,.] ; endif;
    if missrv(idindex,0) ne 0 ; id=bigid[bydx,.]; endif;

   output on ;
   ?; ?;  "   ------------------------------------ ";
    "For ";; format /ld 1,8; $byvar;;" =  ";; format /ld 14,4;bylist[iby];;
     ", # of obs: ";;format /ld 8,0; rows(x);
     ?;
    goto stats_make_tables ;
   OUTPUT OFF ;

ugly1: ;

if doby==0;		@an UGLY hack @
 output off ;
 return;
endif;

if iby== rows(bylist); return; endif;
iby=iby+1;
goto ugly0 ;



return; 


@----------------------------@
@ goto routine to make the tables @
stats_make_tables: ;

 
allmiss=zeros(cols(x),1);        @ rows=obs, col=var @
 
for k (1,cols(x),1);	@ compute stats a column at a time @

  ax=x[.,k];    @ this vars data @
  nax=rows(ax);
  
  if doweight ;
     axp=packr(ax~weight);
     if scalmiss(axp);
          allmiss[k]=2;
          continue ;
     endif;
     awp=axp[.,2]; axp=axp[.,1];
  else ;
     axp=packr(ax);
     if scalmiss(axp);
        allmiss[k]=1;
        continue ;
     endif;
     awp=1;
  endif;
  naxp=rows(axp);

  if dosort>0;
    axps=sortc(axp,1);
  endif;

  didsort=0;			@ set to 1 if sorted @
  if statsdo[1]==1 ;	@ mean@
     if doweight==1 ;
        if wtrep==1 ;
           stat_vals[k,1]=sumc(axp.*awp)/sumc(awp);
        else;
           stat_vals[k,1]=meanc(axp.*awp);
        endif;
    else;
           stat_vals[k,1]=meanc(axp);
    endif;

  endif;

  if statsdo[2]==1 ;	@median @
      aat=ceil(0.5*naxp);
      stat_vals[k,2]=axps[aat];
  endif;

  if statsdo[3]==1 ;	 @ min @
       stat_vals[k,3]=minc(axp);
  endif;

  if statsdo[4]==1 ;	@max @
       stat_vals[k,4]=maxc(axp);
  endif;

  if statsdo[5]==1 ;	@ N @
       stat_vals[k,5]=naxp; 
  endif;

  if statsdo[6]==1 ;	@nmiss @
       stat_vals[k,6]=(nax-naxp);
  endif;

  if statsdo[7]==1 ;	@sd @
       if doweight==0;
          stat_vals[k,7]=stdc(axp);
       else;
          if wtrep==1 ;
              tmean=meanc(axp.*awp)/sumc(awp);
              diff=axp-tmean;
              diff2=diff.^2 ;
              diff2=diff2.*awp;
              avar=sumc(diff2)./sumc(awp);
          else;
              tmean=meanc(axp.*awp);
              diff=(awp.*axp)-tmean;
              diff2=diff.^2 ;
              avar=sumc(diff2)./rows(diff2);
          endif;
          stat_vals[k,7]=sqrt(avar);
       endif;
  endif;

  if statsdo[8]==1 ;	@sum @
       stat_vals[k,8]=sumc(axp.*awp);
  endif;

  if statsdo[9]==1 ;	@ 1 @
      aat=ceil(0.01*naxp);
      stat_vals[k,9]=axps[aat];
  endif;

  if statsdo[10]==1 ;      @ 5 @
      aat=ceil(0.05*naxp);
       stat_vals[k,10]=axps[aat];
  endif;

  if statsdo[11]==1 ;      @ 10 @
      aat=ceil(0.10*naxp);
       stat_vals[k,11]=axps[aat];
  endif;
  
  if statsdo[12]==1 ;      @ 25 @
      aat=ceil(0.25*naxp);
      stat_vals[k,12]=axps[aat];
  endif;

  if statsdo[13]==1 ;      @ 50 @
      aat=ceil(0.50*naxp);
      stat_vals[k,13]=axps[aat];
  endif;

  if statsdo[14]==1 ;      @ 75 @
      aat=ceil(0.75*naxp);
      stat_vals[k,14]=axps[aat];
  endif;

  if statsdo[15]==1 ;      @ 90 @
      aat=ceil(0.90*naxp);
      stat_vals[k,15]=axps[aat];
  endif;

  if statsdo[16]==1 ;      @ 95 @
       aat=ceil(0.95*naxp);
       stat_vals[k,16]=axps[aat];
  endif;

  if statsdo[17]==1 ;      @ 99 @
      aat=ceil(0.99*naxp);
      stat_vals[k,17]=axps[aat];
  endif;


endfor;


@ stats are computed. Now display @


@ do for each line of tabledo @


for j (1,rows(tabledo),1);
  inthis=tabledo[j,6];

  ddat=zeros(rows(stat_vals),inthis);

  if rows(tabledo)>1 ; format /ld 5,0; " Table ";;j;;" of " ;; rows(tabledo); endif;
  " =============================================================================";
   $strsect(hlineS[J],1,79);
  " =============================================================================";
  
  for kx (1,rows(stat_vals),1);
     for jj (1,inthis,1);
        tmp1=tabledo[j,jj];
        ddat[kx,jj]=stat_vals[kx,tmp1];
     endfor;
  endfor;

  newmat=zeros(rows(xname),1)~xname~ddat;
  mask=0~0~1~1~1~1~1;
  let fmt[7,3]=
             "-*.*s" 1 1
             "-*.*s" 8 8
             "*.*lg " 13 5
             "*.*lg " 13 5
             "*.*lg " 13 5
             "*.*lg " 13 5
             "*.*lg " 13 5 ;

   for jj (1,rows(xname),1);
     if allmiss[jj] eq 1; format /ld 8,8 ;
        " ";;$newmat[jj,2];; "    :: no non-missing values ";
     elseif allmiss[jj] eq 2 ;  format /ld 8,8 ;
        " ";;$newmat[jj,2];; "    :: no non-missing weighyted values ";
     else;
        d=printfm(newmat[jj,.],mask[1:inthis+2],fmt[1:inthis+2,.]);?;
     endif;
   endfor;
 ?;
endfor;


 ? ;
 if doweight==1 ; 
    if wtrep==1 ; 
       "  Weight variable (replication): " ;; format /ld 1,8; $wtname ;;
    else;
       "  Weight variable (scale): " ;; format /ld 1,8; $wtname ;; 
   endif;
   "( sum= " ;; format /ld 14,5 ; g_addcomma(sumc(packr(weight)),2) ;;" )";

  endif;


@ now do panel stuff @

if missrv(idindex,0)>0;
  if idtype$=="COUNT";
    panels=g_pnlind_count(id);
  else;
    panels=g_pnlind(id);
  endif;

  ilen=(panels[.,2]-panels[.,1])+1;

  "  ID variable (";; format /ld 1,8;$IDNAME;;")";
  "    Number of panels: ";; format /ld 7,0; rows(panels);

  ii=minindc(ilen);
  ii2=panels[ii,1];
  ii3=id[ii2];
  "    Smallest: ";;minc(ilen);;"(";;ii3;;")";

  ii=maxindc(ilen);
  ii2=panels[ii,1];
  ii3=id[ii2];
  "    Largest: ";;maxc(ilen);;"(";;ii3;;")";
  
  "    Average size: ";;format /ld 8,2 ; meanc(ilen);
endif;


PRINT ;

goto ugly1 ;



/**********************************************************/
@ Show a frequency table


G_FREQ will use the following global variables. 
These are typically set by processing a GRBL2_BATCH command file.

The following globals are required:
   FILENAME : Name of gauss data file. 
  
   Y, YINDEX : The variable to display

   SLCTLIST   :  Observation selection critieria. 0 to use all obs


The following are optional. If ID is specified, some stats on panel sizes
are produced.

  RANGE : The categories to compute within
 		RANGE -DEFAULT  :   Use the default ranges
	        RANGE AUTO=N    :   Program divides into N ranges
	        RANGE -ALL      :  Use all values (yields a count of occurences of each value)
	        RANGE n1 n2 n3  : n1 .. are categories values:

	Note:
		RANGE -DEFAULT
	and
	     RANGE 0 1 2 3 4 5 10 25 50   100 250 500 1000 2000  1000000 ;
        are the same thing.


@

g_freq: ;

open f1=^filename ;

dat=readr(f1,rowsf(f1));
nr_orig=rowsf(f1);

call close(f1);

useme=g_doslct(slctlist,filename);
if scalmiss(useme);
    call grbl2_error(" ERROR: not able to open file " $+ filename,1);
      stop;
endif;

used_obs=g_mkdata("Y",miss(1,1),useme,0);

nr_sel=rows(y);
nr_nomiss=rows(packr(y)) ;

if progdo==-1;		@ -1 means "ALL value @
   range=sortc(unique(y,1),1);
elseif progdo>0 ;
   tmp=sortc(y,1);
   nr=rows(tmp);
   ar=seqa(1,1,progdo)/(progdo+1);
   ar=ceil(ar*nr);
   range=tmp[ar];
   range=sortc(unique(range,1),1);
endif;
rdif=range[2:rows(range)]-range[1:(rows(range)-1)];
if minc(rdif) le 0 ;
  call grbl2_error("Error: Range intervals are NOT strictly increasing ",1);
  stop;
endif;



ctot=zeros(rows(range),1);
min1a=1e+38;max1a=-min1a;
ntot=rows(y);

min1=minc(y);max1=maxc(y);
min1a=minc(min1|min1a);max1a=maxc(max1|max1a) ;

C=COUNTS(y,RANGE);
CTOT=CTOT+C;

output on ;
pcts=ctot/ntot;
ii=rows(pcts);cumpcts=zeros(ii,1);
for ij (1,ii,1);
  cumpcts[ij,1]=sumc(pcts[1:ij,1]);
endfor;



OUTPUT ON;
"***********************************************************";
 foo=g_parse_title(" " $+ titleu,"   ");
 $foo;
"***********************************************************" ;
  format /ld 1,8;   ":  " ;;$datestr(date);; " " ;;$timestr(time);
  ": Dataset= " $filename;
   "   # of rows in dataset: " ;; format 10,0 ; nr_orig ;
   format 8,8;

  if not(missrv(slctlist,0) eq 0);  @ display other selection criteria @
       "Selection criteria: ";
         oo=g_vuslct(slctlist,FILEname,0);
       "   Rows removed using selection criteria: ";;
        format /ld 6,0 ; (nr_orig-nr_sel) ; 
   endif;
   if nr_sel>nr_nomiss;
       "   Rows  with NO missing values: ";;
        format /ld 6,0 ; (nr_nomiss) ; 
   endif;
?;

"  :: Variable: " $yname;;" ; # used observations: " ;; format /ld 10,0; ntot;

let fmt[4,3]=
"*.*lg " 13 7
"*.*lf " 9 0
"*.*lg " 9 4
"*.*lg " 9 4 ;
let mask[1,4]= 1 1 1 1;

if progdo==-1;
"
   Value       : Count  : % of    : Cumlative
               :        : total   :   %
               :        :         :
---------------:--------:---------:----------- " ;


else;
"
    End of     : Count  : % of    : Cumlative
     range     :        : total   :   %
 (inclusive)   :        :         :
---------------:--------:---------:----------- " ;
endif;

d=printfm(range~ctot~pcts~cumpcts,mask,fmt);
 "         =========";  format /LD 13,3 ;
 " Range: " min1a " to " max1a ;
?;?;

output off;


return;


