AbsDistribution1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_ABSDISTRIBUTION1D_HH_
#define NPSTAT_ABSDISTRIBUTION1D_HH_
 
/*!
// \file AbsDistribution1D.hh
//
// \brief Interface definition for 1-d continuous statistical distributions
//
// Author: I. Volobouev
//
// November 2009
*/
 
#include <cassert>
#include <typeinfo>
#include <stdexcept>
 
#include "geners/ClassId.hh"
#include "npstat/nm/SimpleFunctors.hh"
#include "npstat/rng/AbsRandomGenerator.hh"
 
namespace npstat {
    /**
    // All classes derived from this base should have copy constructors
    // (possibly automatic).
    //
    // quantile(0.0) and quantile(1.0) should return the boundaries
    // of the support interval.
    */
    struct AbsDistribution1D
    {
        inline virtual ~AbsDistribution1D() {}
 
        /** Probability density */
        virtual double density(double x) const = 0;
 
        /** Cumulative distribution function */
        virtual double cdf(double x) const = 0;
 
        /**
        // 1 - cdf, known as "survival function" or "exceedance".
        // Implementations should avoid subtractive cancellation.
        */
        virtual double exceedance(double x) const = 0;
 
        /** The quantile function (inverse cdf) */
        virtual double quantile(double x) const = 0;
 
        /**
        // Derived classes should not implement "operator==", implement
        // "isEqual" instead
        */
        inline bool operator==(const AbsDistribution1D& r) const
            {return (typeid(*this) == typeid(r)) && this->isEqual(r);}
 
        /** Logical negation of operator== */
        inline bool operator!=(const AbsDistribution1D& r) const
            {return !(*this == r);}
 
        /** "Virtual copy constructor" */
        virtual AbsDistribution1D* clone() const = 0;
 
        /**
        // Random number generator according to the given distribution.
        // Should return the number of random points used up from the
        // generator.
        */
        inline virtual unsigned random(AbsRandomGenerator& g,
                                       double* generatedRandom) const
        {
            assert(generatedRandom);
            *generatedRandom = quantile(g());
            return 1U;
        }
 
        /**
        // Numerical moment calculation using Gauss-Legendre or
        // Fejer quadrature. If "useFejerQuadrature" argument is
        // set to "false", the number of integration points should
        // be supported by the "GaussLegendreQuadrature" class.
        //
        // No attempt is made to determine whether the moment actually
        // exists. If it does not, the function might return Inf or
        // fail in some other ways. This method might also be highly
        // inaccurate for distributions with singularities or for
        // distributions that have a lot of weight near support boundaries.
        //
        // For better results, it might be useful to calculate the
        // 0th order moment (i.e., the normalization factor for the
        // quadrature scheme chosen) and then divide everything by it.
        */
        long double empiricalMoment(long double center, unsigned order,
                                    unsigned nIntegrationPoints = 1024U,
                                    bool useFejerQuadrature = false) const;
 
        /**
        // Numerical calculation of the expectation value of some functor,
        // using Gauss-Legendre or Fejer quadrature. If "useFejerQuadrature"
        // argument is set to "false", the number of integration points
        // should be supported by the "GaussLegendreQuadrature" class.
        */
        template <class Functor>
        double expectation(const Functor& fcn,
                           unsigned nIntegrationPoints = 1024U,
                           bool useFejerQuadrature = false) const;
 
        //@{
        /** Prototype needed for I/O */
        virtual gs::ClassId classId() const = 0;
        inline virtual bool write(std::ostream&) const {return false;}
        //@}
 
        static inline const char* classname()
            {return "npstat::AbsDistribution1D";}
        static inline unsigned version() {return 1;}
        static AbsDistribution1D* read(const gs::ClassId& id, std::istream&);
 
    protected:
        /** Comparison for equality. To be implemented by derived classes. */
        virtual bool isEqual(const AbsDistribution1D&) const = 0;
 
    private:
        template <class Funct>
        class ExpectationFunctor : public Functor1<double, double>
        {
        public:
            inline ExpectationFunctor(const Funct& fcn,
                                      const AbsDistribution1D& distro)
                : fcn_(fcn), distro_(distro) {}
 
            inline virtual ~ExpectationFunctor() {}
            
            inline double operator()(const double& x) const
                {return fcn_(x)*distro_.density(x);}
        private:
            const Funct& fcn_;
            const AbsDistribution1D& distro_;
        };
 
#ifdef SWIG
    public:
        inline double random2(AbsRandomGenerator& g) const
        {
            double d;
            this->random(g, &d);
            return d;
        }
#endif // SWIG
    };
 
    class UnscaledDensityFunctor1D;
    class UnscaledCdfFunctor1D;
 
    /**
    // This base class is used to model distributions which have
    // "trivial" behavior w.r.t. location and scale parameters. That is,
    // if the distribution density is g(x), the density with location
    // parameter "mu" and scale parameter "sigma" is g((x - mu)/sigma)/sigma.
    //
    // All distributions with infinite support on both sides fall in this
    // category, as well as a large number of distributions with finite
    // support.
    */
    class AbsScalableDistribution1D : public AbsDistribution1D
    {
    public:
        /** Location and scale parameters must be provided */
        inline AbsScalableDistribution1D(const double location,
                                         const double scale) :
            location_(location), scale_(scale)
        {
            if (scale_ <= 0.0) throw std::invalid_argument(
                "In npstat::AbsScalableDistribution1D constructor: "
                "scale parameter must be positive");
        }
 
        inline virtual ~AbsScalableDistribution1D() {}
 
        /** Get the location parameter */
        inline double location() const {return location_;}
 
        /** Get the scale parameter */
        inline double scale() const {return scale_;}
 
        /** Set the location parameter */
        inline void setLocation(const double v) {location_ = v;}
 
        /** Set the scale parameter */
        inline void setScale(const double v)
        {
            if (v <= 0.0) throw std::invalid_argument(
                "In npstat::AbsScalableDistribution1D::setScale: "
                "scale parameter must be positive");
            scale_ = v;
        }
 
        //@{
        /** Method overriden from the AbsDistribution1D base class */
        inline double density(const double x) const
            {return unscaledDensity((x - location_)/scale_)/scale_;}
 
        inline double cdf(const double x) const
            {return unscaledCdf((x - location_)/scale_);}
 
        inline double exceedance(const double x) const
            {return unscaledExceedance((x - location_)/scale_);}
 
        inline double quantile(const double x) const
            {return scale_*unscaledQuantile(x) + location_;}
 
        inline virtual unsigned random(AbsRandomGenerator& g,
                                       double* generatedRandom) const
        {
            const unsigned count = unscaledRandom(g, generatedRandom);
            *generatedRandom *= scale_;
            *generatedRandom += location_;
            return count;
        }
        //@}
 
        /** "Virtual copy constructor" */
        virtual AbsScalableDistribution1D* clone() const = 0;
 
        //@{
        /** Method related to "geners" I/O */
        virtual gs::ClassId classId() const = 0;
        virtual bool write(std::ostream& os) const;
        //@}
 
        /**
        // Pseudo-read function for I/O. Note that this class is abstract,
        // so its instance can not be created.
        */
        static bool read(std::istream& is, double* location, double* scale);
 
    protected:
        /**
        // Derived classes should override the following method as long as
        // they have at least one additional data member. Don't forget to
        // call "isEqual" of the base class inside the derived classes.
        */
        virtual bool isEqual(const AbsDistribution1D& other) const
        {
            const AbsScalableDistribution1D& r = 
                static_cast<const AbsScalableDistribution1D&>(other);
            return location_ == r.location_ && scale_ == r.scale_;
        }
 
    private:
        friend class UnscaledDensityFunctor1D;
        friend class UnscaledCdfFunctor1D;
 
        AbsScalableDistribution1D();
 
        virtual double unscaledDensity(double x) const = 0;
        virtual double unscaledCdf(double x) const = 0;
        virtual double unscaledExceedance(double x) const = 0;
        virtual double unscaledQuantile(double x) const = 0;
        inline virtual unsigned unscaledRandom(AbsRandomGenerator& g,
                                               double* generatedRandom) const
        {
            assert(generatedRandom);
            *generatedRandom = unscaledQuantile(g());
            return 1U;
        }
 
        double location_;
        double scale_;
    };
 
    /**
    // A functor for the density function of the given 1-d distribution.
    // For this and subsequent functors, the distribution is not copied,
    // only a reference is used. It is a responsibility of the user to
    // make sure that the lifetime of the distribution object exceeds
    // the lifetime of the functor.
    */
    class DensityFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit DensityFunctor1D(const AbsDistribution1D& fcn,
                                         const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~DensityFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return norm_*fcn_.density(a);}
 
    private:
        DensityFunctor1D();
        const AbsDistribution1D& fcn_;
        const double norm_;
    };
 
    /** A functor for the density squared of the given 1-d distribution */
    class DensitySquaredFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit DensitySquaredFunctor1D(const AbsDistribution1D& fcn,
                                                const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~DensitySquaredFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
        {
            const double d = fcn_.density(a);
            return norm_*norm_*d*d;
        }
 
    private:
        DensitySquaredFunctor1D();
        const AbsDistribution1D& fcn_;
        const double norm_;
    };
 
    /** 
    // A functor for the cumulative distribution function of
    // the given 1-d distribution
    */
    class CdfFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit CdfFunctor1D(const AbsDistribution1D& fcn,
                                     const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~CdfFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return norm_*fcn_.cdf(a);}
 
    private:
        CdfFunctor1D();
        const AbsDistribution1D& fcn_;
        const double norm_;
    };
 
    /** A functor for the exceedance function of the given 1-d distribution */
    class ExceedanceFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit ExceedanceFunctor1D(const AbsDistribution1D& fcn,
                                            const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~ExceedanceFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return norm_*fcn_.exceedance(a);}
 
    private:
        ExceedanceFunctor1D();
        const AbsDistribution1D& fcn_;
        const double norm_;
    };
 
    /** A functor for the quantile function of the given 1-d distribution */
    class QuantileFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit QuantileFunctor1D(const AbsDistribution1D& fcn,
                                          const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~QuantileFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return fcn_.quantile(a/norm_);}
 
    private:
        QuantileFunctor1D();
        const AbsDistribution1D& fcn_;
        const double norm_;
    };
 
    /** A functor for the unscaled density */
    class UnscaledDensityFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit UnscaledDensityFunctor1D(const AbsScalableDistribution1D& fcn,
                                                 const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~UnscaledDensityFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return norm_*fcn_.unscaledDensity(a);}
 
    private:
        UnscaledDensityFunctor1D();
        const AbsScalableDistribution1D& fcn_;
        const double norm_;
    };
 
    /** A functor for the unscaled cdf */
    class UnscaledCdfFunctor1D : public Functor1<double, double>
    {
    public:
        inline explicit UnscaledCdfFunctor1D(const AbsScalableDistribution1D& fcn,
                                             const double normfactor=1.0)
            : fcn_(fcn), norm_(normfactor) {}
 
        inline virtual ~UnscaledCdfFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
            {return norm_*fcn_.unscaledCdf(a);}
 
    private:
        UnscaledCdfFunctor1D();
        const AbsScalableDistribution1D& fcn_;
        const double norm_;
    };
 
    /**
    // A functor that gives the ratio between the density and its
    // maximum value. Can be used to sample the distributions by
    // acceptance-rejection. Does not copy the distribution, so
    // the functor lifetime should not exceed the lifetime of the
    // distribution.
    */
    class AcceptanceFunctor1D : public Functor1<double, double>
    {
    public:
        inline AcceptanceFunctor1D(const AbsDistribution1D& fcn,
                                   const double maxDensityValue)
            : fcn_(fcn), max_(maxDensityValue)
        {
            if (max_ <= 0.0) throw std::invalid_argument(
                "In npstat::AcceptanceFunctor1D constructor: "
                "maximum density value must be positive");
        }
 
        inline virtual ~AcceptanceFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
        {
            const double ratio = fcn_.density(a)/max_;
            if (ratio > 1.0) throw std::runtime_error(
                "In npstat::AcceptanceFunctor1D::operator(): "
                "the maximum density parameter is too small");
            return ratio;
        }
 
    private:
        AcceptanceFunctor1D();
        const AbsDistribution1D& fcn_;
        const double max_;
    };
 
    /**
    // A functor that gives the ratio between the minimum density
    // value and the density itself. For use in sampling inverse
    // densities.
    */
    class InverseAcceptanceFunctor1D : public Functor1<double, double>
    {
    public:
        inline InverseAcceptanceFunctor1D(const AbsDistribution1D& fcn,
                                          const double minDensityValue)
            : fcn_(fcn), min_(minDensityValue)
        {
            if (min_ <= 0.0) throw std::invalid_argument(
                "In npstat::AcceptanceFunctor1D constructor: "
                "minimum density value must be positive");
        }
 
        inline virtual ~InverseAcceptanceFunctor1D() {}
 
        inline virtual double operator()(const double& a) const
        {
            const double d = fcn_.density(a);
            assert(d > 0.0);
            const double ratio = min_/d;
            if (ratio > 1.0) throw std::runtime_error(
                "In npstat::InverseAcceptanceFunctor1D::operator(): "
                "the minimum density parameter is too large");
            return ratio;
        }
 
    private:
        InverseAcceptanceFunctor1D();
        const AbsDistribution1D& fcn_;
        const double min_;
    };
}
 
#include "npstat/stat/AbsDistribution1D.icc"
 
#endif // NPSTAT_ABSDISTRIBUTION1D_HH_