home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol131 / popstats.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1984-04-29  |  8.4 KB  |  216 lines
  1. PROCEDURE popstats(VAR a          : list;
  2.                    n1st,nlast     : INTEGER;
  3.                    VAR parameters : varrecord);
  4. {$c-    [ control C checking off] }
  5. {
  6. Comment  :  Popstats will calculate the mean , mode , median , variance ,
  7. standard deviation,standard error of the mean,coefficients of skewness 
  8. and kurtosis,standard errors of skewness and kurtosis,and range (smallest
  9. and largest values) of the population of observations in the array "a" (of 
  10. real numbers ) from a[n1st] to a[nlast].
  11. If you wish to analyse the whole array then just set n1st = 1 & nlast = N
  12. where N is the order of the array.
  13. The above statistics ,if no errors occurr ,are returned as the return value
  14. of the variable corresponding to "parameters".This variable is structured as
  15. a variable record of type varrecord which MUST be defined at a higher level
  16. exactly as follows :
  17.  
  18. TYPE
  19.     varrecord    = RECORD
  20.       CASE success : BOOLEAN OF
  21.                TRUE    : ( mean,
  22.                            mostfreq,
  23.                            middle,
  24.                            variance,
  25.                            stddevtn,
  26.                            stderrmn,
  27.                            skewness,
  28.                            kurtosis,
  29.                semedian,
  30.                seskewns,
  31.                sekurtss : REAL;
  32.                range    : ARRAY[1..2] OF REAL);
  33.                FALSE   : ( errmsg1,
  34.                            errmsg2,
  35.                            errmsg3,
  36.                            errmsg4  : BOOLEAN );
  37.                END; 
  38.  
  39. The record field success will indicate if ANY errors occurred.If success
  40. is FALSE then the other fields will indicate the type of error.
  41.                        errmsg1 : TRUE means n < 1;
  42.                        errmsg2 : TRUE means n < 4;
  43.                        errmsg3 : TRUE means array is uniform;
  44.  
  45. Externals required := FUNCTIONS median and select
  46.  
  47. }
  48. CONST
  49.     nbins        = 15; { (MUST be ODD) = 1 : number of bins that
  50.                                data are counted into when calculating mode.
  51.                                                 2 : scaling factor to code
  52.                                the range of observations when calculating
  53.                                sums powers of the deviations from the mean.}
  54.         large           = 150;{ value of n above which simpler calculations
  55.                                 are used for skewness and kurtosis }
  56. VAR
  57.     popular,
  58.     bindex,
  59.     n,i,bin#    : INTEGER;
  60.  
  61.         bin        : ARRAY[1..nbins] OF INTEGER;
  62.  
  63.     temp1,
  64.     temp2,
  65.     ratio1,
  66.     ratio2,
  67.     ratio3,
  68.     ratio4,
  69.     ratio5,
  70.     ratio6,
  71.     ratio7,
  72.     vcoded,
  73.     midval,
  74.     smallest,
  75.     biggest,
  76.         binwidth,
  77.     yc,ycsq,
  78.     ycqd,yc4d,
  79.     syc,sycsq,
  80.     sycqd,syc4d    : REAL;
  81.  
  82. BEGIN
  83. WITH parameters DO BEGIN
  84.  success := FALSE;
  85.  errmsg1 := FALSE;
  86.  errmsg2 := FALSE;
  87.  errmsg3 := FALSE;
  88.  errmsg4 := FALSE;
  89.  IF (median(a,n1st,nlast,midval) = FALSE) { checks list size & gets median }
  90.  THEN BEGIN
  91.        errmsg1 := TRUE;     { To indicate that n < 1 , median failed }
  92.       END
  93.  ELSE BEGIN
  94.        n        := nlast - n1st + 1;
  95.        IF n < 4
  96.        THEN errmsg2 := TRUE
  97.        ELSE BEGIN
  98.          smallest := a[n1st];
  99.          biggest  := smallest;
  100.          FOR  i   := n1st TO nlast DO
  101.            BEGIN
  102.             IF a[i] > biggest
  103.               THEN biggest := a[i]
  104.               ELSE BEGIN
  105.                     IF a[i] < smallest THEN smallest := a[i];
  106.                    END;
  107.            END;
  108.          IF biggest = smallest
  109.          THEN BEGIN
  110.                errmsg3 := TRUE; { To indicate that array is uniform! }
  111.               END
  112.          ELSE BEGIN             { now to calculate the parameters.   }
  113.                success  := TRUE;
  114.                range[1] := smallest;
  115.                range[2] := biggest;
  116.                binwidth := (biggest - smallest) / nbins;
  117.                FOR   i  := 1    TO nbins DO bin[i] := 0;
  118.                FOR   i  := n1st TO nlast DO
  119.                 BEGIN
  120.                  bin#      := TRUNC((a[i] - smallest) / binwidth) + 1;
  121.                  IF bin# > nbins THEN bin# := nbins; { to allow the biggest
  122.                                                        value to go into the 
  123.                                                        correct bin.         }
  124.                  bin[bin#] := bin[bin#] + 1;
  125.                 END;
  126.                bindex  := 1;
  127.                popular := bin[1];
  128.                FOR i := 2 TO nbins DO
  129.                 BEGIN
  130.                  IF bin[i] > popular
  131.                     THEN BEGIN
  132.                           bindex  := i;
  133.                           popular := bin[i];
  134.                          END;
  135.                 END;
  136.                mostfreq := smallest + (binwidth * (bindex - 0.5));{ mode is 
  137.                                                                     done ! }
  138.              { now get the rest of the parameters }
  139.                syc      := 0.0;
  140.                sycsq    := 0.0;
  141.                sycqd    := 0.0;
  142.                syc4d    := 0.0;
  143.                FOR  i := n1st TO nlast DO
  144.                  syc := syc + (a[i] - midval);
  145.                  { comment : by determining the sum of the differrences from
  146.                    the median you will end up with a much smaller number than
  147.                    if you add the raw values.This avoids loss of precision
  148.                    when n is large and the sum of the data can be many orders
  149.                    larger than any individual datum.All you have to do then
  150.                    is the reverse coding to get the true mean.}
  151.                mean := (syc/n) + midval;
  152.                syc  := 0.0;
  153.                FOR  i := n1st TO nlast DO
  154.                 BEGIN
  155.                  yc     := (a[i] - mean)/binwidth;
  156.                            { code deviations so they lie in a small range 
  157.                    about zero, to avoid the problems with raising 
  158.                    large or small #s to their 3rd & 4th powers.}
  159.                  ycsq  := yc * yc;
  160.                  ycqd  := yc * yc * yc;
  161.                  syc   := syc +  yc;
  162.                  sycsq := sycsq +  ycsq;
  163.                  sycqd := sycqd + ycqd;
  164.                  syc4d := syc4d + ycqd * yc;
  165.                 END;
  166.               vcoded     := sycsq /(n - 1);
  167.               variance   := vcoded * binwidth * binwidth;
  168.               stddevtn   := SQRT(variance);
  169.               middle     := midval;
  170.               stderrmn   := stddevtn/(SQRT(n));
  171.               semedian   := 1.2533 * stderrmn;
  172.               IF n > large { use simpler equations for skewness + kurtosis }
  173.                 THEN BEGIN
  174.                      skewness :=  sycqd/(n * vcoded * SQRT(vcoded));
  175.                      kurtosis := (syc4d/(n * SQR(vcoded))) - 3.0;
  176.                      seskewns := SQRT(6.0/n);
  177.                      sekurtss := SQRT(24.0/n)
  178.                      END
  179.                 ELSE BEGIN  {If n is not large then you have to correct the
  180.                              estimates of skewness and kurtosis for sample
  181.                              size bias.This complicates the equations by
  182.                              adding factors of the order of n to both their
  183.                              numerators and their denominators.This in turn
  184.                              creates significant risks of arithmetic errors
  185.                              such as overflow and loss of precision.To avoid
  186.                              these the calculations are done indirectly at a
  187.                              cost of space,time and readability,unfortunately.}
  188.  
  189.                      ratio1   := n/(n - 1);       {First set up these ratios..}
  190.                      ratio2   := (n - 1)/(n - 3); {..all of which will be ... }
  191.                      ratio3   := (n - 1)/(n - 2); {..approx. = 1.0 }
  192.                      ratio4   := ratio2 * ratio3;
  193.                      ratio5   := (n + 1)/(n - 2);
  194.                      ratio6   := n/(n + 1);
  195.                      ratio7   := ratio1 * ratio5;
  196.                      temp2    := sycqd/vcoded;
  197.                      temp1    := ratio1 * temp2;
  198.                      temp1    := temp1/(n - 2);
  199.                      skewness := temp1/(SQRT(vcoded));
  200.                      temp1    := syc4d / (SQR(vcoded));
  201.                      temp2    := ratio7 * temp1 / (n -3);
  202.                      kurtosis := temp2 - (3.0 * ratio4);
  203.                      seskewns := SQRT((6.0*ratio6*ratio3)/(n+3));
  204.                      temp2    := 24.0 / (n + 3);
  205.                      temp1    := ratio4 * temp2;
  206.                      sekurtss := SQRT(temp1/(n +5));
  207.                     END;
  208.             END; { of : calculation of parameters }
  209.       END; { of : if n < 4 }
  210.   END; { of : IF median }
  211.  END;  { of : WITH parameters DO }
  212. END;   { of : PROCEDURE popstats }
  213.  
  214.  
  215.  
  216.