misc_stat.h

Go to the documentation of this file.

#ifndef MISC_STAT_H
#define MISC_STAT_H

#include <itpp/base/math/min_max.h>
#include <itpp/base/mat.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/base/matfunc.h>


namespace itpp
{


class Stat
{
public:
  Stat() {clear();}
  virtual ~Stat() {}

  virtual void clear() {
    _n_overflows = 0;
    _n_samples = 0;
    _n_zeros = 0;
    _max = 0.0;
    _min = 0.0;
    _sqr_sum = 0.0;
    _sum = 0.0;
  }

  virtual void sample(const double s, const bool overflow = false) {
    _n_samples++;
    _sum += s;
    _sqr_sum += s * s;
    if (s < _min) _min = s;
    if (s > _max) _max = s;
    if (overflow) _n_overflows++;
    if (s == 0.0) _n_zeros++;
  }

  int n_overflows() const {return _n_overflows;}
  int n_samples() const {return _n_samples;}
  int n_zeros() const {return _n_zeros;}
  double avg() const {return _sum / _n_samples;}
  double max() const {return _max;}
  double min() const {return _min;}
  double sigma() const {
    double sigma2 = _sqr_sum / _n_samples - avg() * avg();
    return std::sqrt(sigma2 < 0 ? 0 : sigma2);
  }
  double sqr_sum() const {return _sqr_sum;}
  double sum() const {return _sum;}
  vec histogram() const {return vec(0);}

protected:
  int _n_overflows;
  int _n_samples;
  int _n_zeros;
  double _max;
  double _min;
  double _sqr_sum;
  double _sum;
};


double mean(const vec &v);
std::complex<double> mean(const cvec &v);
double mean(const svec &v);
double mean(const ivec &v);
double mean(const mat &m);
std::complex<double> mean(const cmat &m);
double mean(const smat &m);
double mean(const imat &m);

template<class T>
double geometric_mean(const Vec<T> &v)
{
  return std::exp(std::log(static_cast<double>(prod(v))) / v.length());
}

template<class T>
double geometric_mean(const Mat<T> &m)
{
  return std::exp(std::log(static_cast<double>(prod(prod(m))))
                  / (m.rows() * m.cols()));
}

template<class T>
double median(const Vec<T> &v)
{
  Vec<T> invect(v);
  sort(invect);
  return (double)(invect[(invect.length()-1)/2] + invect[invect.length()/2]) / 2.0;
}

double norm(const cvec &v);

template<class T>
double norm(const Vec<T> &v)
{
  double E = 0.0;
  for (int i = 0; i < v.size(); i++)
    E += sqr(static_cast<double>(v[i]));

  return std::sqrt(E);
}

double norm(const cvec &v, int p);

template<class T>
double norm(const Vec<T> &v, int p)
{
  double E = 0.0;
  for (int i = 0; i < v.size(); i++)
    E += std::pow(fabs(static_cast<double>(v[i])), static_cast<double>(p));

  return std::pow(E, 1.0 / p);
}

double norm(const cvec &v, const std::string &s);

template<class T>
double norm(const Vec<T> &v, const std::string &s)
{
  it_assert(s == "fro", "norm(): Unrecognised norm");

  double E = 0.0;
  for (int i = 0; i < v.size(); i++)
    E += sqr(static_cast<double>(v[i]));

  return std::sqrt(E);
}

double norm(const mat &m, int p = 2);

double norm(const cmat &m, int p = 2);

double norm(const mat &m, const std::string &s);

double norm(const cmat &m, const std::string &s);


double variance(const cvec &v);

template<class T>
double variance(const Vec<T> &v)
{
  int len = v.size();
  const T *p = v._data();
  double sum = 0.0, sq_sum = 0.0;

  for (int i = 0; i < len; i++, p++) {
    sum += *p;
    sq_sum += *p * *p;
  }

  return (double)(sq_sum - sum*sum / len) / (len - 1);
}

template<class T>
double energy(const Vec<T> &v)
{
  return sqr(norm(v));
}


inline bool within_tolerance(double x, double xref, double tol = 1e-14)
{
  return (fabs(x -xref) <= tol) ? true : false;
}

inline bool within_tolerance(std::complex<double> x, std::complex<double> xref, double tol = 1e-14)
{
  return (abs(x -xref) <= tol) ? true : false;
}

inline bool within_tolerance(const vec &x, const vec &xref, double tol = 1e-14)
{
  return (max(abs(x -xref)) <= tol) ? true : false;
}

inline bool within_tolerance(const cvec &x, const cvec &xref, double tol = 1e-14)
{
  return (max(abs(x -xref)) <= tol) ? true : false;
}

inline bool within_tolerance(const mat &X, const mat &Xref, double tol = 1e-14)
{
  return (max(max(abs(X -Xref))) <= tol) ? true : false;
}

inline bool within_tolerance(const cmat &X, const cmat &Xref, double tol = 1e-14)
{
  return (max(max(abs(X -Xref))) <= tol) ? true : false;
}

double moment(const vec &x, const int r);

double skewness(const vec &x);


double kurtosisexcess(const vec &x);

inline double kurtosis(const vec &x) {return kurtosisexcess(x) + 3;}


} // namespace itpp

#endif // #ifndef MISC_STAT_H

---