IT++ Logo Newcom Logo

stat.h

Go to the documentation of this file.
00001 
00033 #ifndef STAT_H
00034 #define STAT_H
00035 
00036 #include <cmath>
00037 #include <itpp/base/vec.h>
00038 #include <itpp/base/mat.h>
00039 #include <itpp/base/elmatfunc.h>
00040 #include <itpp/base/matfunc.h>
00041 
00042 
00043 namespace itpp {
00044 
00047 
00051   class Stat {
00052   public:
00054     Stat() {clear();}
00056     virtual ~Stat() {}
00057 
00059     virtual void clear()
00060     {
00061       _n_overflows = 0;
00062       _n_samples = 0;
00063       _n_zeros = 0;
00064       _max = 0.0;
00065       _min = 0.0;
00066       _sqr_sum = 0.0;
00067       _sum = 0.0;
00068     }
00069 
00071     virtual void sample(const double s, const bool overflow=false)
00072     {
00073       _n_samples++;
00074       _sum += s;
00075       _sqr_sum += s*s;
00076       if (s < _min) _min = s;
00077       if (s > _max) _max = s;
00078       if (overflow) _n_overflows++;
00079       if (s == 0.0) _n_zeros++;
00080     }
00081 
00083     int n_overflows() const {return _n_overflows;}
00085     int n_samples() const {return _n_samples;}
00087     int n_zeros() const {return _n_zeros;}
00089     double avg() const {return _sum/_n_samples;}
00091     double max() const {return _max;}
00093     double min() const {return _min;}
00095     double sigma() const
00096     {
00097       double sigma2 = _sqr_sum/_n_samples - avg()*avg();
00098       return std::sqrt(sigma2 < 0 ? 0 : sigma2);
00099     }
00101     double sqr_sum() const {return _sqr_sum;}
00103     double sum() const {return _sum;}
00105     vec histogram() const {return vec(0);}
00106 
00107   protected:
00109     int _n_overflows;
00111     int _n_samples;
00113     int _n_zeros;
00115     double _max;
00117     double _min;
00119     double _sqr_sum;
00121     double _sum;
00122   };
00123 
00124 
00126   template<class T>
00127     T max(const Vec<T> &v)
00128   {
00129     T   maxdata = v(0);
00130     for (int i=1; i<v.length(); i++)
00131       if (v(i)>maxdata)
00132         maxdata=v(i);
00133     return maxdata;
00134   }
00135 
00137   template<class T>
00138     T max(const Vec<T> &v, int& index)
00139   {
00140     T maxdata = v(0);
00141     index = 0;
00142     for (int i=1; i<v.length(); i++)
00143       if (v(i)>maxdata) {
00144         maxdata=v(i);
00145         index = i;
00146       }
00147     return maxdata;
00148   }
00149 
00157   template<class T>
00158     Vec<T> max(const Mat<T> &m, int dim=1)
00159   {
00160     it_assert(dim==1 || dim==2, "max: dimension need to be 1 or 2");
00161     Vec<T> out;
00162 
00163     if (dim == 1) {
00164       out.set_size(m.cols(), false);
00165 
00166       for (int i=0; i<m.cols(); i++)
00167         out(i) = max(m.get_col(i));
00168     } else {
00169       out.set_size(m.rows(), false);
00170 
00171       for (int i=0; i<m.rows(); i++)
00172         out(i) = max(m.get_row(i));
00173     }
00174       
00175     return out;
00176   }
00177 
00188   template<class T>
00189     Vec<T> max(const Mat<T> &m, ivec &index, int dim=1)
00190   {
00191     it_assert(dim==1 || dim==2, "max: dimension need to be 1 or 2");
00192     Vec<T> out;
00193 
00194     if (dim == 1) {
00195       out.set_size(m.cols(), false);
00196       index.set_size(m.cols(), false);
00197 
00198       for (int i=0; i<m.cols(); i++)
00199         out(i) = max(m.get_col(i), index(i));
00200     } else {
00201       out.set_size(m.rows(), false);
00202       index.set_size(m.rows(), false);
00203 
00204       for (int i=0; i<m.rows(); i++)
00205         out(i) = max(m.get_row(i), index(i));
00206     }
00207       
00208     return out;
00209   }
00210 
00212   template<class T>
00213     T min(const Vec<T> &in)
00214   {
00215     T mindata=in[0];
00216     for (int i=1;i<in.length();i++)
00217       if (in[i]<mindata)
00218         mindata=in[i];
00219     return mindata;
00220   }
00221 
00223   template<class T>
00224     T min(const Vec<T> &in, int& index)
00225   {
00226     T mindata=in[0];
00227     index = 0;
00228     for (int i=1;i<in.length();i++)
00229       if (in[i]<mindata) {
00230         mindata=in[i];
00231         index = i;
00232       }
00233     return mindata;
00234   }
00235 
00236 
00244   template<class T>
00245     Vec<T> min(const Mat<T> &m, int dim=1)
00246   {
00247     it_assert(dim==1 || dim==2, "min: dimension need to be 1 or 2");
00248     Vec<T> out;
00249 
00250     if (dim == 1) {
00251       out.set_size(m.cols(), false);
00252 
00253       for (int i=0; i<m.cols(); i++)
00254         out(i) = min(m.get_col(i));
00255     } else {
00256       out.set_size(m.rows(), false);
00257 
00258       for (int i=0; i<m.rows(); i++)
00259         out(i) = min(m.get_row(i));
00260     }
00261     return out;
00262   }
00263 
00264 
00275   template<class T>
00276     Vec<T> min(const Mat<T> &m,  ivec &index, int dim=1)
00277   {
00278     it_assert(dim==1 || dim==2, "min: dimension need to be 1 or 2");
00279     Vec<T> out;
00280 
00281     if (dim == 1) {
00282       out.set_size(m.cols(), false);
00283       index.set_size(m.cols(), false);
00284 
00285       for (int i=0; i<m.cols(); i++)
00286         out(i) = min(m.get_col(i), index(i));
00287     } else {
00288       out.set_size(m.rows(), false);
00289       index.set_size(m.rows(), false);
00290 
00291       for (int i=0; i<m.rows(); i++)
00292         out(i) = min(m.get_row(i), index(i));
00293     }
00294     return out;
00295   }
00296 
00297 
00299   template<class T>
00300     int max_index(const Vec<T> &in)
00301   {
00302     int maxindex=0;
00303     for (int i=0;i<in.length();i++)
00304       if (in[i]>in[maxindex])
00305         maxindex=i;
00306     return maxindex;
00307   }
00308 
00310   template<class T>
00311     void max_index(const Mat<T> &m, int &row, int &col)
00312   {
00313     T maxdata = m(0,0);
00314     int i, j;
00315 
00316     row = col = 0;
00317     for (i=0; i<m.rows(); i++)
00318       for (j=0; j<m.cols(); j++)
00319         if (m(i,j) > maxdata) {
00320           row = i;
00321           col = j;
00322           maxdata = m(i,j);
00323         }
00324   }
00325 
00327   template<class T>
00328     int min_index(const Vec<T> &in)
00329   {
00330     int minindex=0;
00331     for (int i=0;i<in.length();i++)
00332       if (in[i]<in[minindex])
00333         minindex=i;
00334     return minindex;
00335   }
00336 
00338   template<class T>
00339     void min_index(const Mat<T> &m, int &row, int &col)
00340   {
00341     T mindata = m(0,0);
00342     int i, j;
00343 
00344     row = col = 0;
00345     for (i=0; i<m.rows(); i++)
00346       for (j=0; j<m.cols(); j++)
00347         if (m(i,j) < mindata) {
00348           row = i;
00349           col = j;
00350           mindata = m(i,j);
00351         }
00352   }
00353 
00355    double mean(const vec &v);
00357   std::complex<double> mean(const cvec &v);
00359   double mean(const svec &v);
00361   double mean(const ivec &v);
00363   double mean(const mat &m);
00365   std::complex<double> mean(const cmat &m);
00367   double mean(const smat &m);
00368   // The mean value
00369   double mean(const imat &m);
00370 
00372   template<class T>
00373     double geometric_mean(const Vec<T> &v)
00374   {
00375     //return exp(log(static_cast<double>(product(v)))/v.length());
00376     return std::exp(std::log(double(product(v)))/v.length());
00377   }
00378 
00380   template<class T>
00381     double median(const Vec<T> &v)
00382   {
00383     Vec<T> invect(v);
00384     sort(invect);
00385     return (double)(invect[(invect.length()-1)/2]+invect[invect.length()/2])/2.0;
00386   }
00387 
00389   double norm(const cvec &v);
00390 
00392   template<class T>
00393   double norm(const Vec<T> &v)
00394   {
00395     double E = 0.0;
00396     for (int i = 0; i < v.size(); i++)
00397       E += sqr(static_cast<double>(v[i]));
00398 
00399     return std::sqrt(E);
00400   }
00401 
00403   double norm(const cvec &v, int p);
00404 
00406   template<class T>
00407   double norm(const Vec<T> &v, int p)
00408   {
00409     double E = 0.0;
00410     for (int i = 0; i < v.size(); i++)
00411       E += std::pow(fabs(static_cast<double>(v[i])), static_cast<double>(p));
00412 
00413     return std::pow(E, 1.0 / p);
00414   }
00415 
00417   double norm(const cvec &v, const std::string &s);
00418 
00420   template<class T>
00421   double norm(const Vec<T> &v, const std::string &s)
00422   {
00423     it_assert(s == "fro", "norm(): Unrecognised norm");
00424 
00425     double E = 0.0;
00426     for (int i = 0; i < v.size(); i++)
00427       E += sqr(static_cast<double>(v[i]));
00428 
00429     return std::sqrt(E);
00430   }
00431 
00440   double norm(const mat &m, int p = 2);
00441 
00450   double norm(const cmat &m, int p = 2);
00451 
00453   double norm(const mat &m, const std::string &s);
00454   
00456   double norm(const cmat &m, const std::string &s);
00457 
00458 
00460   double variance(const cvec &v);
00461 
00463   template<class T>
00464     double variance(const Vec<T> &v)
00465   {
00466     int len = v.size();
00467     const T *p=v._data();
00468     double sum=0.0, sq_sum=0.0;
00469 
00470     for (int i=0; i<len; i++, p++) {
00471       sum += *p;
00472       sq_sum += *p * *p;
00473     }
00474 
00475     return (double)(sq_sum - sum*sum/len) / (len-1);
00476   }
00477 
00479   template<class T>
00480     double energy(const Vec<T> &v)
00481   {
00482     return sqr(norm(v));
00483   }
00484 
00485 
00487   inline bool within_tolerance(double x, double xref, double tol = 1e-14)
00488   {
00489     return ( fabs(x-xref) <= tol ) ? true : false; 
00490   }
00491 
00493   inline bool within_tolerance(std::complex<double> x, std::complex<double> xref, double tol = 1e-14)
00494   {
00495     return ( abs(x-xref) <= tol ) ? true : false; 
00496   }
00497 
00499   inline bool within_tolerance(const vec &x, const vec &xref, double tol = 1e-14)
00500   {
00501     return ( max(abs(x-xref)) <= tol ) ? true : false; 
00502   }
00503 
00505   inline bool within_tolerance(const cvec &x, const cvec &xref, double tol = 1e-14)
00506   {
00507     return ( max(abs(x-xref)) <= tol ) ? true : false; 
00508   }
00509 
00511   inline bool within_tolerance(const mat &X, const mat &Xref, double tol = 1e-14)
00512   {
00513     return ( max(max(abs(X-Xref))) <= tol ) ? true : false; 
00514   }
00515 
00517   inline bool within_tolerance(const cmat &X, const cmat &Xref, double tol = 1e-14)
00518   {
00519     return ( max(max(abs(X-Xref))) <= tol ) ? true : false; 
00520   }
00521 
00533   double moment(const vec &x, const int r);
00534 
00563   double skewness(const vec &x);
00564 
00565 
00591   double kurtosisexcess(const vec &x);
00592 
00606   inline double kurtosis(const vec &x) {return kurtosisexcess(x)+3;}
00607 
00609 
00610 } // namespace itpp
00611 
00612 #endif // #ifndef STAT_H
SourceForge Logo

Generated on Wed Apr 18 11:45:34 2007 for IT++ by Doxygen 1.5.2