LOrPE1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1D_HH_
#define NPSTAT_LORPE1D_HH_
 
/*!
// \file LOrPE1D.hh
//
// \brief Convenience class which aggregates the kernel and the data
//        for unbinned density estimation by LOrPE
//
// Author: I. Volobouev
//
// May 2022
*/
 
#include <vector>
#include <cassert>
#include <stdexcept>
#include <algorithm>
 
#include "npstat/stat/LOrPE1DSymbetaKernel.hh"
 
namespace npstat {
    template <class Numeric>
    class LOrPE1D
    {
    public:
        typedef Numeric point_type;
 
        inline LOrPE1D(const LOrPE1DSymbetaKernel& i_kernel,
                       const Numeric* i_coords, const unsigned long i_nCoords)
            : kernel_(i_kernel), coords_(i_coords, i_coords + i_nCoords)
        {
            if (coords_.empty()) throw std::invalid_argument(
                "In npstat::LOrPE1D constructor: empty point sample");
            assert(i_coords);
            sortAndNormalize();
        }
 
        inline LOrPE1D(const LOrPE1DSymbetaKernel& i_kernel,
                       const std::vector<Numeric>& i_coords)
            : kernel_(i_kernel), coords_(i_coords)
        {
            if (coords_.empty()) throw std::invalid_argument(
                "In npstat::LOrPE1D constructor: empty point sample");
            sortAndNormalize();
        }
 
        inline ~LOrPE1D() {}
 
        inline void setNormFactor(const double d)
        {
            kernel_.setNormFactor(d);
        }
 
        inline void setSample(const Numeric* i_coords,
                              const unsigned long i_nCoords)
        {
            if (!i_nCoords) throw std::invalid_argument(
                "In npstat::LOrPE1D::setSample: empty point sample");
            assert(i_coords);
            coords_.clear();
            coords_.reserve(i_nCoords);
            std::copy(i_coords, i_coords+i_nCoords, std::back_inserter(coords_));
            sortAndNormalize();
        }
 
        inline void setSample(const std::vector<Numeric>& i_coords)
        {
            if (i_coords.empty()) throw std::invalid_argument(
                "In npstat::LOrPE1D::setSample: empty point sample");
            coords_ = i_coords;
            sortAndNormalize();
        }
 
        /**
        // This method with throw std::runtime_error if Numeric
        // is not an std::pair. Note that this method resets the
        // normalization factor.
        */
        inline void setWeights(const double* weights, unsigned long nWeights)
        {
            if (nWeights != coords_.size()) throw std::invalid_argument(
                "In npstat::LOrPE1D::setWeights: "
                "incompatible length of the weight array");
            Private::SetPairSeconds<Numeric>::set(&coords_[0], weights, nWeights);
            setStandardNormfactor();
        }
 
        /** Is the density estimate auto-normalizing? */
        inline bool normalizingLOrPEEstimate() const {return false;}
 
        inline const LOrPE1DSymbetaKernel& kernel() const {return kernel_;}
        inline double normFactor() const {return kernel_.normFactor();}
        inline const std::vector<Numeric>& coords() const {return coords_;}
        inline unsigned long nCoords() const {return coords_.size();}
        inline Numeric minCoordinate() const {return coords_[0];}
        inline Numeric maxCoordinate() const {return coords_.back();}
        inline double totalSampleWeight() const
            {return Private::SampleWeight<Numeric>(coords_).value;}
        inline double effectiveSampleSize() const
        {
            const double sumw = Private::SampleWeight<Numeric>(coords_).value;
            const double sumWsq = Private::SampleSquaredWeight<Numeric>(coords_).value;
            return sumw/sumWsq*sumw;
        }
 
        inline double pointCoordinate(const unsigned long i) const
        {
            double x, w;
            Private::coordAndWeight(coords_[i], &x, &w);
            return x;
        }
 
        inline double pointWeight(const unsigned long i) const
        {
            double x, w;
            Private::coordAndWeight(coords_[i], &x, &w);
            return w;
        }
 
        //@{
        /** Method included for compatibility with LOrPE1DCV class */
        inline double bandwidthCurve(const double /* x */) const {return 1.0;}
        inline double cvLocalizingWeight(const double /* x */) const {return 1.0;}
        //@}
 
        inline double density(const double x, const double bw) const
            {return kernel_.lorpe(x, bw, &coords_[0], coords_.size(), true);}
 
        /** 
        // The lifetime of the lightweight functor returned by this
        // method must not exceed the lifetime of this object.
        //
        // Fcn can either be a functor which calculates the bandwidth
        // as a function of x or just a double (for the constant
        // bandwidth).
        */
        template <typename Fcn>
        inline LOrPE1DFunctorHelper<Numeric,Fcn> densityFunctor(
            Fcn bandwidthFunctor) const
        {
            return LOrPE1DFunctorHelper<Numeric,Fcn>(
                kernel_, bandwidthFunctor, &coords_[0], coords_.size(), true);
        }
 
        template <class Fcn>
        double densityIntegral(
            unsigned nIntegIntervals, unsigned nIntegPoints,
            Fcn bandwidthFunctor) const;
 
        template <class Fcn>
        double integratedSquaredError(
            const AbsDistribution1D& distro,
            const unsigned nIntegIntervals, const unsigned nIntegPoints,
            Fcn bandwidthFunctor) const;
 
        /** 
        // Set the normalization factor so that the density estimate
        // integrates to 1. The function returns the value of the density
        // integral before normalization.
        */
        template <class Fcn>
        double normalizeDensityEstimate(
            unsigned nIntegIntervals, unsigned nIntegPoints,
            Fcn bandwidthFunctor);
 
        inline bool operator==(const LOrPE1D& r) const
            {return kernel_ == r.kernel_ && coords_ == r.coords_;}
        inline bool operator!=(const LOrPE1D& r) const
            {return !(*this == r);}
 
        /** Is this class capable of cross-validation? */
        inline static bool cvCapable() {return false;}
 
    private:
        inline void sortAndNormalize()
        {
            std::sort(coords_.begin(), coords_.end());
            setStandardNormfactor();
        }
 
        inline void setStandardNormfactor()
        {
            const double w = Private::SampleWeight<Numeric>(coords_).value;
            if (w <= 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1D::setStandardNormfactor: "
                "total sample weight must be positive");
            kernel_.setNormFactor(1.0/w);
        }
 
        LOrPE1DSymbetaKernel kernel_;
        std::vector<Numeric> coords_;
    };
}
 
#include "npstat/stat/LOrPE1D.icc"
 
#endif // NPSTAT_LORPE1D_HH_