NeymanOSDE1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_NEYMANOSDE1D_HH_
#define NPSTAT_NEYMANOSDE1D_HH_
 
/*!
// \file NeymanOSDE1D.hh
//
// \brief OSDE based on minimizing the expected ISE of the empirical
//        comparison density
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <vector>
 
#include "npstat/nm/Matrix.hh"
#include "npstat/stat/LegendreDistro1D.hh"
 
namespace npstat {
    class NeymanOSDE1D
    {
    public:
        template<typename Numeric>
        NeymanOSDE1D(const Numeric* coords, unsigned long lenCoords,
                     double xmin, double xmax);
 
        inline double xmin() const {return xmin_;}
        inline double xmax() const {return xmax_;}
        inline unsigned long sampleSize() const {return nCoords_;}
        inline double coord(unsigned long i) const
            {return scale_*coords_.at(i) + shift_;}
 
        // coeffs      -- array of coefficients, starting with degree 1.
        // nCoeffs     -- length of the array "coeffs".
        //
        // shrinkages  -- shrinkage coefficients for constructing chi-square.
        // nShrinkages -- length of the array "shrinkages". Should be at
        //                least as large as "nCoeffs".
        //
        // The function returns the chi-square value.
        // Also, the following quantities are filled out on exit:
        //
        // polyStatistics -- scaled expansion coefficients of the
        //                   comparison density, all these would have
        //                   mean 0 and standard deviation 1 for oracle
        //                   values of "coeffs". The length of this
        //                   array should be at least "nShrinkages".
        //
        // gradient    -- chi-square gradient w.r.t. coeffs. The length
        //                of this array should be at least "nCoeffs".
        //
        // hessian     -- chi-square Hessian w.r.t. coeffs (if requested).
        //
        double chiSquare(const double* coeffs, unsigned nCoeffs,
                         const double* shrinkages, unsigned nShrinkages,
                         double* polyStatistics,
                         double* gradient,
                         Matrix<double>* hessian = nullptr) const;
 
        // "Sensitivity" matrix. The returned matrix will be
        // dimensioned nShrinkages x nCoeffs.
        Matrix<double> sensitivityMatrix(
            const double* coeffs, unsigned nCoeffs,
            unsigned nShrinkages) const;
 
        // Get the expectation values of the Legendre polynomials
        // for the sample, starting with degree 1. The length of
        // the array "coeffs" should be maxdeg.
        void sampleCoeffs(double* coeffs, unsigned maxdeg) const;
 
        // Once the coefficients have been determined, get the
        // density estimate for the original sample. Note that
        // this density estimate will not be guaranteed positive.
        // The coefficients start with degree 1.
        LegendreDistro1D densityEstimate(const double* coeffs,
                                         unsigned nCoeffs) const;
 
        // Check whether the given set of coefficients creates
        // a series which can be used as a bona fide density.
        // The coefficients start with degree 1.
        inline bool isPositive(const double* coeffs,
                               const unsigned nCoeffs) const
            {return densityEstimate(coeffs, nCoeffs).isPositive();}
 
    private:
        void moveCoordsToStandardInterval();
        void fillIntegralTable(unsigned maxdeg) const;
        void fillDerivativeTables(unsigned maxdeg) const;
 
        // Coefficients start with degree 1 (not 0)
        double cdf(const double* coeffs, unsigned ncoeffs,
                   unsigned long ipt) const;
 
        // Coordinates will be shifted and scaled to the [-1, 1] interval
        std::vector<double> coords_;
        unsigned long nCoords_;
        double xmin_;
        double xmax_;
        double scale_;
        double shift_;
 
        // integs_ is a table of integrals int_{-1}^coords_[i] P_m(x) dx,
        // with m starting from 1. For the given m and i, use
        // integs_[(m-1)*nCoords_ + i]. For m = 0, the corresponding
        // integral is just coords_[i] + 1.0. The table should be grown
        // as needed by the "fillIntegralTable" method.
        mutable std::vector<double> integs_;
        mutable unsigned currentMaxDeg_;
 
        // Series coefficients for the derivatives of the
        // Legendre polynomials on [0, 1]
        mutable std::vector<std::vector<double> > firstDeriv_;
 
        // Series coefficients for the second derivatives
        // of the Legendre polynomials on [0, 1]
        mutable std::vector<std::vector<double> > secondDeriv_;
    };
}
 
#include "npstat/stat/NeymanOSDE1D.icc"
 
#endif // NPSTAT_NEYMANOSDE1D_HH_