home *** CD-ROM | disk | FTP | other *** search
/ norge.freeshell.org (192.94.73.8) / 192.94.73.8.tar / 192.94.73.8 / pub / computers / cpm / alphatronic / JRTPAS40.ZIP / JSTAT.PAS < prev    next >
Pascal/Delphi Source File  |  1999-04-03  |  3KB  |  129 lines
  1. (*  JSTAT ver 4.0    JRT Systems        *)
  2. (*                        *)
  3. (* jstat computes several basic statistics on    *)
  4. (* an input array.                *)
  5. (*                        *)
  6. (* parameters:                    *)
  7. (*    n - the number of data items in the    *)
  8. (*        input array            *)
  9. (*    x - the input array of real numbers,    *)
  10. (*        may be up to 1000 elements,    *)
  11. (*        actual variable in calling pgm    *)
  12. (*        may be much smaller array    *)
  13. (*    r - the computed statistics are stored    *)
  14. (*        in this record            *)
  15.  
  16.  
  17. EXTERN
  18.  
  19. TYPE
  20. jstat_interface = RECORD
  21.        mean, standard_deviation, variance, skewness, kurtosis, m1, 
  22.               m2, m3, m4 : real;
  23.        END;
  24. jstat_array =  ARRAY [1..1000] OF real;
  25.  
  26.  
  27. {=================================================================}
  28. PROCEDURE jstat (n : integer; VAR x : jstat_array; VAR r : jstat_interface);
  29.        
  30. VAR
  31. i : integer;
  32. total_x, total_x2, total_x3, total_x4 : real;
  33.  
  34.  
  35. {=================================================================}
  36. FUNCTION cube (x : real) : real;
  37.  
  38. BEGIN
  39. cube := x * sqr(x);
  40. END;
  41.  
  42.  
  43. {=================================================================}
  44. FUNCTION sqrt (x : real) : real;
  45. VAR
  46. sq, a, b : real;
  47. exponent, i : integer;
  48. zap : RECORD
  49.        
  50.        CASE integer OF
  51.        0 : (num : real);
  52.        1 : (ch8 :  ARRAY [1..8] OF char);
  53.        END;
  54.  
  55. BEGIN
  56. IF x = 0.0 THEN
  57.        sqrt := 0.0
  58. ELSE
  59.        BEGIN
  60.        sq := abs(x);
  61.        zap.num := sq;
  62.        exponent := ord(zap.ch8[1]);
  63.        exponent := (exponent DIV 2) + 32;
  64.        zap.ch8[1] := chr(exponent);
  65.        a := zap.num;
  66.        b := 0;
  67.        i := 0;
  68.        WHILE a <> b DO
  69.               BEGIN
  70.               b := sq / a;
  71.               a := (a + b) / 2;
  72.               i := i + 1;
  73.               IF i > 4 THEN
  74.                      BEGIN
  75.                      i := 0;
  76.                      IF abs(a - b) < (1.0e-12 * a) THEN
  77.                             a := b;
  78.                      END;
  79.               END;
  80.        sqrt := a;
  81.        END;  (* else *)
  82. END;  (* sqrt *)
  83.  
  84.  
  85. {=================================================================}
  86. PROCEDURE totals;
  87. VAR
  88. i : integer;
  89. tx, tx2, tx3, tx4 : real;
  90. sum_x, mean : real;
  91.  
  92. BEGIN  (* totals *)
  93.  
  94. total_x := 0;
  95. total_x2 := 0;
  96. total_x3 := 0;
  97. total_x4 := 0;
  98. sum_x := 0;
  99. FOR i := 1 TO n DO
  100.        sum_x := sum_x + x[i];
  101. mean := sum_x / n;
  102. r.mean := mean;
  103. FOR i := 1 TO n DO
  104.        BEGIN
  105.        tx := x[i] - mean;
  106.        tx2 := sqr(tx);
  107.        tx3 := tx * tx2;
  108.        tx4 := tx * tx3;
  109.        total_x := total_x + tx;
  110.        total_x2 := total_x2 + tx2;
  111.        total_x3 := total_x3 + tx3;
  112.        total_x4 := total_x4 + tx4;
  113.        END;
  114. END;  (* totals *)
  115.  
  116. BEGIN  (* jstat *)
  117.  
  118. totals;
  119. r.m1 := total_x / n;
  120. r.m2 := total_x2 / n;
  121. r.m3 := total_x3 / n;
  122. r.m4 := total_x4 / n;
  123.  
  124. r.standard_deviation := sqrt(r.m2);
  125. r.variance := r.m2;
  126. r.kurtosis := r.m4 / sqr(r.m2);
  127. r.skewness := r.m3 / sqrt(cube(r.m2));
  128. END;  (* jstat *)  .
  129.