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
Generated on Thu Apr 19 14:43:44 2007 for IT++ by Doxygen 1.5.1