robuststats.hpp

linux的gps应用

#pragma ident "$Id: RobustStats.hpp 341 2006-12-08 18:41:42Z btolman $"

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

/**
 * @file RobustStats.hpp
 * Namespace Robust includes basic robust statistical computations, including median,
 * median average deviation, quartiles and m-estimate, as well as implementation of
 * of stem-and-leaf plots, quantile plots and robust least squares estimation of a
 * polynomial.
 * Reference: Mason, Gunst and Hess, "Statistical Design and
 *            Analysis of Experiments," Wiley, New York, 1989.
 */
 
#ifndef GPSTK_ROBUSTSTATS_HPP
#define GPSTK_ROBUSTSTATS_HPP

//------------------------------------------------------------------------------------
// system includes
#include <string>

// GPSTk
#include "Exception.hpp"
#include "MiscMath.hpp"
#include "StringUtils.hpp"

//------------------------------------------------------------------------------------
#define ABSOLUTE(x) ((x) < T() ? -(x) : (x))
/// tuning constant used in M-estimate and Robust least squares (SRIFilter.cpp)
#define RobustTuningT (1.5)         // or 1.345
/// tuning constant used in robust estimate of variance
#define RobustTuningA (0.778)       // or 0.67
/// tuning constant used in MAD
#define RobustTuningE (0.6745)

//------------------------------------------------------------------------------------
namespace gpstk
{
   /** @addtogroup math */
   //@{

   //--------------------------------------------------------------------------------
   // quick sort, for use by robust statistics routines

   /// Comparison function for sorting.
   /// Default comparision function int comp(T a, T b) returns
   /// 1 if a > b, -1 if a < b, and 0 if a==b. A user defined comparison
   /// function may be passed as a calling argument to the sort routines.
   /// @param a and b objects of type T to be compared
   /// @return 1 if a > b, -1 if a < b, or 0 if a==b.
   template <typename T>
   int Qsort_compare(const T& a, const T& b) {
      if(a > b) return 1;
      else if(a < b) return -1;
      else return 0;
   }
   
   /// Insert sort. operator>() and operator<() must be defined for T,
   /// and a user comparison function comp(T,T) may be passed to override
   /// the default Qsort_compare().
   /// @param sa is the array of type T to be sorted.
   /// @param na length of the array to be sorted.
   template <typename T>
   void insert(T *sa,
               int na,
               int (*comp)(const T&, const T&) = Qsort_compare)
   {
      int i,j;
      T stemp;
      for(i=1; i < na; i=i+1) { // insert the i-th element into the sorted array
         stemp = sa[i];
         j = i - 1;             // find where it goes
         while((j >= 0) && (comp(stemp,sa[j]) < 0)) { 
            sa[j+1] = sa[j];
            j = j - 1;
         }
         sa[j+1] = stemp;
      }
   }  // end insert sort

   /// Quick sort in memory, with insert sort for small arrays.
   /// operator>() and operator<() must be defined for T,
   /// and a user comparison function comp(T,T) may be passed to
   /// override the default Qsort_compare().
   /// @param sa is the array of type T to be sorted.
   /// @param na length of the array to be sorted.
   /// @param comp (optional) the comparison function to be used.
   template <typename T>
   void QSort(T *sa,
              int na,
              int (*comp)(const T&, const T&) = Qsort_compare)
   {
      int i,j,nr;
      T stemp, spart;
      while(1) { 
         if(na < 8) {                     // use insert sort for small arrays
            insert(sa,na,comp);
            return;
         }
         spart = sa[na/2];                // pick middle element as pivot
         i = -1; 
         j = na;
         while(1) {
            do {                          // find first element to move right
               i = i + 1;
            } while(comp(sa[i],spart) < 0);
            do {                          // find first element to move left
               j = j - 1;
            } while(comp(sa[j], spart) > 0);
                                          // if the boundaries have met,
                                          // through paritioning,
            if(i >= j) break;
                                          // swap i and j elements
            stemp = sa[i];
            sa[i] = sa[j];
            sa[j] = stemp;
         }
         nr = na - i;
         if(i < (na/2)) {                 // now sort each partition
            QSort(sa, i, comp);           // sort left side 
            sa = &sa[i];                  // sort right side here
            na = nr;                      // memsort(&(sa[i]),nr,comp);
         }
         else { 
            QSort(&(sa[i]), nr, comp);    // sort right side
            na = i;
         }
      }
   }  // end QSort
   
   /// Insert sort one vector, keeping a second parallel.
   /// See the single-vector version of insert.
   /// @param sa is the array of type T to be sorted.
   /// @param pa is the array of type S to be kept parallel to the first.
   /// @param na length of the array to be sorted.
   template <typename T, typename S>
   void insert(T *sa,
               S *pa,
               int na,
               int (*comp)(const T&, const T&) = Qsort_compare)
   {
      int i,j;
      T stemp;
      S ptemp;
      for(i=1; i < na; i=i+1) { // insert the i-th element into the sorted array
         stemp = sa[i];
         ptemp = pa[i];
         j = i - 1;             // find where it goes
         while((j >= 0) && (comp(stemp,sa[j]) < 0)) { 
            sa[j+1] = sa[j];
            pa[j+1] = pa[j];
            j = j - 1;
         }
         sa[j+1] = stemp;
         pa[j+1] = ptemp;
      }
   }  // end insert sort

   /// Quick sort of one vector, keeping another parallel.
   /// See the single-vector version of QSort.
   /// @param sa is the array of type T to be sorted.
   /// @param pa is the array of type S to be kept parallel to the first.
   /// @param na length of the array to be sorted.
   /// @param comp (optional) the comparison function to be used.
   template <typename T, typename S>
   void QSort(T *sa,
              S *pa,
              int na,
              int (*comp)(const T&, const T&) = Qsort_compare)
   {
      int i,j,nr;
      T stemp, spart;
      S ptemp, ppart;
      while(1) { 
         if(na < 8) {                     // use insert sort for small arrays
            insert(sa,pa,na,comp);
            return;
         }
         spart = sa[na/2];                // pick middle element as pivot
         ppart = pa[na/2];
         i = -1; 
         j = na;
         while(1) {
            do {                          // find first element to move right
               i = i + 1;
            } while(comp(sa[i],spart) < 0);
            do {                          // find first element to move left
               j = j - 1;
            } while(comp(sa[j], spart) > 0);
                                          // if the boundaries have met,
                                          // through paritioning,
            if(i >= j) break;
                                          // swap i and j elements
            stemp = sa[i]; ptemp = pa[i];
            sa[i] = sa[j]; pa[i] = pa[j];
            sa[j] = stemp; pa[j] = ptemp;
         }
         nr = na - i;
         if(i < (na/2)) {                 // now sort each partition
            QSort(sa, pa, i, comp);       // sort left side 
            sa = &sa[i];                  // sort right side here
            pa = &pa[i];
            na = nr;                      // memsort(&(sa[i]),nr,comp);
         }
         else { 
            QSort(&(sa[i]), &(pa[i]), nr, comp);    // sort right side
            na = i;
         }
      }
   }  // end QSort
   
   //--------------------------------------------------------------------------------
   /// Robust statistics.
   namespace Robust
   {
   /// Compute median of an array of length nd;
   /// array xd is returned sorted, unless save_flag is true.
   /// @param xd         array of data.
   /// @param nd         length of array xd.
   /// @param save_flag  if true (default), array xd will NOT be trashed.
   /// @return median of the data in array xd.