GPS ToolKit Software Library: BivarStats.hpp Source File

00001 #pragma ident "$Id: BivarStats.hpp 3140 2012-06-18 15:03:02Z susancummins $"
00002 
00003 
00004 
00010 #ifndef INCLUDE_GPSTK_BIVARSTATS_HPP
00011 #define INCLUDE_GPSTK_BIVARSTATS_HPP
00012 
00013 //============================================================================
00014 //
00015 //  This file is part of GPSTk, the GPS Toolkit.
00016 //
00017 //  The GPSTk is free software; you can redistribute it and/or modify
00018 //  it under the terms of the GNU Lesser General Public License as published
00019 //  by the Free Software Foundation; either version 2.1 of the License, or
00020 //  any later version.
00021 //
00022 //  The GPSTk is distributed in the hope that it will be useful,
00023 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00024 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00025 //  GNU Lesser General Public License for more details.
00026 //
00027 //  You should have received a copy of the GNU Lesser General Public
00028 //  License along with GPSTk; if not, write to the Free Software Foundation,
00029 //  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
00030 //  
00031 //  Copyright 2004, The University of Texas at Austin
00032 //
00033 //============================================================================
00034 
00035 #include "MiscMath.hpp"
00036 #include "Vector.hpp"
00037 #include "Exception.hpp"
00038 #include "Stats.hpp"
00039 
00040 namespace gpstk
00041 {
00042 
00049    template <class T>
00050    class BivarStats
00051    {
00052    public:
00059       BivarStats(bool scale=false);
00060       BivarStats(const T& x, const T&, bool scale=false);
00061       BivarStats(const std::vector<T>& x, const std::vector<T>& y, 
00062                  bool scale=false);
00063       BivarStats(const std::vector< std::pair<T, T> >& d, bool scale=false);
00064       BivarStats(const Vector<T>& x, const Vector<T>& y, bool scale=false);
00070       void add(const T& x, const T& y);
00071       void add(const std::vector<T>& x, const std::vector<T>& y);
00072       void add(const std::vector< std::pair<T, T> >& d);
00073       void add(const Vector<T>& x, const Vector<T>& y);
00079       void subtract(const T& x, const T& y);
00080       void subtract(const std::vector<T>& x, const std::vector<T>& y);
00081       void subtract(const std::vector< std::pair<T, T> >& d);
00082       void subtract(const Vector<T>& x, const Vector<T>& y);
00085       void clear(void); 
00086       size_t n(void) const; 
00087 
00088       T minimumX(void) const;
00089       T maximumX(void) const;
00090       T minimumY(void) const;
00091       T maximumY(void) const;
00092 
00093       T averageX(void) const;
00094       T averageY(void) const;
00095 
00096       T varianceX(void) const;
00097       T varianceY(void) const;
00098       T stdDevX(void) const;
00099       T stdDevY(void) const;
00100 
00102       T slope(void) const;
00104       T intercept(void) const; 
00106       T sigmaSlope(void) const;
00107 
00108       T correlation(void) const;
00109 
00111       T sigmaYX(void) const;
00112 
00114       T eval(const T& x) const {return intercept() + x * slope();};
00115 
00118       BivarStats<T>& operator+=(BivarStats<T>& S);
00119 
00120       Stats<T> estimateDeviation(const std::vector< std::pair<T, T> >& d) const;
00121 
00122    private:
00123 
00125       size_t ns;
00126 
00127       T xMin, xMax, yMin, yMax;
00128       T scaleX, scaleY;
00129       bool scaled;
00130       T sumX, sumY, sumX2, sumY2, sumXY;
00131    }; // end class BivarStats
00132 
00134    template <class T>
00135    std::ostream& operator<<(std::ostream& s, const BivarStats<T>& BVS) 
00136    {
00137       s << " N       = " << BVS.n() << std::endl
00138         << " Minimum: X = " << BVS.minimumX()
00139         << "  Y = " << BVS.minimumY()
00140         << "  Maximum: X = " << BVS.maximumX()
00141         << "  Y = " << BVS.maximumY() << std::endl
00142         << " Average: X = " << BVS.averageX()
00143         << "  Y = " << BVS.averageY()
00144         << "  Std Dev: X = " << BVS.stdDevX()
00145         << "  Y = " << BVS.stdDevY() << std::endl
00146         << " Intercept = " << BVS.intercept()
00147         << "  Slope = " << BVS.slope()
00148         << " with uncertainty = " << BVS.sigmaSlope() << std::endl
00149         << " Conditional uncertainty (sigma y given x) = " << BVS.sigmaYX()
00150         << "  Correlation = "  << BVS.correlation() << std::endl;
00151       return s;
00152    }
00153 
00154    template<class T>
00155    BivarStats<T>::BivarStats(bool s)
00156          :ns(0), scaled(s)
00157    {}
00158 
00159    template<class T>
00160    BivarStats<T>::BivarStats(const T& x, const T& y, bool s)
00161          :ns(0), scaled(s)
00162    {
00163       add(x,y);
00164    }
00165    
00166    template<class T>
00167    BivarStats<T>::BivarStats(const std::vector<T>& x, const std::vector<T>& y, 
00168                              bool s )
00169          :ns(0), scaled(s)
00170    {
00171       add(x,y);
00172    }
00173 
00174    template<class T>
00175    BivarStats<T>::BivarStats(const std::vector< std::pair<T, T> >& d, bool s)
00176          :ns(0), scaled(s)
00177    {
00178       add(d);
00179    }
00180 
00181    template<class T>
00182    BivarStats<T>::BivarStats(const Vector<T>& x, const Vector<T>& y, bool s)
00183       :ns(0), scaled(s)
00184    {
00185       add(x,y);
00186    }
00187 
00188    template<class T>
00189    void BivarStats<T>::add(const T& x, const T& y)
00190    {
00191       if (ns == 0)
00192       {
00193          sumX = sumY = sumX2 = sumY2 = sumXY = T(0);
00194          xMin = xMax = x;
00195          yMin = yMax = y;
00196          scaleX = scaleY = T(1);
00197       }
00198 
00199       if (scaled)
00200       {
00201          if (scaleX==T(1) && x!=T()) scaleX=ABS(x);
00202          if (scaleY==T(1) && y!=T()) scaleY=ABS(y);
00203          T tx(x/scaleX);
00204          T ty(y/scaleY);
00205          sumX += tx;
00206          sumY += ty;
00207          sumX2 += tx*tx;
00208          sumY2 += ty*ty;
00209          sumXY += tx*ty;
00210       }
00211       else
00212       {
00213          sumX += x;
00214          sumY += y;
00215          sumX2 += x*x;
00216          sumY2 += y*y;
00217          sumXY += x*y;
00218       }
00219 
00220       if(x < xMin) xMin=x;
00221       if(x > xMax) xMax=x;
00222       if(y < yMin) yMin=y;
00223       if(y > yMax) yMax=y;
00224       ns++;
00225    }
00226 
00227    template<class T>
00228    void BivarStats<T>::add(const std::vector<T>& x, const std::vector<T>& y)
00229    {
00230       size_t m = x.size() < y.size() ? x.size() : y.size();
00231       if(m==0)
00232          return;
00233       for (size_t i=0; i<m; i++)
00234          add(x[i], y[i]);
00235    }
00236 
00237    template<class T>
00238    void BivarStats<T>::add(const std::vector< std::pair<T, T> >& d)
00239    {
00240       size_t max( d.size() );
00241       for (size_t i=0; i<max; i++)
00242          add(d[i].first, d[i].second);
00243    }
00244 
00245    template<class T>
00246    void BivarStats<T>::add(const Vector<T>& x, const Vector<T>& y)
00247    {
00248       size_t m = x.size() < y.size() ? x.size() : y.size();
00249       if (m==0)
00250          return;
00251       for (size_t i=0; i<m; i++)
00252          add(x(i), y(i));
00253    }
00254 
00255    template<class T>
00256    void BivarStats<T>::subtract(const T& x, const T& y)
00257    {
00258       if (ns < 2)
00259       {
00260          ns = 0;
00261          return;
00262       }
00263 
00264       if (scaled)
00265       {
00266          T tx(x/scaleX);
00267          T ty(y/scaleY);
00268          sumX -= tx;
00269          sumY -= ty;
00270          sumX2 -= tx*tx;
00271          sumY2 -= ty*ty;
00272          sumXY -= tx*ty;
00273       }
00274       else
00275       {
00276          sumX -= x;
00277          sumY -= y;
00278          sumX2 -= x*x;
00279          sumY2 -= y*y;
00280          sumXY -= x*y;
00281       }
00282 
00283       ns--;
00284    }
00285 
00286    template<class T>
00287    void BivarStats<T>::subtract(const std::vector<T>& x, const std::vector<T>& y)
00288    {
00289       size_t m = x.size()<y.size() ? x.size() : y.size();
00290       if(m==0)
00291          return;
00292       for (size_t i=0; i<m; i++)
00293          subtract(x[i], y[i]);
00294    }
00295 
00296    template<class T>
00297    void BivarStats<T>::subtract(const std::vector< std::pair<T, T> >& d)
00298    {
00299       size_t max( d.size() );
00300       for (size_t i=0; i<max; d++)
00301          subtract(d[i].first, d[i].second);
00302    }
00303 
00304    template<class T>
00305    void BivarStats<T>::subtract(const Vector<T>& x, const Vector<T>& y)
00306    {
00307       size_t m = x.size()<y.size() ? x.size() : y.size();
00308       if (m==0)
00309          return;
00310       for (size_t i=0; i<m; i++)
00311          subtract(x(i), y(i));
00312    }
00313 
00315    template<class T>
00316    void BivarStats<T>::clear(void) { ns=0; }
00317 
00318    template<class T>
00319    inline size_t BivarStats<T>::n(void) const { return ns; }
00320 
00321    template<class T>
00322    T BivarStats<T>::minimumX(void) const { return ns>0 ? xMin : T(0); }
00323    template<class T>
00324    T BivarStats<T>::maximumX(void) const { return ns>0 ? xMax : T(0); }
00325    template<class T>
00326    T BivarStats<T>::minimumY(void) const { return ns>0 ? yMin : T(0); }
00327    template<class T>
00328    T BivarStats<T>::maximumY(void) const { return ns>0 ? yMax : T(0); }
00329 
00330    template<class T>
00331    T BivarStats<T>::averageX(void) const
00332    { return ns>0 ? scaleX*sumX/T(ns) : T(0); }
00333    template<class T>
00334    T BivarStats<T>::averageY(void) const 
00335    { return ns>0 ? scaleY*sumY/T(ns) : T(0); }
00336 
00337    
00338    template<class T>
00339    T BivarStats<T>::varianceX(void) const
00340    {  
00341       return (ns>1) ? scaleX*scaleX * (sumX2 - sumX*sumX/T(ns)) / T(ns-1) : T(0); 
00342    }
00343 
00344    template<class T>
00345    T BivarStats<T>::varianceY(void) const
00346    {
00347       return (ns>1) ? scaleY*scaleY * (sumY2 - sumY*sumY/T(ns)) / T(ns-1) : T(0);
00348    }
00349 
00350    template<class T>
00351    T BivarStats<T>::stdDevX(void) const { return SQRT(varianceX()); }
00352    template<class T>
00353    T BivarStats<T>::stdDevY(void) const { return SQRT(varianceY()); }
00354 
00355    template<class T>
00356    T BivarStats<T>::slope(void) const
00357    {
00358       if (ns>0)
00359          return (scaleY/scaleX) * (sumXY - sumX*sumY/T(ns)) / 
00360             (sumX2 - sumX*sumX/T(ns));
00361       else
00362          return T();
00363    }
00364 
00365    template<class T>
00366    T BivarStats<T>::intercept(void) const
00367    {
00368       if (ns>0)
00369          return averageY() - slope() * averageX();
00370       else
00371          return T();
00372    }
00373    
00374    template<class T>
00375    T BivarStats<T>::sigmaSlope(void) const
00376    {
00377       if (ns>2)
00378          return sigmaYX() / (stdDevX() * SQRT(T(ns-1)));
00379       else
00380          return T();
00381    }
00382    
00383    template<class T>
00384    T BivarStats<T>::correlation(void) const
00385    {
00386       if (ns>1)
00387          return scaleX*scaleY * (sumXY - sumX*sumY/T(ns)) /
00388             (stdDevX() * stdDevY() * T(ns-1));
00389       else
00390          return T();
00391    }
00392 
00393    template<class T>
00394    T BivarStats<T>::sigmaYX(void) const
00395    {
00396       if (ns>2)
00397          return (stdDevY() * SQRT(T(ns-1) / T(ns-2))
00398                  * SQRT(T(1) - correlation() * correlation()) );
00399       else return T();
00400    }
00401 
00404    template<class T>
00405    BivarStats<T>& BivarStats<T>::operator+=(BivarStats<T>& S)
00406    {
00407       if(ns + S.ns == 0) return *this;
00408       xMin = std::min(xMin, S.xMin);
00409       xMax = std::max(xMax, S.xMax);
00410       yMin = std::min(yMin, S.yMin);
00411       yMax = std::max(yMax, S.yMax);
00412       T xscaler( S.scaleX/scaleX ), yscaler( S.scaleY/scaleY ); 
00413       sumX += xscaler * S.sumX;
00414       sumY += yscaler * S.sumY;
00415       sumX2 += xscaler * xscaler * S.sumX2;
00416       sumY2 += yscaler * yscaler * S.sumY2;
00417       sumXY += xscaler * yscaler * S.sumXY;
00418       ns += S.ns;
00419       return *this;
00420    } // end of  @addtogroup math
00422 
00423    template<class T>
00424    Stats<T> BivarStats<T>::estimateDeviation(const std::vector< std::pair<T, T> >& d) const
00425    {
00426       Stats<T> estats;
00427       size_t max( d.size() );
00428       for (size_t i=0; i<max; i++)
00429          estats.Add(std::abs(d[i].second - eval(d[i].first)));
00430       return estats;
00431    }
00432 
00433 }  // namespace
00434 
00435 #endif