/*---------------------------------------------------------------------------*
* IT++ *
*---------------------------------------------------------------------------*
* Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger, *
* Tobias Ringström, and Johan Bergman. *
* *
* Permission to use, copy, modify, and distribute this software and its *
* documentation under the terms of the GNU General Public License is hereby *
* granted. No representations are made about the suitability of this *
* software for any purpose. It is provided "as is" without expressed or *
* implied warranty. See the GNU General Public License for more details. *
*---------------------------------------------------------------------------*/
/*!
\file
\brief Functions and classes for calculation of statistics
\author Tony Ottosson and Johan Bergman
1.10
2004/09/30 16:17:17
*/
#ifndef __stat_h
#define __stat_h
#include <cmath>
#include "base/vec.h"
#include "base/mat.h"
#include "base/scalfunc.h"
#include "base/matfunc.h"
namespace itpp {
//! \addtogroup fixtypes
//!@{
/*!
\brief A class for sampling a signal and calculating statistics
*/
class Stat {
public:
//! Default constructor
Stat() {clear();}
//! Destructor
virtual ~Stat() {}
//! Clear statistics
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;
}
//! Register a sample and flag for overflow
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++;
}
//! Number of reported overflows
int n_overflows() const {return _n_overflows;}
//! Number of samples
int n_samples() const {return _n_samples;}
//! Number of zero samples
int n_zeros() const {return _n_zeros;}
//! Average over all samples
double avg() const {return _sum/_n_samples;}
//! Maximum sample
double max() const {return _max;}
//! Minimum sample
double min() const {return _min;}
//! Standard deviation of all samples
double sigma() const
{
double sigma2 = _sqr_sum/_n_samples - avg()*avg();
return std::sqrt(sigma2 < 0 ? 0 : sigma2);
}
//! Squared sum of all samples
double sqr_sum() const {return _sqr_sum;}
//! Sum of all samples
double sum() const {return _sum;}
//! Histogram over all samples (not implemented yet)
vec histogram() const {return vec(0);}
protected:
//! Number of reported overflows
int _n_overflows;
//! Number of samples
int _n_samples;
//! Number of zero samples
int _n_zeros;
//! Maximum sample
double _max;
//! Minimum sample
double _min;
//! Squared sum of all samples
double _sqr_sum;
//! Sum of all samples
double _sum;
};
//!@}
/*! \addtogroup Statistics
*/
//!@{
//! Maximum value of vector
template<class T>
T max(const Vec<T> &v)
{
T maxdata = v(0);
for (int i=1; i<v.length(); i++)
if (v(i)>maxdata)
maxdata=v(i);
return maxdata;
}
//! Maximum value of vector, also returns the index position of max value
template<class T>
T max(const Vec<T> &v, int index)
{
T maxdata = v(0);
index = 0;
for (int i=1; i<v.length(); i++)
if (v(i)>maxdata) {
maxdata=v(i);
index = i;
}
return maxdata;
}
/*!
\brief max of elements in the matrix \c m
sum(m)=sum(m,1) returns a vector where the elements are max over each column
sum(m,2) returns a vector where the elements are max over each row
*/
template<class T>
Vec<T> max(const Mat<T> &m, int dim=1)
{
it_assert(dim==1 || dim==2, "max: dimension need to be 1 or 2");
Vec<T> out;
if (dim == 1) {
out.set_size(m.cols(), false);
for (int i=0; i<m.cols(); i++)
out(i) = max(m.get_col(i));
} else {
out.set_size(m.rows(), false);
for (int i=0; i<m.rows(); i++)
out(i) = max(m.get_row(i));
}
return out;
}
/*!
\brief max of elements in the matrix \c m
sum(m)=sum(m,1) returns a vector where the elements are max over each column
sum(m,2) returns a vector where the elements are max over each row
Also returns an vector of indices with positions of max value within column/row
*/
template<class T>
Vec<T> max(const Mat<T> &m, ivec &index, int dim=1)
{
it_assert(dim==1 || dim==2, "max: dimension need to be 1 or 2");
Vec<T> out;
if (dim == 1) {
out.set_size(m.cols(), false);
index.set_size(m.cols(), false);
for (int i=0; i<m.cols(); i++)
out(i) = max(m.get_col(i), index(i));
} else {
out.set_size(m.rows(), false);
index.set_size(m.rows(), false);
for (int i=0; i<m.rows(); i++)
out(i) = max(m.get_row(i), index(i));
}
return out;
}
//! Minimum value of vector
template<class T>
T min(const Vec<T> &in)
{
T mindata=in[0];
for (int i=1;i<in.length();i++)
if (in[i]<mindata)
mindata=in[i];
return mindata;
}
//! Minimum value of vector, also returns the index position of min value
template<class T>
T min(const Vec<T> &in, int index)
{
T mindata=in[0];
index = 0;
for (int i=1;i<in.length();i++)
if (in[i]<mindata) {
mindata=in[i];
index = i;
}
return mindata;
}
/*!
\brief min of elements in the matrix \c m
sum(m)=sum(m,1) returns a vector where the elements are min over each column
sum(m,2) returns a vector where the elements are min over each row
*/
template<class T>
Vec<T> min(const Mat<T> &m, int dim=1)
{
it_assert(dim==1 || dim==2, "min: dimension need to be 1 or 2");
Vec<T> out;
if (dim == 1) {
out.set_size(m.cols(), false);
for (int i=0; i<m.cols(); i++)
out(i) = min(m.get_col(i));
} else {
out.set_size(m.rows(), false);
for (int i=0; i<m.rows(); i++)
out(i) = min(m.get_row(i));
}
return out;
}
/*!
\brief min of elements in the matrix \c m
sum(m)=sum(m,1) returns a vector where the elements are min over each column
sum(m,2) returns a vector where the elements are min over each row
Also returns an vector of indices with positions of max value within column/row
*/
template<class T>
Vec<T> min(const Mat<T> &m, ivec &index, int dim=1)
{
it_assert(dim==1 || dim==2, "min: dimension need to be 1 or 2");
Vec<T> out;
if (dim == 1) {
out.set_size(m.cols(), false);
index.set_size(m.cols(), false);
for (int i=0; i<m.cols(); i++)
out(i) = min(m.get_col(i), index(i));
} else {
out.set_size(m.rows(), false);
index.set_size(m.rows(), false);
for (int i=0; i<m.rows(); i++)
out(i) = min(m.get_row(i), index(i));
}
return out;
}
//! Return the postion of the maximum element in the vector
template<class T>
int max_index(const Vec<T> &in)
{
int maxindex=0;
for (int i=0;i<in.length();i++)
if (in[i]>in[maxindex])
maxindex=i;
return maxindex;
}
//! Return the postion of the maximum element in the matrix
template<class T>
void max_index(const Mat<T> &m, int &row, int &col)
{
T maxdata = m(0,0);
int i, j;
row = col = 0;
for (i=0; i<m.rows(); i++)
for (j=0; j<m.cols(); j++)
if (m(i,j) > maxdata) {
row = i;
col = j;
maxdata = m(i,j);
}
}
//! Return the postion of the minimum element in the vector
template<class T>
int min_index(const Vec<T> &in)
{
int minindex=0;
for (int i=0;i<in.length();i++)
if (in[i]<in[minindex])
minindex=i;
return minindex;
}
//! Return the postion of the minimum element in the matrix
template<class T>
void min_index(const Mat<T> &m, int &row, int &col)
{
T mindata = m(0,0);
int i, j;
row = col = 0;
for (i=0; i<m.rows(); i++)
for (j=0; j<m.cols(); j++)
if (m(i,j) < mindata) {
row = i;
col = j;
mindata = m(i,j);
}
}
//! The mean value
double mean(const vec &v);
//! The mean value
complex<double> mean(const cvec &v);
//! The mean value
double mean(const svec &v);
//! The mean value
double mean(const ivec &v);
//! The mean value
double mean(const mat &m);
//! The mean value
complex<double> mean(const cmat &m);
//! The mean value
double mean(const smat &m);
// The mean value
double mean(const imat &m);
//! The geometric mean value
template<class T>
double geometric_mean(const Vec<T> &v)
{
//return exp(log(static_cast<double>(product(v)))/v.length());
return std::exp(std::log(double(product(v)))/v.length());
}
//! The median
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;
}
//! Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
double norm(const cvec &v);
//! Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
template<class T>
double norm(const Vec<T> &v)
{
int i;
double E=0.0;
for (i=0; i<v.size(); i++)
E += sqr(double(v[i]));
return std::sqrt(E);
}
//! Calculate the p-norm: norm(v,p)=sum(abs(v).^2)^(1/p)
double norm(const cvec &v, int p);
//! Calculate the p-norm: norm(v,p)=sum(abs(v).^2)^(1/p)
template<class T>
double norm(const Vec<T> &v, int p)
{
int i;
double E=0;
for (i=0;i<v.size();i++)
E += std::pow(fabs((double)v[i]),(double)p);
return std::pow(E,1.0/(double)p);
}
/*! Calculate the p-norm of a real matrix
p=1: max(svd(m))
p=2: max(sum(abs(X)))
Default if no p is given is the 2-norm
*/
double norm(const mat &m, int p=2);
/*! Calculate the p-norm of a complex matrix
p=1: max(svd(m))
p=2: max(sum(abs(X)))
Default if no p is given is the 2-norm
*/
double norm(const cmat &m, int p=2);
//! The variance of the elements in the vector. Normalized with N-1 to be unbiased.
double variance(const cvec &v);
//! The variance of the elements in the vector. Normalized with N-1 to be unbiased.
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);
}
//! Calculate the energy: squared 2-norm. energy(v)=sum(abs(v).^2)
template<class T>
double energy(const Vec<T> &v)
{
return sqr(norm(v));
}
/*!
\brief Calculate the central moment of vector x
The \f$r\f$th sample central moment of the samples in the vector
\f$ \mathbf{x} \f$ is defined as
\f[
m_r = \mathrm{E}[x-\mu]^r = \frac{1}{n} \sum_{i=0}^{n-1} (x_i - \mu)^r
\f]
where \f$\mu\f$ is the sample mean.
*/
double moment(const vec &x, const int r);
/*!
\brief Calculate the skewness excess of the input vector x
The skewness is a measure of the degree of asymmetry of distribution. Negative
skewness means that the distribution is spread more to the left of the mean than to
the right, and vice versa if the skewness is positive.
The skewness of the samples in the vector \f$ \mathbf{x} \f$ is
\f[
\gamma_1 = \frac{\mathrm{E}[x-\mu]^3}{\sigma^3}
\f]
where \f$\mu\f$ is the mean and \f$\sigma\f$ the standard deviation.
The skewness is estimated as
\f[
\gamma_1 = \frac{k_3}{{k_2}^{3/2}}
\f]
where
\f[
k_2 = \frac{n}{n-1} m_2
\f]
and
\f[
k_3 = \frac{n^2}{(n-1)(n-2)} m_3
\f]
Here \f$m_2\f$ is the sample variance and \f$m_3\f$ is the 3rd sample
central moment.
*/
double skewness(const vec &x);
// This is not the kurtosis defined by matlab (don't subtract 3)
/*!
\brief Calculate the kurtosis excess of the input vector x
The kurtosis is a measure of peakedness of a distribution. Several definitions
of kurtosis exist, but here it is assumed the so called kurtosis excess defined
as
\f[
\gamma_2 = \frac{\mathrm{E}[x-\mu]^4}{\sigma^4} - 3
\f]
where \f$\mu\f$ is the mean and \f$\sigma\f$ the standard deviation.
The kurtosis excess is estimated as
\f[
\gamma_2 = \frac{k_4}{{k_2}^2}
\f]
where
\f[
k_2 = \frac{n}{n-1} m_2
\f]
and
\f[
k_4 = \frac{n^2 [(n+1)m_4 - 3(n-1){m_2}^2]}{(n-1)(n-2)(n-3)}
\f]
Here \f$m_2\f$ is the sample variance and \f$m_4\f$ is the 4th sample
central moment.
This is not the kurtosis defined by matlab which do not subtract 3!
*/
double kurtosis(const vec &x);
} //namespace itpp
#endif // __stat_h
syntax highlighted by Code2HTML, v. 0.9.1