ModulatedDistribution1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_MODULATEDDISTRIBUTION1D_HH_
#define NPSTAT_MODULATEDDISTRIBUTION1D_HH_
 
/*!
// \file ModulatedDistribution1D.hh
//
// \brief 1-d continuous statistical distributions multiplied
//        by an arbitrary smooth non-negative function
//
// Author: I. Volobouev
//
// June 2022
*/
 
#include <climits>
#include <algorithm>
 
#include "npstat/stat/AbsDistribution1D.hh"
 
#include "npstat/nm/GaussLegendreQuadrature.hh"
#include "npstat/nm/findRootUsingBisections.hh"
 
namespace npstat {
    template<class Functor>
    class ModulatedDistribution1D : public AbsDistribution1D
    {
    public:
        /**
        // Constructor arguments are as follows:
        //
        //  distro             -- distribution whose density we want to
        //                        modulate
        //
        //  fcn                -- functor that will be multiplied by
        //                        the density
        //
        //  fcnXmin, fcnXmax   -- limits for the functor support
        //
        //  nIntegPoints       -- number of points to use for various
        //                        integrations. Should be supported by
        //                        the GaussLegendreQuadrature class.
        //
        //  nIntegIntervals    -- number of intervals to use for various
        //                        integrations. The support will be
        //                        split into this number of intervals
        //                        and then the quadrature with nIntegPoints
        //                        points will be used on each interval.
        */
        inline ModulatedDistribution1D(const AbsDistribution1D& distro,
                                       const Functor& fcn,
                                       const double fcnXmin = -DBL_MAX,
                                       const double fcnXmax = DBL_MAX,
                                       const unsigned nIntegPoints = 1024,
                                       const unsigned nIntegIntervals = 1)
            : distro_(0), fcn_(fcn), norm_(1.0),
              nIntegPoints_(nIntegPoints), nIntegIntervals_(nIntegIntervals)
        {
            if (!(fcnXmin < fcnXmax)) throw std::invalid_argument(
                "In npstat::ModulatedDistribution1D constructor: "
                "invalid support specification for the modulating function");
            xmin_ = std::max(distro.quantile(0.0), fcnXmin);
            xmax_ = std::min(distro.quantile(1.0), fcnXmax);
            if (!(xmin_ < xmax_)) throw std::invalid_argument(
                "In npstat::ModulatedDistribution1D constructor: "
                "supports of the density and of the modulating function "
                "do not overlap");
            GaussLegendreQuadrature glq(nIntegPoints);
            if (!nIntegIntervals_) throw std::invalid_argument(
                "In npstat::ModulatedDistribution1D constructor: "
                "number of integration intervals must be positive");
            distro_ = distro.clone();
            norm_ = glq.integrate(DensityFunctor1D(*this),
                                  xmin_, xmax_, nIntegIntervals_);
            assert(norm_ > 0.0);
        }
 
        inline ModulatedDistribution1D(const ModulatedDistribution1D& r)
            : distro_(r.distro_->clone()), fcn_(r.fcn_),
              xmin_(r.xmin_), xmax_(r.xmax_), norm_(r.norm_),
              nIntegPoints_(r.nIntegPoints_),
              nIntegIntervals_(r.nIntegIntervals_) {}
 
        inline ModulatedDistribution1D& operator=(const ModulatedDistribution1D& r)
        {
            if (&r != this)
            {
                AbsDistribution1D* p = r.distro_->clone();
                delete distro_;
                distro_ = p;
                fcn_ = r.fcn_;
                xmin_ = r.xmin_;
                xmax_ = r.xmax_;
                norm_ = r.norm_;
                nIntegPoints_ = r.nIntegPoints_;
                nIntegIntervals_ = r.nIntegIntervals_;
            }
            return *this;
        }
 
        inline virtual ModulatedDistribution1D* clone() const
            {return new ModulatedDistribution1D(*this);}
 
        inline virtual ~ModulatedDistribution1D() {delete distro_;}
 
        /** Distribution density */
        inline virtual double density(const double x) const
        {
            if (x >= xmin_ && x <= xmax_)
            {
                const double d = distro_->density(x)*fcn_(x)/norm_;
                assert(d >= 0.0);
                return d;
            }
            else
                return 0.0;
        }
 
        /** Cumulative distribution function */
        inline virtual double cdf(const double x) const
        {
            if (x <= xmin_)
                return 0.0;
            else if (x >= xmax_)
                return 1.0;
            else
            {
                GaussLegendreQuadrature glq(nIntegPoints_);
                return glq.integrate(DensityFunctor1D(*this),
                                     xmin_, x, nIntegIntervals_);
            }
        }
 
        /** 1 - cdf, avoiding subtractive cancellation */
        inline virtual double exceedance(const double x) const
        {
            if (x <= xmin_)
                return 1.0;
            else if (x >= xmax_)
                return 0.0;
            else
            {
                GaussLegendreQuadrature glq(nIntegPoints_);
                return glq.integrate(DensityFunctor1D(*this),
                                     x, xmax_, nIntegIntervals_);
            }
        }
 
        /** The quantile function */
        inline virtual double quantile(const double r1) const
        {
            if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
                "In npstat::ModulatedDistribution1D::quantile: "
                "cdf argument outside of [0, 1] interval");
            if (r1 == 0.0)
                return xmin_;
            else if (r1 == 1.0)
                return xmax_;
            else
            {
                const double tol = 2.0*std::numeric_limits<double>::epsilon();
                double q;
                const bool status = findRootUsingBisections(
                    CdfFunctor1D(*this), r1, xmin_, xmax_, tol, &q);
                assert(status);
                return q;
            }
        }
 
        /**
        // Read/write functions are not implemented so that we can
        // use this class with functors that do not support I/O
        */
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
 
        static inline const char* classname()
            {return "npstat::ModulatedDistribution1D";}
        static inline unsigned version() {return 1;}
 
    protected:
        inline virtual bool isEqual(const AbsDistribution1D& other) const
            {return false; /* Can't compare functors for equality */}
 
    private:
        ModulatedDistribution1D();
 
        AbsDistribution1D* distro_;
        Functor fcn_;
        double xmin_;
        double xmax_;
        double norm_;
        unsigned nIntegPoints_;
        unsigned nIntegIntervals_;
    };
 
    /** Convenience function for creating ModulatedDistribution1D objects */
    template<class Functor>
    inline ModulatedDistribution1D<Functor>
    make_ModulatedDistribution1D(const AbsDistribution1D& distro,
                                 const Functor& fcn,
                                 const double fcnXmin = -DBL_MAX,
                                 const double fcnXmax = DBL_MAX,
                                 const unsigned nIntegPoints = 1024,
                                 const unsigned nIntegIntervals = 1)
    {
        return ModulatedDistribution1D<Functor>(
            distro, fcn, fcnXmin, fcnXmax, nIntegPoints, nIntegIntervals);
    }
}
 
#endif // NPSTAT_MODULATEDDISTRIBUTION1D_HH_