home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snipps97.zip / MATHSTAT.C < prev    next >
C/C++ Source or Header  |  1997-07-04  |  12KB  |  588 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4. **  MATHSTAT.C - Statistical analysis in C and C++
  5. **
  6. **  Public domain by Bob Stout
  7. */
  8.  
  9. #include <stddef.h>
  10. #include <float.h>
  11. #include "mathstat.h"
  12.  
  13. /*
  14. **  Initialize a stat structure
  15. */
  16.  
  17. #if !(__cplusplus)
  18.  
  19. void stat_init(Stat_T *ptr)
  20. {
  21.       ptr->count   = 0;
  22.       ptr->total   = 0.0;
  23.       ptr->total2  = 0.0;
  24.       ptr->recip   = 0.0;
  25.       ptr->product = 1.0;
  26.       ptr->min1 = ptr->min2 = ptr->oldmin = DBL_MAX;
  27.       ptr->max1 = ptr->max2 = ptr->oldmax = DBL_MIN;
  28. }
  29.  
  30. #else /* C++ */
  31.  
  32. inline Stat::Stat()
  33. {
  34.       count   = 0;
  35.       total   = 0.0;
  36.       total2  = 0.0;
  37.       recip   = 0.0;
  38.       product = 1.0;
  39.       min1 = min2 = oldmin = DBL_MAX;
  40.       max1 = max2 = oldmax = DBL_MIN;
  41. }
  42.  
  43. #endif
  44.  
  45. /*
  46. **  Return the number of data points
  47. */
  48.  
  49. #if !(__cplusplus)
  50.  
  51. size_t stat_count(Stat_T *ptr)
  52. {
  53.       return ptr->count;
  54. }
  55.  
  56. #else /* C++ */
  57.  
  58. inline size_t Stat::Count()
  59. {
  60.       return count;
  61. }
  62.  
  63. #endif
  64.  
  65. /*
  66. **  Add a value for statistical analysis
  67. */
  68.  
  69. #if !(__cplusplus)
  70.  
  71. size_t stat_add(double datum, Stat_T *ptr)
  72. {
  73.       ptr->total   += datum;
  74.       ptr->total2  += (datum * datum);
  75.       ptr->recip   += 1.0 / datum;
  76.       ptr->product *= datum;
  77.       ++(ptr->count);
  78.       if (datum < ptr->min1)
  79.       {
  80.             ptr->oldmin = ptr->min1;
  81.             ptr->min1 = datum;
  82.       }
  83.       else if (datum < ptr->min2)
  84.             ptr->min2 = datum;
  85.       if (datum > ptr->max1)
  86.       {
  87.             ptr->oldmax = ptr->max1;
  88.             ptr->max1 = datum;
  89.       }
  90.       else if (datum < ptr->max2)
  91.             ptr->max2 = datum;
  92.       return ptr->count;
  93. }
  94.  
  95. #else /* C++ */
  96.  
  97. size_t Stat::Add(double datum)
  98. {
  99.       total   += datum;
  100.       total2  += (datum * datum);
  101.       recip   += 1.0 / datum;
  102.       product *= datum;
  103.       ++(count);
  104.       if (datum < min1)
  105.       {
  106.             oldmin = min1;
  107.             min1 = datum;
  108.       }
  109.       else if (datum < min2)
  110.             min2 = datum;
  111.       if (datum > max1)
  112.       {
  113.             oldmax = max1;
  114.             max1 = datum;
  115.       }
  116.       else if (datum < max2)
  117.             max2 = datum;
  118.       return count;
  119. }
  120.  
  121. #endif
  122.  
  123. /*
  124. **  Delete a value from a statistical analysis
  125. */
  126.  
  127. #if !(__cplusplus)
  128.  
  129. size_t stat_delete(double datum, Stat_T *ptr)
  130. {
  131.       ptr->total   -= datum;
  132.       ptr->total2  -= (datum * datum);
  133.       ptr->recip   -= 1.0 / datum;
  134.       ptr->product /= datum;
  135.       --(ptr->count);
  136.       if (datum == ptr->min1)
  137.             ptr->min1 = (DBL_MAX == ptr->oldmin) ? ptr->min2 : ptr->oldmin;
  138.       if (datum == ptr->max1)
  139.             ptr->max1 = (DBL_MIN == ptr->oldmax) ? ptr->max2 : ptr->oldmax;
  140.       return ptr->count;
  141. }
  142.  
  143. #else /* C++ */
  144.  
  145. size_t Stat::Delete(double datum)
  146. {
  147.       total   -= datum;
  148.       total2  -= (datum * datum);
  149.       recip   -= 1.0 / datum;
  150.       product /= datum;
  151.       --(count);
  152.       if (datum == min1)
  153.             min1 = (DBL_MAX == oldmin) ? min2 : oldmin;
  154.       if (datum == max1)
  155.             max1 = (DBL_MIN == oldmax) ? max2 : oldmax;
  156.       return count;
  157. }
  158.  
  159. #endif
  160.  
  161. /*
  162. **  "Olympic filter" - toss out the high and low data
  163. */
  164.  
  165. #if !(__cplusplus)
  166.  
  167. Boolean_T stat_olympic(Stat_T *ptr)
  168. {
  169.       if (ptr->count < 3)
  170.             return Error_;
  171.       stat_delete(ptr->min1, ptr);
  172.       stat_delete(ptr->max1, ptr);
  173.       return Success_;
  174. }
  175.  
  176. #else /* C++ */
  177.  
  178. Boolean_T Stat::Olympic()
  179. {
  180.       if (count < 3)
  181.             return Error_;
  182.       Delete(min1);
  183.       Delete(max1);
  184.       return Success_;
  185. }
  186.  
  187. #endif
  188.  
  189. /*
  190. **  Return the minimum datum
  191. */
  192.  
  193. #if !(__cplusplus)
  194.  
  195. double stat_min(Stat_T *ptr)
  196. {
  197.       return ptr->min1;
  198. }
  199.  
  200. #else /* C++ */
  201.  
  202. inline double Stat::Min()
  203. {
  204.       return min1;
  205. }
  206.  
  207. #endif
  208.  
  209. /*
  210. **  Return the maximum datum
  211. */
  212.  
  213. #if !(__cplusplus)
  214.  
  215. double stat_max(Stat_T *ptr)
  216. {
  217.       return ptr->max1;
  218. }
  219.  
  220. #else /* C++ */
  221.  
  222. inline double Stat::Max()
  223. {
  224.       return max1;
  225. }
  226.  
  227. #endif
  228.  
  229. /*
  230. **  Return the error (%) for the minimum datum
  231. */
  232.  
  233. #if !(__cplusplus)
  234.  
  235. double stat_minerror(Stat_T *ptr)
  236. {
  237.       double Mean = stat_mean(ptr);
  238.  
  239.       return 100.0 * ((ptr->min1 - Mean) / Mean) ;
  240. }
  241.  
  242. #else /* C++ */
  243.  
  244. double Stat::Minerror()
  245. {
  246.       double mean_ = Mean();
  247.  
  248.       return 100.0 * ((min1 - mean_) / mean_) ;
  249. }
  250.  
  251. #endif
  252.  
  253. /*
  254. **  Return the error (%) for the maximum datum
  255. */
  256.  
  257. #if !(__cplusplus)
  258.  
  259. double stat_maxerror(Stat_T *ptr)
  260. {
  261.       double Mean = stat_mean(ptr);
  262.  
  263.       return 100.0 * ((ptr->max1 - Mean) / Mean) ;
  264. }
  265.  
  266. #else /* C++ */
  267.  
  268. double Stat::Maxerror()
  269. {
  270.       double mean_ = Mean();
  271.  
  272.       return 100.0 * ((max1 - mean_) / mean_) ;
  273. }
  274.  
  275. #endif
  276.  
  277. /*
  278. **  Compute the arithmetic mean
  279. */
  280.  
  281. #if !(__cplusplus)
  282.  
  283. double stat_mean(Stat_T *ptr)
  284. {
  285.       if (ptr->count)
  286.             return (ptr->total / ptr->count);
  287.       else  return 0.0;
  288. }
  289.  
  290. #else /* C++ */
  291.  
  292. double Stat::Mean()
  293. {
  294.       if (count)
  295.             return (total / count);
  296.       else  return 0.0;
  297. }
  298.  
  299. #endif
  300.  
  301. /*
  302. **  Compute the geometric mean
  303. */
  304.  
  305. #if !(__cplusplus)
  306.  
  307. double stat_gmean(Stat_T *ptr)
  308. {
  309.       if (ptr->count)
  310.             return pow(ptr->product, 1.0 / ptr->count);
  311.       else  return 0.0;
  312. }
  313.  
  314. #else /* C++ */
  315.  
  316. double Stat::Gmean()
  317. {
  318.       if (count)
  319.             return pow(product, 1.0 / count);
  320.       else  return 0.0;
  321. }
  322.  
  323. #endif
  324.  
  325. /*
  326. **  Compute the harmonic mean
  327. */
  328.  
  329. #if !(__cplusplus)
  330.  
  331. double stat_hmean(Stat_T *ptr)
  332. {
  333.       if (ptr->count)
  334.             return ptr->count / ptr->recip;
  335.       else  return 0.0;
  336. }
  337.  
  338. #else /* C++ */
  339.  
  340. double Stat::Hmean()
  341. {
  342.       if (count)
  343.             return count / recip;
  344.       else  return 0.0;
  345. }
  346.  
  347. #endif
  348.  
  349. /*
  350. **  Compute the standard deviation of the population
  351. */
  352.  
  353. #if !(__cplusplus)
  354.  
  355. double stat_stddevP(Stat_T *ptr)
  356. {
  357.       double tmp;
  358.  
  359.       tmp  = ptr->total / ptr->count;
  360.       tmp *= tmp;
  361.       return sqrt((ptr->total2 / ptr->count) - tmp);
  362. }
  363.  
  364. #else /* C++ */
  365.  
  366. double Stat::StddevP()
  367. {
  368.       double tmp;
  369.  
  370.       tmp  = total / count;
  371.       tmp *= tmp;
  372.       return sqrt((total2 / count) - tmp);
  373. }
  374.  
  375. #endif
  376.  
  377. /*
  378. **  Compute the standard deviation of the sampled data
  379. */
  380.  
  381. #if !(__cplusplus)
  382.  
  383. double stat_stddevS(Stat_T *ptr)
  384. {
  385.       double mean, total, total2, N, tmp;
  386.  
  387.       if (ptr->count < 2)
  388.             return 0.0;
  389.  
  390.       mean   = ptr->total / ptr->count;
  391.       total  = ptr->total - mean;
  392.       total2 = ptr->total2 - (mean * mean);
  393.       N      = ptr->count - 1;
  394.       tmp    = total / N;
  395.       tmp   *= tmp;
  396.       return sqrt((total2 / N) - tmp);
  397. }
  398.  
  399. #else /* C++ */
  400.  
  401. double Stat::StddevS()
  402. {
  403.       double Mean, Total, Total2, N, tmp;
  404.  
  405.       if (count < 2)
  406.             return 0.0;
  407.  
  408.       Mean   = total / count;
  409.       Total  = total - Mean;
  410.       Total2 = total2 - (Mean * Mean);
  411.       N      = count - 1;
  412.       tmp    = Total / N;
  413.       tmp   *= tmp;
  414.       return sqrt((Total2 / N) - tmp);
  415. }
  416.  
  417. #endif
  418.  
  419. /*
  420. **  Compute the variance
  421. */
  422.  
  423. #if !(__cplusplus)
  424.  
  425. double stat_var(Stat_T *ptr)
  426. {
  427.       return stat_stddevS(ptr) * stat_stddevS(ptr);
  428. }
  429.  
  430. #else /* C++ */
  431.  
  432. inline double Stat::Var()
  433. {
  434.       return StddevS() * StddevS();
  435. }
  436.  
  437. #endif
  438.  
  439. /*
  440. **  Compute the coefficient of variation (percentage) for the population
  441. */
  442.  
  443. #if !(__cplusplus)
  444.  
  445. double stat_varcoeffP(Stat_T *ptr)
  446. {
  447.       return (stat_stddevP(ptr) / stat_mean(ptr)) * 100.0;
  448. }
  449.  
  450. #else /* C++ */
  451.  
  452. double Stat::VarcoeffP()
  453. {
  454.       return (StddevP() / Mean()) * 100.0;
  455. }
  456.  
  457. #endif
  458.  
  459. /*
  460. **  Compute the coefficient of variation (percentage) for the sample
  461. */
  462.  
  463. #if !(__cplusplus)
  464.  
  465. double stat_varcoeffS(Stat_T *ptr)
  466. {
  467.       return (stat_stddevS(ptr) / stat_mean(ptr)) * 100.0;
  468. }
  469.  
  470. #else /* C++ */
  471.  
  472. double Stat::VarcoeffS()
  473. {
  474.       return (StddevS() / Mean()) * 100.0;
  475. }
  476.  
  477. #endif
  478.  
  479.  
  480. /*
  481. **  Test harness begins here. Define TEST macro to build stand-alone.
  482. */
  483.  
  484. #ifdef TEST
  485.  
  486. #include <stdlib.h>
  487.  
  488. #if !(__cplusplus)
  489.  
  490. #include <stdio.h>
  491.  
  492. int main(int argc, char *argv[])
  493. {
  494.       Stat_T data;
  495.  
  496.       stat_init(&data);
  497.       while (--argc)
  498.       {
  499.             double ftmp;
  500.  
  501.             ftmp = atof(*(++argv));
  502.             stat_add(ftmp, &data);
  503.       }
  504.       puts("\nBefore \"Olympic\" filtering\n");
  505.       printf("Minimum datum             = %g\n", stat_min(&data));
  506.       printf("Maximum datum             = %g\n", stat_max(&data));
  507.       printf("Number of samples         = %d\n", stat_count(&data));
  508.       printf("Arithmetic mean           = %g\n", stat_mean(&data));
  509.       printf("Geometric mean            = %g\n", stat_gmean(&data));
  510.       printf("Harmonic mean             = %g\n", stat_hmean(&data));
  511.       printf("Standard deviation (N)    = %g\n", stat_stddevP(&data));
  512.       printf("Standard deviation (N-1)  = %g\n", stat_stddevS(&data));
  513.       printf("Variance                  = %g\n", stat_var(&data));
  514.       printf("Population coeff. of var. = %g%%\n", stat_varcoeffP(&data));
  515.       printf("Sample coeff. of var.     = %g%%\n", stat_varcoeffS(&data));
  516.  
  517.       puts("\nAfter \"Olympic\" filtering\n");
  518.       printf("stat_olympic() returned %s\n", stat_olympic(&data) ?
  519.             "ERROR" : "SUCCESS");
  520.       printf("Minimum datum             = %g\n", stat_min(&data));
  521.       printf("Maximum datum             = %g\n", stat_max(&data));
  522.       printf("Number of samples         = %d\n", stat_count(&data));
  523.       printf("Arithmetic mean           = %g\n", stat_mean(&data));
  524.       printf("Geometric mean            = %g\n", stat_gmean(&data));
  525.       printf("Harmonic mean             = %g\n", stat_hmean(&data));
  526.       printf("Standard deviation (N)    = %g\n", stat_stddevP(&data));
  527.       printf("Standard deviation (N-1)  = %g\n", stat_stddevS(&data));
  528.       printf("Variance                  = %g\n", stat_var(&data));
  529.       printf("Population coeff. of var. = %g%%\n", stat_varcoeffP(&data));
  530.       printf("Sample coeff. of var.     = %g%%\n", stat_varcoeffS(&data));
  531.  
  532.       return EXIT_SUCCESS;
  533. }
  534.  
  535. #else /* C++ */
  536.  
  537. #include <iostreams.h>
  538.  
  539. int main(int argc, char *argv[])
  540. {
  541.       class Stat data;
  542.       char *str;
  543.  
  544.       while (--argc)
  545.       {
  546.             double ftmp;
  547.  
  548.             ftmp = atof(*(++argv));
  549.             data.Add(ftmp);
  550.       }
  551.       cout << endl << "Before \"Olympic\" filtering\n" << endl << endl;
  552.       cout << "Minimum datum             = " << data.Min()       << endl;
  553.       cout << "Maximum datum             = " << data.Max()       << endl;
  554.       cout << "Number of samples         = " << data.Count()     << endl;
  555.       cout << "Arithmetic mean           = " << data.Mean()      << endl;
  556.       cout << "Geometric mean            = " << data.Gmean()     << endl;
  557.       cout << "Harmonic mean             = " << data.Hmean()     << endl;
  558.       cout << "Standard deviation (N)    = " << data.StddevP()   << endl;
  559.       cout << "Standard deviation (N-1)  = " << data.StddevS()   << endl;
  560.       cout << "Variance                  = " << data.Var()       << endl;
  561.       cout << "Population coeff. of var. = " << data.VarcoeffP() << "%" << endl;
  562.       cout << "Sample coeff. of var.     = " << data.VarcoeffS() << "%" << endl;
  563.  
  564.       if (Success_ == data.Olympic())
  565.             str = "SUCCESS";
  566.       else  str = "ERROR";
  567.  
  568.       cout << endl << "After \"Olympic\" filtering" << endl << endl;
  569.       cout << "data.Olympic() returned " << str << endl;
  570.       cout << "Minimum datum             = " << data.Min()       << endl;
  571.       cout << "Maximum datum             = " << data.Max()       << endl;
  572.       cout << "Number of samples         = " << data.Count()     << endl;
  573.       cout << "Arithmetic mean           = " << data.Mean()      << endl;
  574.       cout << "Geometric mean            = " << data.Gmean()     << endl;
  575.       cout << "Harmonic mean             = " << data.Hmean()     << endl;
  576.       cout << "Standard deviation (N)    = " << data.StddevP()   << endl;
  577.       cout << "Standard deviation (N-1)  = " << data.StddevS()   << endl;
  578.       cout << "Variance                  = " << data.Var()       << endl;
  579.       cout << "Population coeff. of var. = " << data.VarcoeffP() << "%" << endl;
  580.       cout << "Sample coeff. of var.     = " << data.VarcoeffS() << "%" << endl;
  581.  
  582.       return EXIT_SUCCESS;
  583. }
  584.  
  585. #endif /* C/C++ */
  586.  
  587. #endif
  588.