stats.hpp

一个gps小工具包

#pragma ident "$Id: Stats.hpp 208 2006-10-11 16:35:25Z btolman $"

/**
 * @file Stats.hpp
 * One and two-sample statistics
 */
 
#ifndef INCLUDE_GPSTK_STATS_HPP
#define INCLUDE_GPSTK_STATS_HPP

//============================================================================
//
//  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
//
//============================================================================

#include "MiscMath.hpp"
#include "Vector.hpp"
#include "Exception.hpp"

namespace gpstk
{
   /** @addtogroup math */
   //@{
 
/** Conventional statistics for one sample.  Constructor does the same as
 * Reset(); use this when starting a new series of input samples.
 * Results are available at any time by calling N(), Minimum(), Maximum(),
 * Average(), Variance() and StdDev().
 */
   template <class T>

   class Stats
   {
   public:
      /// constructor
      explicit Stats() { n=0; weighted=false; }

      /// constructor given a vector of data
      Stats(Vector<T>& X, Vector<T>& W)
      {
         n = 0;
         weighted = false;
         Add(X,W);
      }

      /// reset, i.e. ignore earlier data and restart sampling
      inline void Reset(void) { n=0; weighted=false; W=T(); }

      /// access the sample size
      inline unsigned int N(void) const { return n; }

      /// return minimum value
      inline T Minimum(void) const { if(n) return min; else return T(); }

      /// return maximum value
      inline T Maximum(void) const { if(n) return max; else return T(); }

      /// return computed average
      inline T Average(void) const
      {
         // normalization constant is (W=sum wts)/n, -> 1 when all wts=1
         if(n == 0)
            return T();
         if(weighted)
            return (T(n)*ave/W);
         return ave;
      }

      /// return computed variance
      inline T Variance(void) const
      {
         if(n == 0)
            return T();
         if(weighted) {
            T wn=W/T(n);
            return (var/wn/wn/wn/wn);
         }
         return var;
      }

      /// return computed standard deviation
      inline T StdDev(void) const
      {
         if(n == 0)
            return T();
         if(weighted) {
            T wn=W/T(n);
            return SQRT(ABS(var))/wn/wn;
         }
         return SQRT(ABS(var));
      }

      /// return the normalization constant = sum weights
      inline T Normalization(void) const
      { return W; }

      /// return weight flag
      inline bool Weighted(void) const
      { return weighted; }

      /// add a single sample to the computation of statistics, with optional weight
      void Add(const T& x, const T& wt=T())
      {
         if(wt != T()) weighted=true;

         T xx(x);
         n++;
         if(n == 1) {
            min = max = ave = xx;
            if(weighted) ave *= wt;
            var = T();
            W = T();
         }
         else {
            if(xx < min) min=xx;
            if(xx > max) max=xx;
         }

         if(weighted) { xx *= wt;  W += wt; }
         ave += (xx-ave)/T(n);
         if(n > 1)
            var = (var*T(n-2) + T(n)*(xx-ave)*(xx-ave)/T(n-1))/T(n-1);
      }

      /// add a Vector<T> of samples to the computation of statistics,
      /// with optional weights
      inline void Add(Vector<T>& X, Vector<T> w = Vector<T>())
      {
         if(w.size() > 0 && w.size() < X.size()) {
            Exception e("Inconsistent input: weight vector short");
            GPSTK_THROW(e);
         }

         size_t i;
         if(w.size() > 0) 
            for(i=0; i<X.size(); i++) Add(X(i),w(i));
         else
            for(i=0; i<X.size(); i++) Add(X(i));
      }

      /// remove a sample from the computation of statistics (can't do min and max).
      void Subtract(T x)
      {
         T dn=T(n);
         if(n > 2) var=(var*(dn-T(1))-dn*(x-ave)*(x-ave)/(dn-T(1)))/(dn-T(2));
         else var=T();
        
         if(n > 1) ave=(ave*dn-x)/(dn-T(1));
         else if(n==1) ave=x;
         else ave=T();

         n--;
      }

      /// remove a Vector<T> of samples to the computation of statistics
      inline void Subtract(Vector<T>& X)
      {
         for(size_t i=0; i<X.size(); i++) Subtract(X(i));
      }

      /// define private members directly; useful in continuing with an object
      /// that was earlier saved (e.g. to a file).
      void Load(unsigned int in_n, T in_min, T in_max, T in_ave, T in_var,
                  bool wtd=false, T norm=1.0)
      {
         n = in_n;
         min = in_min;
         max = in_max;
         var = in_var;
         ave = in_ave;
         weighted = wtd;
         W = norm;
      }
      
      /// combine two Stats (assumed taken from the same or equivalent samples);
      /// both must be either weighted or unweighted.
      Stats<T>& operator+=(const Stats<T>& S)
      {
         if(S.n == 0) return *this;
         if(n==0) {
            *this = S;
            return *this;
         }
         if((weighted && !S.weighted) || (!weighted && S.weighted)) {
            Exception e("Stats::operator+= : operands have inconsistent weighting");
            GPSTK_THROW(e);
         }

         if(S.min < min) min=S.min;
         if(S.max > max) max=S.max;

         if(weighted) {
            T W1 = W/n;
            T W2 = S.W/S.n;
            var  = W1*W1*((n-T(1))*var + n*ave*ave)
                 + W2*W2*((S.n-T(1))*S.var + S.n*S.ave*S.ave);
            var *= (n+S.n)/(W+S.W);
            var *= (n+S.n)/(W+S.W);
            W += S.W;
         }
         else {
            var  = ((n-T(1))*var + n*ave*ave);
            var += ((S.n-T(1))*S.var + S.n*S.ave*S.ave);
         }
         ave = (n*ave + S.n*S.ave)/(n+S.n);
         n += S.n;
         var -= n*ave*ave;
         var /= (n-T(1));

         return *this;

      }  // end Stats operator+=

      /// dump the data stored in the class

   private:
      /// Number of samples added to the statistics so far
      unsigned int n;
      /// Minimum value
      T min;
      /// Maximum value
      T max;
      /// Average value
      T var;
      /// Variance (square of the standard deviation)
      T ave;
      /// Normalization constant = sum weights
      T W;
      /// Flag weighted input; ALL input must be consistently weighted or not
      bool weighted;
   }; // end class Stats

   /// Output operator for Stats class
   template <class T>
   std::ostream& operator<<(std::ostream& s, const Stats<T>& ST) 
   {
      std::ofstream savefmt;
      savefmt.copyfmt(s);
      s << " N       = " << ST.N() << (ST.Weighted() ? " ":" not") << " weighted\n";
      s << " Minimum = "; s.copyfmt(savefmt); s << ST.Minimum();