EdgeworthSeries1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_EDGEWORTHSERIES1D_HH_
#define NPSTAT_EDGEWORTHSERIES1D_HH_
 
/*!
// \file EdgeworthSeries1D.hh
//
// \brief Distribution defined by Edgeworth series w.r.t. normal
//
// Author: I. Volobouev
//
// June 2019
*/
 
#include <vector>
 
#include "npstat/nm/SimpleFunctors.hh"
 
#include "npstat/stat/AbsDistribution1D.hh"
#include "npstat/stat/EdgeworthSeriesMethod.hh"
 
namespace npstat {
    class EdgeworthSeries1D : public AbsDistribution1D
    {
    public:
        /**
        // Constructor arguments are as follows:
        //
        // cumulants  -- The vector of cumulants. The first element
        //               of the vector (with index 0) is the first
        //               cumulant, the second element is the second, etc.
        //
        // m          -- The method used to construct the distribition.
        //               See header EdgeworthSeriesMethod.hh for further
        //               explanations.
        //
        // order      -- The expansion order in powers of 1/sqrt(n).
        //               O(1/sqrt(n)) is the 0th order, O(1/n) is the first
        //               order, O(1/(n sqrt(n)) is the second order, etc.
        //               Currently, the highest supported order is 14 for
        //               EDGEWORTH_CLASSICAL method in a non-SLR mode, and
        //               4 in all other situations.
        //
        // slrMode    -- If true, we will assume that the kth cumulant
        //               vanishes (for k > 2) to O(n^{-k/2}) instead of
        //               the usual O(n^{-(k-2)/2}). This is the case for
        //               the SLR statistic.
        */
        EdgeworthSeries1D(const std::vector<double>& cumulants,
                          EdgeworthSeriesMethod m, unsigned order,
                          bool slrMode = false);
 
        inline virtual EdgeworthSeries1D* clone() const
            {return new EdgeworthSeries1D(*this);}
 
        inline virtual ~EdgeworthSeries1D() {}
 
        virtual double density(double x) const;
        virtual double cdf(double x) const;
        virtual double exceedance(double x) const;
        virtual double quantile(double x) const;
 
        /**
        // Inverse exceedance function. More precise than "quantile(1 - x)"
        // for small x.
        */
        double inverseExceedance(double x) const;
 
        //@{
        /** Inspector */
        inline EdgeworthSeriesMethod method() const {return m_;}
        inline bool slrMode() const {return slrMode_;}
        inline unsigned order() const {return order_;}
        //@}
 
        /**
        // This function returns the constructor argument, not the
        // real cumulant which depends on the expansion order used.
        // Use the "empiricalCentralMoment" function and then convert
        // central moments to cumulants if you need real cumulants.
        // The real mean is 0 if order == 0 and cum(0) if order > 0.
        */
        inline double cum(const unsigned k) const {return cumulants_.at(k);};
 
        /** Mean of the function used to build the expansion */
        double edgMean() const;
 
        /** Standard deviation of the function used to build the expansion */
        double edgStdev() const;
 
        /** Edgeworth expansion factor for the density */
        double densityFactor(double x) const;
 
        /** Edgeworth expansion factor for the CDF */
        double cdfFactor(double x) const;
 
        /** 
        // Actual central moment of the density (evaluated numerically
        // for degree > 1 using Gauss-Hermite quadrature)
        */
        double empiricalCentralMoment(unsigned degree) const;
 
        //@{
        /** Prototype needed for I/O */
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream&) const;
        //@}
 
        static inline const char* classname()
            {return "npstat::EdgeworthSeries1D";}
        static inline unsigned version() {return 1;}
        static EdgeworthSeries1D* read(
            const gs::ClassId& id, std::istream& is);
 
    protected:
        virtual bool isEqual(const AbsDistribution1D&) const;
 
    private:
        void hermiteCoeffsClassical(long double* buf, unsigned* maxdeg) const;
        void hermiteCoeffsNormal(long double* buf, unsigned* maxdeg) const;
        void hermiteCoeffsSLR(long double* buf, unsigned* maxdeg) const;
        double normalizedCoord(double x) const;
        double normalizedSigma() const;
        unsigned minNCumulants() const;
 
        std::vector<double> cumulants_;
        EdgeworthSeriesMethod m_;
        unsigned order_;
        bool slrMode_;
    };
 
    /** This class is mostly useful for testing purposes */
    class EdgeworthSeriesMomentFcn : public Functor1<double,double>
    {
    public:
        inline EdgeworthSeriesMomentFcn(const EdgeworthSeries1D& series,
                                        const unsigned degree)
            : series_(series), mean_(0.0), degree_(degree)
            {if (series_.order()) {mean_ = series_.cum(0);}}
 
        inline virtual ~EdgeworthSeriesMomentFcn() {}
 
        virtual double operator()(const double&) const;
 
    private:
        const EdgeworthSeries1D& series_;
        double mean_;
        unsigned degree_;
    };
}
 
#endif // NPSTAT_EDGEWORTHSERIES1D_HH_