home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / sigm / vols000 / vol082 / jstat.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1984-04-29  |  2.1 KB  |  115 lines
  1. (*  JSTAT ver 1.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. extern
  17.  
  18. type
  19. jstat_interface =
  20.     record
  21.     mean, standard_deviation,
  22.     variance, skewness, kurtosis,
  23.     m1, m2, m3, m4 : real;
  24.     end;
  25. jstat_array = array [1..1000] of real;
  26.  
  27. procedure jstat ( n : integer; var x : jstat_array;
  28.         var r : jstat_interface );
  29. var
  30. i : integer;
  31. total_x,total_x2,total_x3,total_x4 : real;
  32.  
  33. function cube ( x : real ): real;
  34. begin
  35. cube:= x * sqr(x);
  36. end;
  37.  
  38. function sqrt ( x : real ): real;
  39. var
  40. sq,a,b : real;
  41. exponent,i : integer;
  42. zap : record
  43.     case integer of
  44.     0 : (num : real);
  45.     1 : (ch8 : array [1..8] of char);
  46.     end;
  47.  
  48. begin
  49. if x = 0.0 then sqrt:=0.0 
  50. else
  51. begin
  52. sq:=abs(x);
  53.  
  54. zap.num:=sq;
  55. exponent:=ord(zap.ch8[1]);
  56. exponent:=(exponent div 2) + 32;
  57. zap.ch8[1]:=chr(exponent);
  58. a:=zap.num;
  59.  
  60. b:=0;
  61. i:=0;
  62.  
  63. while a <> b do
  64.   begin
  65.   b:=sq/a;
  66.   a:=(a+b)/2;
  67.   i:=i+1;
  68.   if i > 4 then
  69.     begin
  70.     i:=0;
  71.     if abs(a-b) < (1.0e-12 * a) then a:=b;
  72.     end;
  73.   end;
  74.  
  75. sqrt:=a;
  76. end; (* else *)
  77. end; (* sqrt *)
  78.  
  79. procedure totals;
  80. var
  81. i : integer;
  82. tx,tx2,tx3,tx4 : real;
  83. sum_x, mean : real;
  84. begin (* totals *)
  85. total_x:=0; total_x2:=0; total_x3:=0; total_x4:=0;
  86. sum_x:=0;
  87. for i:=1 to n do sum_x:=sum_x + x[i];
  88. mean:=sum_x / n;
  89. r.mean:=mean;
  90. for i:=1 to n do
  91.     begin
  92.     tx := x[i] - mean;
  93.     tx2 := sqr(tx);
  94.     tx3 := tx * tx2;
  95.     tx4 := tx * tx3;
  96.     total_x := total_x + tx;
  97.     total_x2 := total_x2 + tx2;
  98.     total_x3 := total_x3 + tx3;
  99.     total_x4 := total_x4 + tx4;
  100.     end;
  101. end; (* totals *)
  102.  
  103. begin (* jstat *)
  104. totals;
  105. r.m1 := total_x / n;
  106. r.m2 := total_x2 / n;
  107. r.m3 := total_x3 / n;
  108. r.m4 := total_x4 / n;
  109.  
  110. r.standard_deviation := sqrt(r.m2);
  111. r.variance := r.m2;
  112. r.kurtosis := r.m4 / sqr(r.m2);
  113. r.skewness := r.m3 / sqrt( cube(r.m2));
  114. end; (* jstat *).
  115.