robuststats.hpp

1 gps基本算法

   T Median(T *xd, const int nd, bool save_flag=true)
      throw(Exception)
   {
      if(!xd || nd < 2) {
         Exception e("Invalid input");
         GPSTK_THROW(e);
      }

      try {
         int i;
         T med, *save;

         if(save_flag) {
            save = new T[nd];
            if(!save) {
               Exception e("Could not allocate temporary array");
               GPSTK_THROW(e);
            }
            for(i=0; i<nd; i++) save[i]=xd[i];
         }

         QSort(xd,nd);

         if(nd%2)
            med = xd[(nd+1)/2-1];
         else
            med = (xd[nd/2-1]+xd[nd/2])/T(2);

            // restore original data from temporary
         if(save_flag) {
            for(i=0; i<nd; i++) xd[i]=save[i];
            delete[] save;
         }

         return med;
      }
      catch(Exception& e) { GPSTK_RETHROW(e); }

   }  // end Median

   /// Compute the quartiles Q1 and Q3 of an array of length nd.
   /// Array is assumed sorted in ascending order.
   /// Quartile are values such that one fourth of the
   /// samples are larger (smaller) than Q3(Q1).
   /// @param xd array of data.
   /// @param nd length of array xd.
   /// @param Q1 (output) first quartile of data in array xd.
   /// @param Q3 (output) third quartile of data in array xd.
   template <typename T>
   void Quartiles(const T *xd, const int nd, T& Q1, T& Q3)
      throw(Exception)
   {
      if(!xd || nd < 2) {
         Exception e("Invalid input");
         GPSTK_THROW(e);
      }

      int q;
      if(nd % 2) q = (nd+1)/2;
      else       q = nd/2;

      if(q % 2) {
         Q1 = xd[(q+1)/2-1];
         Q3 = xd[nd-(q+1)/2];
      }
      else {
         Q1 = (xd[q/2-1]+xd[q/2])/T(2);
         Q3 = (xd[nd-q/2]+xd[nd-q/2-1])/T(2);
      }
   }  // end Quartiles

   /// Compute the median absolute deviation of a double array of length nd,
   /// as well as the median (M = Median(xd,nd));
   /// NB this routine will trash the array xd unless save_flag is true (default).
   /// @param xd array of data (input).
   /// @param nd length of array xd (input).
   /// @param M median of data in array xd (output).
   /// @param save_flag if true (default), array xd will NOT be trashed.
   /// @return median absolute deviation of data in array xd.
   template <typename T>
   T MedianAbsoluteDeviation(T *xd, int nd, T& M, bool save_flag=true)
      throw(Exception)
   {
      int i;
      T mad, *save;

      if(!xd || nd < 2) {
         Exception e("Invalid input");
         GPSTK_THROW(e);
      }

         // store data in a temporary array
      if(save_flag) {
         save = new T[nd];
         if(!save) {
            Exception e("Could not allocate temporary array");
            GPSTK_THROW(e);
         }
         for(i=0; i<nd; i++) save[i]=xd[i];
      }

         // get the median (don't care if xd gets sorted...)
      M = Median(xd, nd, false);

         // compute xd=abs(xd-M)
      for(i=0; i<nd; i++) xd[i] = ABSOLUTE(xd[i]-M);

         // sort xd in ascending order
      QSort(xd, nd);

         // find median and normalize to get mad
      mad = Median(xd, nd, false) / T(RobustTuningE);

         // restore original data from temporary
      if(save_flag) {
         for(i=0; i<nd; i++) xd[i]=save[i];
         delete[] save;
      }

      return mad;

   }  // end MedianAbsoluteDeviation

      /// Compute the median absolute deviation of a double array of length nd;
      /// see MedianAbsoluteDeviation().
   template <typename T>
   T MAD(T *xd, int nd, T& M, bool save_flag=true)
      throw(Exception)
   { return MedianAbsoluteDeviation(xd,nd,M,save_flag); }

   /// Compute the m-estimate. Iteratively determine the m-estimate, which
   /// is a measure of mean or median, but is less sensitive to outliers.
   /// M is the median (M=Median(xd,nd)), and MAD is the
   /// median absolute deviation (MAD=MedianAbsoluteDeviation(xd,nd,M)).
   /// Optionally, a pointer to an array, which will contain the weights
   /// on output, may be provided.
   /// @param xd input array of data.
   /// @param nd input length of array xd.
   /// @param M input median of data in array xd.
   /// @param MAD input median absolute deviation of data in array xd.
   /// @param w output array of length nd to contain weights on output.
   /// @return m-estimate of data in array xd.
   template <typename T>
   T MEstimate(T *xd, int nd, const T& M, const T& MAD, T *w=NULL)
      throw(Exception)
   {
      try {
         T tv, m, mold, sum, sumw, *wt, weight, *t;
         T tol=0.000001;
         int i, n, N=10;      // N is the max number of iterations

         if(!xd || nd < 2) {
            Exception e("Invalid input");
            GPSTK_THROW(e);
         }

         tv = T(RobustTuningT)*MAD;
         n = 0;
         m = M;
         do {
            mold = m;
            n++;
            sum = sumw = T();
            for(i=0; i<nd; i++) {
               if(w) wt=&w[i];
               else wt=&weight;
               *wt = T(1);
               if(xd[i]<m-tv)      *wt = -tv/(xd[i]-m);
               else if(xd[i]>m+tv) *wt =  tv/(xd[i]-m);
               sumw += (*wt);
               sum += (*wt)*xd[i];
            }
            m = sum / sumw;

         } while(T(ABSOLUTE((m-mold)/m)) > tol && n < N);

         return m;
      }
      catch(Exception& e) { GPSTK_RETHROW(e); }

   }  // end MEstimate

   /// Fit a polynomial of degree n to data xd, with independent variable td,
   /// using robust techniques. The post-fit residuals are returned in the data
   /// vector, and the computed weights in the result may be output as well.
   /// Specifically, the equation describing the fit is
   /// c0 + c[1]*t(j) + c[2]*t(j)*t(j) + ... c[n-1]*pow(t(j),n-1) = xd[j],
   /// where the zero-th coefficient and the independent variable are debiased
   /// by the first value; i.e. c0 = c[0]+xd[0] and t(j) = td[i]-td[0].
   /// Specifically, to evaluate the polynomial at t, eval = f(t), do the following.
   /// xd0 = xd[0];
   /// Robust::RobustPolyFit(xd,td,nd,n,c);
   /// eval = xd0+c[0]; tt = 1.0;
   /// for(j=1; j<nd; j++) { tt *= (t-td[0]); eval += c[j]*tt; }
   /// @param xd (input) array of data, of length nd; contains residuals on output.
   /// @param td (input) array of independent variable, length nd (parallel to xd).
   /// @param nd (input) length of arrays xd and td.
   /// @param n (input) degree of polynomial and dimension of coefficient array.
   /// @param c (output) array of coefficients (dimension n).
   /// @param w (output, if non-null) array of length nd to contain weights.
   /// @return 0 for success, -1 for singular problem, -2 failure to converge.
   int RobustPolyFit(double *xd, const double *td, int nd,
                     int n, double *c, double *w=NULL)
      throw(Exception);

   /// Print 'stem and leaf' plot of the data in the double array xd of length nd,
   /// with an optional message, on the given output ostream. It is assumed that
   /// the input array xd is sorted in ascending order.
   /// @param os ostream on which to write output.
   /// @param xd array of data.
   /// @param nd length of array xd.
   /// @param msg string containing optional message for output.
   void StemLeafPlot(std::ostream& os, double *xd, long nd,
         std::string msg=std::string(""))
      throw(Exception);

   /// Generate data for a quantile-quantile plot. Given an array of data yd,
   /// of length nd (sorted in ascending order), and another array xd of the
   /// same length, fill the xd array with data such that (xd,yd) give a
   /// quantile-quantile plot. The distribution of yd is a normal distribution
   /// to the extent that this plot is a straight line, with y-intercept and
   /// slope identified with mean and standard deviation, respectively, of the
   /// distribution.
   /// @param yd array of data, sorted in ascending order.
   /// @param nd length of array xd.
   /// @param xd array of length nd containing quantiles of yd on output.
   void QuantilePlot(double *yd, long nd, double *xd)
      throw(Exception);

   }  // end Robust namespace

   //@}

}  // end namespace gpstk

#undef ABSOLUTE

#endif