OrthoPolyGOFTest1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_ORTHOPOLYGOFTEST1D_HH_
#define NPSTAT_ORTHOPOLYGOFTEST1D_HH_
 
/*!
// \file OrthoPolyGOFTest1D.hh
//
// \brief Orthogonal polynomial goodness-of-fit test
//
// Author: I. Volobouev
//
// November 2020
*/
 
#include <memory>
 
#include "npstat/nm/Matrix.hh"
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
#include "npstat/nm/matrixIndexPairs.hh"
 
#include "npstat/stat/AbsUnbinnedGOFTest1D.hh"
 
namespace npstat {
    class OrthoPolyGOFTest1D : public AbsUnbinnedGOFTest1D
    {
    public:
        typedef std::pair<double,UUPair> EigenInfo;
 
        // The mask, if provided, should contain 2*maxdeg elements.
        // Put some number in the mask in order to discard the
        // corresponding eigenvector and put 0 there to keep it.
        template <class Quadrature>
        OrthoPolyGOFTest1D(const AbsDistribution1D& distro,
                           const Quadrature& q,
                           const AbsClassicalOrthoPoly1D& poly,
                           unsigned maxdeg,
                           const unsigned char* mask=0, unsigned lenMask=0);
 
        inline virtual ~OrthoPolyGOFTest1D() {}
 
        virtual std::string shortName() const;
 
        inline unsigned maxDegree() const {return maxdeg_;}
        inline unsigned numDeviations() const {return nUnmasked_;}
        inline bool isEigenvectorUsed(const unsigned i) const
            {return i >= 2U*maxdeg_ ? false : !mask_[i];}
 
        /**
        // Get the uncertainty (i.e., the square root of the covariance
        // matrix eigenvalue) associated with a particular eigenvector.
        // This function ignores the mask. The return value is 0 if the
        // eigenvector number is out of range.
        */
        double eigenUncertainty(unsigned i) const;
 
        /**
        // Get the covariance matrix eigenvector with the given number. If
        // the vector number is out of range, an exception will be thrown.
        // Interpretation of these eigenvectors is for experts only. Make
        // sure that you know what you are doing.
        */
        std::vector<double> covmatEigenvector(unsigned i) const;
 
        /**
        // Collect the information about the given eigenvector.
        // The first element in the "EigenInfo" pair is the
        // eigenvector value and the second element is the
        // pair of polynomials corresponding to that value.
        // Components are included in the result only if the
        // magnitude of the eigenvector value is above the threshold.
        */
        std::vector<EigenInfo> eigenComponents(
            unsigned i, double threshold) const;
 
        /**
        // Get the value of the test statistic with all unmasked
        // terms weighted equally
        */
        inline virtual double testStatistic(
            const double* data, const unsigned long sz, bool) const
            {return testStat(data, sz);}
 
        inline virtual double testStatistic(
            const float* data, const unsigned long sz, bool) const
            {return testStat(data, sz);}
 
        inline virtual bool hasAnalyticPValue() const {return true;}
 
        virtual double analyticPValue(double stat, unsigned long sz) const;
 
        virtual double inverseExceedance(double pvalue, unsigned long sz,
                                         double smin, double smax) const;
 
        /**
        // Get individual normalized deviations from which
        // the test statistic is constructed. Only unmasked
        // deviations are returned. The return value of this
        // function is the number of such deviations.
        */
        template<typename Numeric>
        unsigned normalizedDeviations(
            const Numeric* data, unsigned long lenData,
            double* deviations, unsigned lenDeviations) const;
 
    private:
        template <class Quadrature>
        void calculateNormalizingMatrix(const Quadrature& q,
                                        Matrix<lapack_double>* m) const;
        template<typename Numeric>
        double testStat(const Numeric* data, unsigned long lenData) const;
 
        std::shared_ptr<const AbsClassicalOrthoPoly1D> poly_;
        Matrix<lapack_double> invsqr_;
        std::vector<double> buf_;
        std::vector<double> statBuf_;
        std::vector<double> sqrBuf_;
        std::vector<unsigned char> mask_;
        unsigned maxdeg_;
        unsigned nUnmasked_;
    };
}
 
#include "npstat/stat/OrthoPolyGOFTest1D.icc"
 
#endif // NPSTAT_ORTHOPOLYGOFTEST1D_HH_