AbsKDE1DKernel.hh

Go to the documentation of this file.

#ifndef NPSTAT_ABSKDE1DKERNEL_HH_
#define NPSTAT_ABSKDE1DKERNEL_HH_
 
/*!
// \file AbsKDE1DKernel.hh
//
// \brief Brute-force non-discretized KDE in 1-d. No boundary correction.
//
// Author: I. Volobouev
//
// November 2018
*/
 
#include <cmath>
#include <stdexcept>
 
#include "npstat/stat/AbsDistribution1D.hh"
#include "npstat/nm/SimpleFunctors.hh"
 
namespace npstat {
    /** Base class for 1-d kernels and KDE */
    class AbsKDE1DKernel
    {
    public:
        inline virtual ~AbsKDE1DKernel() {}
 
        virtual void setNormFactor(double normfactor) = 0;
 
        virtual double xmin() const = 0;
        virtual double xmax() const = 0;
        virtual double normFactor() const = 0;
 
        virtual double operator()(double x) const = 0;
 
        /** "Virtual copy constructor" */
        virtual AbsKDE1DKernel* clone() const = 0;
 
        //@{
        /**
        // "nIntegPoints" argument in the following method is the
        // number of points to use in the quadrature rule. If this
        // method is not overriden, Gauss-Legendre quadrature is used.
        */
        virtual double momentIntegral(unsigned k, unsigned nIntegPoints) const;
        virtual double squaredIntegral(unsigned nIntegPoints) const;
        //@}
 
        /**
        // Summing K((x - x[i])/h), as in the standard KDE definition.
        // Setting "coordinatesSorted" flag to "true" if coordinates
        // are sorted in the increasing order can potentially increase
        // the speed of KDE evaluation.
        */
        template <typename Numeric>
        double kde(double x, double bandwidth,
                   const Numeric* coords, unsigned long nCoords,
                   bool coordinatesSorted = false) const;
 
        /** Summing K((x[i] - x)/h) */
        template <typename Numeric>
        double reverseKde(double x, double bandwidth,
                          const Numeric* coords, unsigned long nCoords,
                          bool coordinatesSorted = false) const;
 
        /**
        // Leave-one-out density estimate. The original density estimate
        // should be obtained earlier with either "kde" or "reverseKde"
        // method used with the same values of "bandwidth" and "nCoords"
        // arguments.
        */
        double looKde(double bandwidth, unsigned long nCoords,
                      double originalDensityEstimate) const;
 
        /**
        // Integrated squared error w.r.t. a known density. Parameter
        // "nIntegPoints" is the number of integration points per interval,
        // and it must be supported by the GaussLegendreQuadrature class.
        // The integration interval is defined by the support of the density.
        */
        template <typename Numeric>
        double integratedSquaredError(const AbsDistribution1D& distro,
                   unsigned nIntegIntervals, unsigned nIntegPoints,
                   double bandwidth, bool useReverseKde,
                   const Numeric* coords, unsigned long nCoords,
                   bool coordinatesSorted = false) const;
 
        /**
        // Integral of KDE squared. Parameter "nIntegPoints" is the number
        // of integration points per interval, and it must be supported by
        // the GaussLegendreQuadrature class.
        */
        template <typename Numeric>
        double integratedKdeSquared(
                   double xmin, double xmax,
                   unsigned nIntegIntervals, unsigned nIntegPoints,
                   double bandwidth, bool useReverseKde,
                   const Numeric* coords, unsigned long nCoords,
                   bool coordinatesSorted = false) const;
    protected:
        class MomentFcn : public Functor1<double, double>
        {
        public:
            inline MomentFcn(const AbsKDE1DKernel& ref, const unsigned power)
                : kernel(ref), xpow(power) {}
 
            inline virtual ~MomentFcn() {}
 
            inline double operator()(const double& x) const
                {return kernel(x)*(xpow ? pow(x, xpow) : 1.0);}
        private:
            const AbsKDE1DKernel& kernel;
            unsigned xpow;
        };
 
        class SquareFcn : public Functor1<double, double>
        {
        public:
            inline explicit SquareFcn(const AbsKDE1DKernel& ref)
                : kernel(ref) {}
 
            inline virtual ~SquareFcn() {}
 
            inline double operator()(const double& x) const
                {const double k = kernel(x); return k*k;}
        private:
            const AbsKDE1DKernel& kernel;
        };
 
#ifdef SWIG
    // Simplified non-template methods for the python API. Not for use in C++.
    public:
        inline double kde2(double x, double bandwidth,
                           const double* coords, unsigned nCoords,
                           bool coordinatesSorted = false) const
        {
            return kde(x, bandwidth, coords, nCoords, coordinatesSorted);
        }
 
        inline double reverseKde2(double x, double bandwidth,
                                  const double* coords, unsigned nCoords,
                                  bool coordinatesSorted = false) const
        {
            return reverseKde(x, bandwidth, coords, nCoords, coordinatesSorted);
        }
 
        inline double integratedSquaredError2(
            const AbsDistribution1D& distro,
            unsigned nIntegIntervals, unsigned nIntegPoints,
            double bandwidth, bool useReverseKde,
            const double* coords, unsigned nCoords,
            bool coordinatesSorted = false) const
        {
            return integratedSquaredError(distro, nIntegIntervals, nIntegPoints,
                                          bandwidth, useReverseKde, coords, nCoords,
                                          coordinatesSorted);
        }
 
        inline double integratedKdeSquared2(
                   double xmin, double xmax,
                   unsigned nIntegIntervals, unsigned nIntegPoints,
                   double bandwidth, bool useReverseKde,
                   const double* coords, unsigned nCoords,
                   bool coordinatesSorted = false) const
        {
            return integratedKdeSquared(xmin, xmax, nIntegIntervals,  nIntegPoints,
                                        bandwidth, useReverseKde, coords, nCoords,
                                        coordinatesSorted);
        }
#endif // SWIG
    };
 
    /** Arbitrary 1-d densities used as kernel smoothers */
    class KDE1DDensityKernel : public AbsKDE1DKernel
    {
    public:
        inline explicit KDE1DDensityKernel(const AbsDistribution1D& d,
                                           const double normfactor=1.0)
            : distro_(0), norm_(0.0)
        {
            setNormFactor(normfactor);
            distro_ = d.clone();
        }
 
        inline KDE1DDensityKernel(const KDE1DDensityKernel& r)
            : distro_(0), norm_(r.norm_) {distro_ = r.distro_->clone();}
 
        inline KDE1DDensityKernel& operator=(const KDE1DDensityKernel& r)
        {
            if (this != &r)
            {
                delete distro_;
                distro_ = 0;
                distro_ = r.distro_->clone();
                norm_ = r.norm_;
            }
            return *this;
        }
 
        inline virtual KDE1DDensityKernel* clone() const
            {return new KDE1DDensityKernel(*this);}
 
        inline virtual ~KDE1DDensityKernel() {delete distro_;}
 
        inline virtual void setNormFactor(const double normfactor)
        {
            if (normfactor <= 0.0) throw std::invalid_argument(
                "In npstat::KDE1DDensityKernel::setNormFactor: "
                "normalization factor must be positive");
            norm_ = normfactor;
        }
 
        inline double xmin() const {return distro_->quantile(0.0);}
        inline double xmax() const {return distro_->quantile(1.0);}
        inline double normFactor() const {return norm_;}
 
        inline double operator()(const double x) const
            {return norm_*distro_->density(x);}
 
    private:
        AbsDistribution1D* distro_;
        double norm_;
    };
 
    /**
    // A lightweight functor for getting KDE values. It will not copy
    // the kernel and will not own the data. It is the responsibility
    // of the user of this class to ensure that the lifetimes of the
    // kernel and the data exceed the lifetime of the functor.
    //
    // Intended usage is via the "KDE1DFunctor" utility function.
    */
    template <typename Numeric>
    class KDE1DFunctorHelper : public Functor1<double, double>
    {
    public:
        inline KDE1DFunctorHelper(const AbsKDE1DKernel& kernel,
                                  const double bandwidth,
                                  const Numeric* coords,
                                  const unsigned long nCoords,
                                  const bool coordinatesSorted,
                                  const bool useReverseKde)
            : kernel_(kernel), bandwidth_(bandwidth),
              coords_(coords), nCoords_(nCoords),
              coordinatesSorted_(coordinatesSorted),
              useReverseKde_(useReverseKde) {}
 
        inline virtual ~KDE1DFunctorHelper() {}
 
        inline double operator()(const double& x) const
        {
            double kde;
            if (useReverseKde_)
                kde = kernel_.reverseKde(x, bandwidth_, coords_, nCoords_,
                                         coordinatesSorted_);
            else
                kde = kernel_.kde(x, bandwidth_, coords_, nCoords_,
                                  coordinatesSorted_);
            return kde;
        }
 
    private:
        const AbsKDE1DKernel& kernel_;
        double bandwidth_;
        const Numeric* coords_;
        unsigned long nCoords_;
        bool coordinatesSorted_;
        bool useReverseKde_;
    };
 
    /** A convenience function for making lightweight KDE functors */
    template <typename Numeric>
    inline KDE1DFunctorHelper<Numeric> KDE1DFunctor(
        const AbsKDE1DKernel& kernel, const double bandwidth,
        const Numeric* coords, const unsigned long nCoords,
        const bool coordinatesSorted = false, const bool useReverseKde = false)
    {
        return KDE1DFunctorHelper<Numeric>(kernel, bandwidth, coords, nCoords,
                                           coordinatesSorted, useReverseKde);
    }
 
    /** A functor of (x,y) for use with kernelSensitivityMatrix, etc */
    class KDE1DFunctor2 : public Functor2<double, double, double>
    {
    public:
        inline KDE1DFunctor2(const AbsKDE1DKernel& kernel,
                             const double bandwidth,
                             const bool useReverseKde = false)
            : kernel_(kernel), bandwidth_(bandwidth),
              useReverseKde_(useReverseKde)
        {
            if (bandwidth_ <= 0.0) throw std::invalid_argument(
                "In npstat::KDE1DFunctor2 constructor: bandwidth must be positive");
        }
 
        inline virtual ~KDE1DFunctor2() {}
 
        inline double operator()(const double& x, const double& y) const
        {
            const double delta = useReverseKde_ ? x - y : y - x;
            return kernel_(delta/bandwidth_)/bandwidth_;
        }
 
    private:
        const AbsKDE1DKernel& kernel_;
        double bandwidth_;
        bool useReverseKde_;
    };
}
 
#include "npstat/stat/AbsKDE1DKernel.icc"
 
#endif // NPSTAT_ABSKDE1DKERNEL_HH_