statistic.h

Go to the documentation of this file.

 
#include "elementa/license.inc"
#include "elementa/checks.inc"
 
#ifndef ELEMENTA_MATH_STATS_STATISTIC_H
#define ELEMENTA_MATH_STATS_STATISTIC_H
 
#include <limits>
#include <type_traits>
#include <string>
#include <cmath>
#include <vector>
#include <functional>
#include "elementa/base/exceptions.h"
 
namespace elementa
{
 
namespace math
{
 
namespace stats
{
 
/* *************************************************************************/
/* ***************************************************************************
*
*   Types: SampleSize, StaDatum, GetStaDatumFun
*   Constexpr: getstadatum_direct
*
***************************************************************************** */
 
using SampleSize = size_t;
 
using StaDatum = double;
 
template <class DatumIterator>
using GetStaDatumFun = std::function< StaDatum(const DatumIterator &) >;
 
template <class DatumIterator>
constexpr StaDatum getstadatum_direct(const DatumIterator & ti)
    { return(static_cast<StaDatum>(*ti)); }
 
 
/* *****************************************************************************
*
*   Base policy class: StatPol_BaseOnInc
*   Policy classes: StaPol_Min, StaPol_Max, StaPol_Mean, StaPol_Var
*
*******************************************************************************/
 
// forward declaration
template <class Policy>
class Statistic;
 
 
 
class StaPol_BaseOnInc
{
 public:
 
    template <class TIt> 
    StaDatum calc(const TIt & first, const TIt & pastend, SampleSize & n,
                  const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>);
 
    virtual StaDatum incCalc(StaDatum prev, StaDatum v, SampleSize n) = 0;
};
 
 
 
class StaPol_Min: public StaPol_BaseOnInc
{
 public:
 
    StaDatum incCalc(StaDatum prev, StaDatum v, SampleSize n) 
        { return(v < prev ? v : prev); }
};
 
 
 
class StaPol_Max: public StaPol_BaseOnInc
{
 public:
 
    StaDatum incCalc(StaDatum prev, StaDatum v, SampleSize n) 
        { return(v > prev ? v : prev); }
};
 
 
 
class StaPol_Mean
{
 public:
 
    template <class TIt> 
    StaDatum calc(const TIt & first, const TIt & pastend, SampleSize & n,
                  const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>);
 
    StaDatum incCalc(StaDatum prev, StaDatum v, SampleSize n) 
        { // here n is from 0 to N-1, but we need from 1 to N, thus we add 1
          ++n; return(1.0/static_cast<StaDatum>(n) * 
                        (v + static_cast<StaDatum>(n-1) * prev)); }
};
 
 
 
class StaPol_Var
{
 public:
 
 
    StaPol_Var(Statistic<StaPol_Mean> & meanstat):meanstat_{&meanstat} {}
 
 
    template <class TIt> 
    StaDatum calc(const TIt & first, const TIt & pastend, SampleSize & n,
                  const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>);
 
    StaDatum incCalc(StaDatum prev, StaDatum v, SampleSize n);
 
 private:
 
    Statistic<StaPol_Mean> * meanstat_;
    StaDatum prev_;
};
 
 
/* *****************************************************************************
*
*   Template class: Statistic
*
*******************************************************************************/
 
 
template <class Calc >
class Statistic
{
 public:
 
    Statistic(const Calc c):calculator_{c} { reset(); }
 
 
    void reset(void) { calculated_ = false; domsize_ = 0; }
 
    bool calculated(void) const noexcept { return(calculated_); }
 
    StaDatum value(void) const { checkCalculated(); return(value_); }
 
 
    StaDatum previous(void) const;
 
 
    SampleSize domainSize(void) const noexcept { return(domsize_); }
 
    template <class TIt>
    StaDatum calculate(const TIt & first, const TIt & pastend,
                       const GetStaDatumFun<TIt> & getdatum = 
                            getstadatum_direct<TIt>);
 
 
    StaDatum update(StaDatum v);
 
 
 private:
 
    Calc calculator_;
    bool calculated_;
    StaDatum value_,prev_;
    SampleSize domsize_;
 
    void checkCalculated(void) const
        { if (!calculated_) ELE_CODE_INVSTATE("Statistic not calculated yet"); }
};
 
 
/* *****************************************************************************
*
*   Classes: StaGroup_MinMax, StaGroup_MeanVar, StaGroup_MinMaxMeanVar
*
*******************************************************************************/
 
 
class StaGroup_MinMax
{
 public:
 
    Statistic<StaPol_Min> stamin{StaPol_Min{}}; 
    Statistic<StaPol_Max> stamax{StaPol_Max{}}; 
 
    template <class TIt>
    void calculate(const TIt & first, const TIt & pastend,
                   const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>)
        { stamin.calculate(first,pastend,getdatum);
          stamax.calculate(first,pastend,getdatum); }
 
    void update(StaDatum v) { stamin.update(v); stamax.update(v); }
 
    void reset(void) { stamin.reset(); stamax.reset(); }
};
 
 
 
class StaGroup_MeanVar
{
 public:
 
    Statistic<StaPol_Mean> stamean{StaPol_Mean{}}; 
    Statistic<StaPol_Var> stavar{StaPol_Var{stamean}}; 
 
    template <class TIt>
    void calculate(const TIt & first, const TIt & pastend,
                   const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>)
        { stamean.calculate(first,pastend,getdatum);
          stavar.calculate(first,pastend,getdatum); }
 
    void update(StaDatum v) { stamean.update(v); stavar.update(v); }
 
    void reset(void) { stamean.reset(); stavar.reset(); }
};
 
 
 
class StaGroup_MinMaxMeanVar
{
 public:
 
    Statistic<StaPol_Min> stamin{StaPol_Min{}}; 
    Statistic<StaPol_Max> stamax{StaPol_Max{}}; 
    Statistic<StaPol_Mean> stamean{StaPol_Mean{}}; 
    Statistic<StaPol_Var> stavar{StaPol_Var{stamean}}; 
 
    template <class TIt>
    void calculate(const TIt & first, const TIt & pastend,
                   const GetStaDatumFun<TIt> & getdatum = 
                        getstadatum_direct<TIt>)
        { stamin.calculate(first,pastend,getdatum);
          stamax.calculate(first,pastend,getdatum); 
          stamean.calculate(first,pastend,getdatum);
          stavar.calculate(first,pastend,getdatum); }
 
    void update(StaDatum v) 
        { stamin.update(v); stamax.update(v); 
          stamean.update(v); stavar.update(v); }
 
    void reset(void) 
        { stamin.reset(); stamax.reset(); 
          stamean.reset(); stavar.reset(); }
};
 
 
 
/* ===========================================================================
 
        ========== IMPLEMENTATION OF TEMPLATES ==========
 
============================================================================= */
 
 
/* *****************************************************************************
*
*   Template class: Statistic
*
*******************************************************************************/
 
template <class Calc >
StaDatum Statistic<Calc>::previous(void) const
{ 
    checkCalculated();
    if (domsize_ < 2) ELE_CODE_INVSTATE("No previous value");
    return(prev_); 
}
 
template <class Calc>
template <class TIt>
StaDatum Statistic<Calc>::calculate(const TIt & first, const TIt & pastend,
                                  const GetStaDatumFun<TIt> & getdatum)
{ 
    if (first == pastend) ELE_CODE_INVARG("Cannot calculate statistic on empty "
                                          "sample");
    value_ = calculator_.calc(first,pastend,domsize_,getdatum); 
    prev_ = value_;
    calculated_ = true; 
    return(value_); 
}
 
template <class Calc>
StaDatum Statistic<Calc>::update(StaDatum v)
{ 
    if (domsize_ == std::numeric_limits<SampleSize>::max())
        ELE_CODE_INVSTATE("Maximum size of sample reached");
    if (domsize_ == 0) 
    { 
        value_ = v; 
        prev_ = value_; 
    }
    else 
    { 
        prev_ = value_; 
        value_= calculator_.incCalc(prev_,v,domsize_);
    } 
    ++domsize_; 
    calculated_ = true; 
    return(value_); 
}
 
 
/* *****************************************************************************
*
*   Base policy class: StatPol_BaseOnInc
*
*******************************************************************************/
 
template <class TIt> 
StaDatum StaPol_BaseOnInc::calc(const TIt & first, const TIt & pastend, 
                              SampleSize & n, 
                              const GetStaDatumFun<TIt> & getdatum)
{
    StaDatum v = getdatum(first); 
    auto it{first};
    ++it; 
    n = 1;
    while (it != pastend) 
    { 
        v = incCalc(v,getdatum(it),n); 
        ++n; 
        ++it; 
    }
    return(v); 
}
 
 
/* *****************************************************************************
*
*   Policy class: StaPol_Mean
*
*******************************************************************************/
 
template <class TIt> 
StaDatum StaPol_Mean::calc(const TIt & first, const TIt & pastend, 
                           SampleSize & n, const GetStaDatumFun<TIt> & getdatum)
{ 
    StaDatum s{getdatum(first)}; 
    n = 1; 
    auto it{first};
    ++it; 
    while (it != pastend) 
    { 
        s += getdatum(it); 
        ++it; 
        ++n; 
    }
    return(s / static_cast<StaDatum>(n)); 
}
 
 
/* *****************************************************************************
*
*   Policy class: StaPol_Var
*
*******************************************************************************/
 
template <class TIt> 
StaDatum StaPol_Var::calc(const TIt & first, const TIt & pastend, 
                          SampleSize & n, const GetStaDatumFun<TIt> & getdatum)
{ 
    n = meanstat_->domainSize();
    if (n < 2) return(0); // no variance
    StaDatum mu = meanstat_->value();
    auto it0{first}; 
    StaDatum s{getdatum(it0) - mu},aux; 
    s *= s; 
    ++it0;
    while (it0 != pastend) 
    { 
        aux = getdatum(it0) - mu; 
        s += aux*aux; 
        ++it0; 
    }
    prev_ = s / static_cast<StaDatum>(n-1); return(prev_); 
}
 
 
 
} // end stats namespace
 
} // end math namespace
 
} // end elementa namespace
 
 
#endif