mathstat.c

国外网站上的一些精典的C程序

**  Returns: Coefficient of variation of the sampled data
**
**  Note: If no data, return STAT_NaN
*/

#if !(__cplusplus)

double stat_varcoeffS(Stat_T *ptr)
{
      if (ptr->count >= 2)
            return (stat_stddevS(ptr) / stat_mean(ptr)) * 100.0;
      else  return STAT_NaN;
}

#else /* C++ */

double Stat::VarcoeffS()
{
      if (ptr->count >= 2)
            return (StddevS() / Mean()) * 100.0;
      else  return STAT_NaN;
}

#endif

/*
**  Notes about filters:
**
**  Sampled data are always subject to spurious random effects from noise or simple
**  sampling errors. To improve the quality of the data population, you use a filter.
**  A filter is simply a device which removes spurious data. After running a filter,
**  the quantity of your data may be reduced, but its quality is improved.
**
**  So it is with the filters here. Each acts on the stored data to selectively remove
**  specific data with the goal of reducing the variance of the remaining population.
**  After running the filter, new statistics may be run on the remaining data set, the
**  most common of which are arithmetic mean and standard deviation/variance.
**
**  Three specific filter types are currently implemented, although more are planned:
**
**  1.  An "Olympic" filter derives its nickname from a scoring technique used in the
**      Olympic games where the high and low data are removed from consideration. This
**      is a time-tested and amazingly effective technique which is especially well
**      adapted to resource-restricted embedded systems.
**
**  2.  A median filter simply removes all but a fixed-width window of data surrounding
**      the mid-point of an array of the sampled data. For example, if you have a
**      population of 10 samples and impose a median filter with a window width of 4,
**      the lowest 3 and highest 3 data would be removed.
**
**  3.  Combining elements of both the above, an "Olympiad" filter is an iterative
**      Olympic filter, with  smart windowing. You call an Olympiad filter with a
**      target value for the maximum tolerable coefficient of variation (expressed as
**      a percentage). The Olympiad filter will iteratively remove high and low data
**      until either the target is achieved or no more data may be removed. The only
**      disadvantage to the Olympiad filter is that by using the coefficient of
**      variation, it inherently uses the arithmetic mean of the data. Should you wish
**      to filter using the geometric or harmonic mean, you would have to create your
**      own function which iteratively called the Olympic filter.
**
**
**  Implementation notes:
**
**  As written, these functions are meant to be general purpose. As such, they use
**  techniques which might not be suitable for specific environments. In particular,
**  many embedded systems cannot afford to use dynamic memory allocation. It is left
**  as an exercise to the reader to adapt these as required for a specific application.
**
**  The original Olympic filter module used fixed registers for min and max data to
**  date. This is extremely fast and efficient, but inadequate for more advanced
**  filtering techniques. The current scheme saves all data entered into the analysis
**  in sorted order, so the min and max are always available at the top and bottom of
**  the array. The price is that the data array must be maintained for each datum
**  entry and deletion. To return to the original system, you must replace the data
**  array with min and max registers. Without the data array, only the Olympic filter
**  should be used and even it can no longer be used iteratively.
**
**  The min and max data are still saved if no data array is defined. However, an
**  analysis with a NULL data array pointer will not support filtering or removal of
**  data from the anlysis.
*/

/*** Begin filters ******************************************************************/

/*
**  "Olympic filter" - toss out the high and low data
**
**  Parameters: 1 - Pointer to Stat_T structure type
**
**  Returns: Error_ if less than 3 data, else Success_
*/

#if !(__cplusplus)

Boolean_T stat_olympic(Stat_T *ptr)
{
      Boolean_T retval;

      if (!ptr->filt || ptr->count < 3)
            return Error_;

      stat_delete(ptr->filt[0], ptr);
      stat_delete(ptr->filt[ptr->count - 1], ptr);

      return Success_;
}

#else /* C++ */

Boolean_T Stat::Olympic()
{
      if (!filt || count < 3)
            return Error_;

      Delete(filt[0]);
      Delete(filt[count - 1]);

      return Success_;
}

#endif

/*
**  Iterative "Olympic" filter - toss out the high and low data to reach a desired
**  coefficient of variation for the population. Functionally, it may be considered
**  as a "smart" median filter which auto-adjusts its window.
**
**  Parameters: 1 - Pointer to Stat_T structure type
**              2 - Target coefficient of variation
**
**  Returns: Error_ if target is not achieved, else Success_
**
**  Note: Since this filter uses the coefficient of variation, it uses the arithmetic
**        mean for its target. If you want to use either the geometric or harmonic
**        mean instead, you can do it by explicitly iterating calls to the Opympic
**        filter instead.
*/

#if !(__cplusplus)

Boolean_T stat_olympiad(Stat_T *ptr, double CoefVar)
{
      double v;

      do
      {
            if (!ptr->filt || ptr->count < 3)
                  return Error_;

            stat_delete(ptr->filt[0], ptr);
            stat_delete(ptr->filt[ptr->count - 1], ptr);

            v = stat_varcoeffP(ptr);
      } while (v > CoefVar);

      return Success_;
}

#else /* C++ */

Boolean_T Stat::Olympiad(double CoefVar)
{
      double v;

      do
      {
            if (!filt || count < 3)
                  return Error_;

            Delete(filt[0]);
            Delete(filt[count - 1]);

            v = VarcoeffP();
      } while (v > CoefVar);

      return Success_;
}

#endif

/*
**  Median filter
**
**  Parameters: 1 - Pointer to Stat_T structure type
**
**  Returns: Nothing
**
**  Side effects: Removes all data outside the specified window from the data array
*/

#if !(__cplusplus)

Boolean_T stat_median(Stat_T *ptr, size_t window)
{
      size_t i, n;

      if (!ptr->filt || (ptr->count < 1) || (ptr->count < window))
            return Error_;

      for (i = 0, n = (ptr->count - window) / 2; i < n; ++i)
            stat_delete(ptr->filt[0], ptr);

      for (i = window; ptr->count >= i; ++i)
            stat_delete(ptr->filt[window], ptr);

      return Success_;
}

#else /* C++ */

Boolean_T Stat::Median(Stat_T *ptr, size_t window)
{
      size_t i, n;

      if (!filt || (count < 1) || (count < window))
            return Error_;

      for (i = 0, n = (count - window) / 2; i < n; ++i)
            stat_delete(filt[0], ptr);

      for (i = window; count >= i; ++i)
            stat_delete(filt[window], ptr);

      return Success_;
}

#endif

/*** End filters ********************************************************************/

/*
**  Debugging functions
*/

/*
**  Debugging function to dump the Stat_T registers
*/

#if !(__cplusplus)

void stat_dump(Stat_T *ptr)
{
      printf("stat_dump(): count   = %d\n", ptr->count);
      printf("             total   = %d\n", ptr->total);
      printf("             total2  = %f\n", ptr->total2);
      printf("             recip   = %f\n", ptr->recip);
      printf("             product = %f\n", ptr->product);
      printf("             filt    = %p\n", ptr->filt);
      printf("             min     = %p\n", ptr->min);
      printf("             max     = %p\n", ptr->max);
}

#else /* C++ */

void Stat::dump()
{
      cout << "dump(): count   = " << count << endl;
      cout << "        total   = " << total << endl;
      cout << "        total2  = " << total2 << endl;
      cout << "        recip   = " << recip << endl;
      cout << "        product = " << product << endl;
      cout << "        filt    = " << filt << endl;
      cout << "        min     = " << min << endl;
      cout << "        max     = " << min << endl;
}

#endif

/*
**  Debugging function to dump the data array
*/

#if !(__cplusplus)

void stat_qdump(Stat_T *ptr)
{
      size_t i;

      puts(">>DATA ARRAY");
      for (i = 0; i < ptr->count; ++i)
            printf("filt[%02d]=%f\n", i, ptr->filt[i]);
}

#else /* C++ */

void Stat::qdump()
{
      size_t i;

      cout << ">>DATA ARRAY" << endl;
      for (i = 0; i < ptr->count; ++i)
            cout << "filt[" << i << "]=" << i, ptr->filt[i] << endl;
}

#endif

/*
**  Test harness begins here. Define TEST macro to build stand-alone.
*/

#ifdef TEST

#include <stdlib.h>

#if !(__cplusplus)

#include <stdio.h>

int main(int argc, char *argv[])
{
      Stat_T data;
      size_t i;

      stat_init(&data, 20);
      while (--argc)
      {
            double ftmp;

            ftmp = atof(*(++argv));
            stat_add(ftmp, &data);
      }

      puts("\nBefore \"Olympic\" filtering\n");
      printf("Minimum datum             = %g\n", stat_min(&data));
      printf("Maximum datum             = %g\n", stat_max(&data));
      printf("Number of samples         = %d\n", stat_count(&data));
      printf("Arithmetic mean           = %g\n", stat_mean(&data));
      printf("Geometric mean            = %g\n", stat_gmean(&data));
      printf("Harmonic mean             = %g\n", stat_hmean(&data));
      printf("Standard deviation (N)    = %g\n", stat_stddevP(&data));
      printf("Standard deviation (N-1)  = %g\n", stat_stddevS(&data));
      printf("Variance (population)     = %g\n", stat_varP(&data));
      printf("Variance (sample)         = %g\n", stat_varS(&data));
      printf("Population coeff. of var. = %g%%\n", stat_varcoeffP(&data));
      printf("Sample coeff. of var.     = %g%%\n", stat_varcoeffS(&data));
/*      qdump(&data);*/

      i = data.count - 2;
      if (Success_ == stat_median(&data, i))
            printf("\nMedian filter (window=%d)  = %g\n", i, stat_mean(&data));
      else  printf("Median filter failed, window=%d\n", i);
/*      qdump(&data);*/

      printf("\nDeleting sample #%d\n", data.count / 2);
      stat_delete(data.filt[data.count / 2], &data);
/*      qdump(&data);*/

      puts("\nAfter \"Olympic\" filtering\n");
      printf("stat_olympic() returned %s\n", stat_olympic(&data) ?
            "ERROR" : "SUCCESS");
      printf("Minimum datum             = %g\n", stat_min(&data));
      printf("Maximum datum             = %g\n", stat_max(&data));
      printf("Number of samples         = %d\n", stat_count(&data));
      printf("Arithmetic mean           = %g\n", stat_mean(&data));
      printf("Geometric mean            = %g\n", stat_gmean(&data));
      printf("Harmonic mean             = %g\n", stat_hmean(&data));
      printf("Standard deviation (N)    = %g\n", stat_stddevP(&data));
      printf("Standard deviation (N-1)  = %g\n", stat_stddevS(&data));
      printf("Variance (population)     = %g\n", stat_varP(&data));
      printf("Variance (sample)         = %g\n", stat_varS(&data));
      printf("Population coeff. of var. = %g%%\n", stat_varcoeffP(&data));
      printf("Sample coeff. of var.     = %g%%\n", stat_varcoeffS(&data));
/*      qdump(&data);*/

      puts("\nAfter \"Olympiad\" filtering, CoeffVar = 5.0\n");
      printf("stat_olympiad() returned %s\n", stat_olympiad(&data, 5.0) ?
             "ERROR" : "SUCCESS");
      printf("Minimum datum             = %g\n", stat_min(&data));
      printf("Maximum datum             = %g\n", stat_max(&data));
      printf("Number of samples         = %d\n", stat_count(&data));
      printf("Arithmetic mean           = %g\n", stat_mean(&data));
      printf("Geometric mean            = %g\n", stat_gmean(&data));
      printf("Harmonic mean             = %g\n", stat_hmean(&data));
      printf("Standard deviation (N)    = %g\n", stat_stddevP(&data));
      printf("Standard deviation (N-1)  = %g\n", stat_stddevS(&data));
      printf("Variance (population)     = %g\n", stat_varP(&data));
      printf("Variance (sample)         = %g\n", stat_varS(&data));
      printf("Population coeff. of var. = %g%%\n", stat_varcoeffP(&data));
      printf("Sample coeff. of var.     = %g%%\n", stat_varcoeffS(&data));
/*      qdump(&data);*/

      return EXIT_SUCCESS;
}

#else /* C++ */

#include <iostreams.h>

int main(int argc, char *argv[])
{
      class Stat data(20);
      char *str;

      while (--argc)
      {
            double ftmp;

            ftmp = atof(*(++argv));
            data.Add(ftmp);
      }

      cout << endl << "Before \"Olympic\" filtering\n" << endl << endl;
      cout << "Minimum datum             = " << data.Min()       << endl;
      cout << "Maximum datum             = " << data.Max()       << endl;
      cout << "Number of samples         = " << data.Count()     << endl;
      cout << "Arithmetic mean           = " << data.Mean()      << endl;
      cout << "Geometric mean            = " << data.Gmean()     << endl;
      cout << "Harmonic mean             = " << data.Hmean()     << endl;
      cout << "Standard deviation (N)    = " << data.StddevP()   << endl;
      cout << "Standard deviation (N-1)  = " << data.StddevS()   << endl;
      cout << "Variance (population)     = " << data.VarP()       << endl;
      cout << "Variance (sample)         = " << data.VarS()       << endl;
      cout << "Population coeff. of var. = " << data.VarcoeffP() << "%" << endl;
      cout << "Sample coeff. of var.     = " << data.VarcoeffS() << "%" << endl;

      if (Success_ == data.Olympic())
            str = "SUCCESS";
      else  str = "ERROR";

      cout << endl << "After \"Olympic\" filtering" << endl << endl;
      cout << "data.Olympic() returned " << str << endl;
      cout << "Minimum datum             = " << data.Min()       << endl;
      cout << "Maximum datum             = " << data.Max()       << endl;
      cout << "Number of samples         = " << data.Count()     << endl;
      cout << "Arithmetic mean           = " << data.Mean()      << endl;
      cout << "Geometric mean            = " << data.Gmean()     << endl;
      cout << "Harmonic mean             = " << data.Hmean()     << endl;
      cout << "Standard deviation (N)    = " << data.StddevP()   << endl;
      cout << "Standard deviation (N-1)  = " << data.StddevS()   << endl;
      cout << "Variance (population)     = " << data.VarP()       << endl;
      cout << "Variance (sample)         = " << data.VarS()       << endl;
      cout << "Population coeff. of var. = " << data.VarcoeffP() << "%" << endl;
      cout << "Sample coeff. of var.     = " << data.VarcoeffS() << "%" << endl;

      return EXIT_SUCCESS;
}

#endif /* C/C++ */

#endif