DistributionsND.hh

Go to the documentation of this file.

#ifndef NPSTAT_DISTRIBUTIONSND_HH_
#define NPSTAT_DISTRIBUTIONSND_HH_
 
/*!
// \file DistributionsND.hh
//
// \brief A number of useful multivariate continuous statistical distributions
//
// Author: I. Volobouev
//
// March 2010
*/
 
#include <cassert>
#include <stdexcept>
 
#include "npstat/nm/Matrix.hh"
#include "npstat/nm/MathUtils.hh"
 
#include "npstat/stat/AbsDistributionND.hh"
#include "npstat/stat/GridRandomizer.hh"
#include "npstat/stat/Distributions1D.hh"
 
namespace npstat {
    /**
    // A generic product distribution. The argument distributions
    // will be cloned internally.
    */
    class ProductDistributionND : public AbsDistributionND
    {
    public:
        ProductDistributionND(const AbsDistribution1D** marginals,
                              unsigned nMarginals);
 
        // Convenience constructors. They are also good for use from Python.
        ProductDistributionND(const AbsDistribution1D& m0);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4,
                              const AbsDistribution1D& m5);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4,
                              const AbsDistribution1D& m5,
                              const AbsDistribution1D& m6);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4,
                              const AbsDistribution1D& m5,
                              const AbsDistribution1D& m6,
                              const AbsDistribution1D& m7);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4,
                              const AbsDistribution1D& m5,
                              const AbsDistribution1D& m6,
                              const AbsDistribution1D& m7,
                              const AbsDistribution1D& m8);
 
        ProductDistributionND(const AbsDistribution1D& m0,
                              const AbsDistribution1D& m1,
                              const AbsDistribution1D& m2,
                              const AbsDistribution1D& m3,
                              const AbsDistribution1D& m4,
                              const AbsDistribution1D& m5,
                              const AbsDistribution1D& m6,
                              const AbsDistribution1D& m7,
                              const AbsDistribution1D& m8,
                              const AbsDistribution1D& m9);
 
        ProductDistributionND(const ProductDistributionND& r);
        ProductDistributionND& operator=(const ProductDistributionND& r);
 
        virtual ~ProductDistributionND();
 
        inline virtual ProductDistributionND* clone() const
            {return new ProductDistributionND(*this);}
 
        double density(const double* x, unsigned dim) const;
        void unitMap(const double* rnd, unsigned bufLen, double* x) const;
        inline bool mappedByQuantiles() const {return true;}
 
        /** 
        // Check whether all marginals are instances of
        // AbsScalableDistribution1D
        */
        bool isScalable() const;
 
        /** Get the marginal distribution with the given index */
        inline AbsDistribution1D* getMarginal(const unsigned i)
            {return marginals_.at(i);}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname()
            {return "npstat::ProductDistributionND";}
        static inline unsigned version() {return 1;}
        static ProductDistributionND* read(const gs::ClassId& id,
                                           std::istream& in);
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        inline explicit ProductDistributionND(const unsigned n)
            : AbsDistributionND(n) {}
        void cleanup();
 
        std::vector<AbsDistribution1D*> marginals_;
    };
 
    /**
    // Multivariate Gaussian distribution parameterized
    // by its covariance matrix
    */
    class GaussND : public AbsDistributionND
    {
    public:
        GaussND(const double* location1, unsigned dim1,
                const Matrix<double>& covmat);
 
        inline virtual ~GaussND() {}
        inline virtual GaussND* clone() const {return new GaussND(*this);}
 
        double density(const double* x, unsigned dim) const;
        void unitMap(const double* rnd, unsigned bufLen, double* x) const;
        inline bool mappedByQuantiles() const {return false;}
 
        inline double getLocation(const unsigned index) const
            {return location_.at(index);}
        inline const Matrix<double>& getCovMat() const {return covmat_;}
        const Matrix<double>& getCorrMat() const;
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname() {return "npstat::GaussND";}
        static inline unsigned version() {return 1;}
        static GaussND* read(const gs::ClassId& id, std::istream&);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        void initialize();
 
        std::vector<double> location_;
        mutable std::vector<double> buf_;
        mutable Matrix<double> sqrt_;
        Matrix<double> covmat_;
        Matrix<double> invsqrt_;
        mutable Matrix<double> corrmat_;
        double norm_;
 
#ifdef SWIG
    public:
        inline GaussND(const double* location1, unsigned dim1,
                       const double* covmat, unsigned nrows, unsigned ncols)
            : AbsDistributionND(dim1), location_(location1, location1+dim1),
              buf_(dim1), covmat_(nrows, ncols, covmat)
        {
            assert(location1);
            assert(covmat);
            covmat_ = (covmat_ + covmat_.T())/2.0;
            initialize();
        }
#endif // SWIG
    };
 
    /** Multivariate uniform distribution */
    class UniformND : public AbsScalableDistributionND
    {
    public:
        inline UniformND(const double* location, const double* scale,
                         const unsigned dim)
            : AbsScalableDistributionND(location, scale, dim) {}
 
        inline virtual ~UniformND() {}
 
        inline virtual UniformND* clone() const {return new UniformND(*this);}
        inline bool mappedByQuantiles() const {return true;}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& of) const;
 
        static inline const char* classname() {return "npstat::UniformND";}
        static inline unsigned version() {return 1;}
        static UniformND* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        virtual bool isEqual(const AbsDistributionND& r) const
            {return AbsScalableDistributionND::isEqual(r);}
 
    private:
        double unscaledDensity(const double* x) const;
        void unscaledUnitMap(const double* rnd, unsigned bufLen,
                             double* x) const;
#ifdef SWIG
    public:
        inline UniformND(const double* location1, unsigned dim1,
                         const double* scale2, unsigned dim2)
            : AbsScalableDistributionND(location1, scale2, dim2)
        {
            if (dim1 != dim2) throw std::invalid_argument(
                "In npstat::UniformND constructor: "
                "sizes of locations and scales are incompatible");
        }
#endif // SWIG
    };
 
    /**
    // Multivariate symmetric distribution proportional to pow(1.0 - R^2, n)
    // with the support region defined by 0 <= R <= 1, where R is the
    // distance (after scaling) to the distribution location point.
    // Can be scaled separately in each dimension.
    */
    class ScalableSymmetricBetaND : public AbsScalableDistributionND
    {
    public:
        ScalableSymmetricBetaND(const double* location, const double* scale,
                                unsigned dim, double n);
        inline virtual ScalableSymmetricBetaND* clone() const
            {return new ScalableSymmetricBetaND(*this);}
 
        inline virtual ~ScalableSymmetricBetaND() {}
 
        inline bool mappedByQuantiles() const {return false;}
        inline double power() const {return power_;}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname()
            {return "npstat::ScalableSymmetricBetaND";}
        static inline unsigned version() {return 1;}
        static ScalableSymmetricBetaND* read(const gs::ClassId& id, std::istream&);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        void initialize();
        double unscaledDensity(const double* x) const;
        void unscaledUnitMap(const double* rnd, unsigned bufLen,
                             double* x) const;
 
        double radialQuantile(double cdf) const;
 
        double power_;
        double norm_;
 
#ifdef SWIG
    public:
        inline ScalableSymmetricBetaND(const double* location1, unsigned dim1,
                                       const double* scale2, unsigned dim2,
                                       const double n)
            : AbsScalableDistributionND(location1, scale2, dim2), power_(n)
        {
            if (dim1 != dim2) throw std::invalid_argument(
                "In npstat::ScalableSymmetricBetaND constructor: "
                "sizes of locations and scales are incompatible");
            initialize();
        }
#endif // SWIG
    };
 
    /**
    // Multidimensional Huber-like function. Gaussian near the
    // coordinate origin and exponential when radius is above
    // certain threshold. "tailWeight" parameter is the weight
    // of the exponential tail. It must be between 0 and 1.
    */
    class ScalableHuberND : public AbsScalableDistributionND
    {
    public:
        ScalableHuberND(const double* location, const double* scale,
                        unsigned dim, double tailWeight);
        inline virtual ScalableHuberND* clone() const
            {return new ScalableHuberND(*this);}
 
        inline virtual ~ScalableHuberND() {}
 
        inline bool mappedByQuantiles() const {return false;}
 
        inline double tailWeight() const {return tWeight_;}
        inline double transition() const {return a_;}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname()
            {return "npstat::ScalableHuberND";}
        static inline unsigned version() {return 1;}
        static ScalableHuberND* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        void initialize();
        double unscaledDensity(const double* x) const;
        void unscaledUnitMap(const double* rnd, unsigned bufLen,
                             double* x) const;
 
        // Weight of the Gaussian part
        double norm1(double a) const;
 
        // Weight of the exponential part
        double norm2(double a) const;
 
        double radialQuantile(double cdf) const;
 
        double tWeight_;
        double a_;
        double normfactor_;
        double const1_;
        double const2_;
 
#ifdef SWIG
    public:
        inline ScalableHuberND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const double alpha)
            : AbsScalableDistributionND(location1, scale2, dim2), tWeight_(alpha)
        {
            if (dim1 != dim2) throw std::invalid_argument(
                "In npstat::ScalableHuberND constructor: "
                "sizes of locations and scales are incompatible");
            initialize();
        }
#endif // SWIG
    };
 
    /** Product of symmetric beta distributions in each dimension */
    struct ProductSymmetricBetaND : 
        public HomogeneousProductDistroND<SymmetricBeta1D>
    {
        ProductSymmetricBetaND(const double* location, const double* scale,
                               unsigned dim, double power);
        inline virtual ProductSymmetricBetaND* clone() const
            {return new ProductSymmetricBetaND(*this);}
 
        inline virtual ~ProductSymmetricBetaND() {}
 
        inline double power() const {return marginals_[0].power();}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname()
            {return "npstat::ProductSymmetricBetaND";}
        static inline unsigned version() {return 1;}
        static ProductSymmetricBetaND* read(const gs::ClassId& id,
                                            std::istream& in);
#ifdef SWIG
    public:
        inline ProductSymmetricBetaND(const double* location1, unsigned dim1,
                                      const double* scale2, unsigned dim2,
                                      const double power)
            : HomogeneousProductDistroND<SymmetricBeta1D>(dim2)
        {
            if (dim1 != dim2) throw std::invalid_argument(
                "In npstat::ProductSymmetricBetaND constructor: "
                "sizes of locations and scales are incompatible");
            assert(location1);
            assert(scale2);
            for (unsigned i=0; i<dim2; ++i)
                marginals_.emplace_back(location1[i], scale2[i], power);
        }
#endif // SWIG
    };
 
    /** Distribution defined by the interpolation table of its radial profile */
    class RadialProfileND : public AbsScalableDistributionND
    {
    public:
        /**
        // "location", "scale", and "dim" are the usual multivariate
        // density parameters which define the distribution location,
        // scale in each dimension, and the number of dimensions.
        //
        // The "data" array gives density values, (d Prob)/(dx_1 ... dx_n),
        // at equidistant radial intervals. data[0] is density at r = 0.0,
        // and data[dataLen-1] is density at r = 1.0. If "dataLen" is
        // less than 2, uniform distribution will be created. The total
        // probability content will be automatically normalized to 1, so
        // the "data" values do not have to be normalized by the user.
        // Internally, the profile is kept in double precision.
        //
        // "interpolationDegree" is the order of the radial profile
        // interpolation polynomial. It must be less than 4 and less
        // than "dataLen".
        */
        template <typename Real>
        RadialProfileND(const double* location, const double* scale,
                        unsigned dim,
                        const Real* data, unsigned dataLen,
                        unsigned interpolationDegree);
 
        inline virtual RadialProfileND* clone() const
            {return new RadialProfileND(*this);}
 
        inline virtual ~RadialProfileND() {}
 
        inline bool mappedByQuantiles() const {return false;}
 
        inline unsigned interpolationDegree() const {return deg_;}
        inline unsigned profileLength() const {return len_;}
        inline const double* profileData() const {return &profile_[0];}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& of) const;
 
        static inline const char* classname()
            {return "npstat::RadialProfileND";}
        static inline unsigned version() {return 1;}
        static RadialProfileND* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        inline RadialProfileND(const double* location, const double* scale,
                               const unsigned dim)
            : AbsScalableDistributionND(location, scale, dim) {}
 
        template <typename Real>
        void initialize(const Real* data, unsigned dataLen);
 
        double unscaledDensity(const double* x) const;
        void unscaledUnitMap(const double* rnd, unsigned bufLen,
                             double* x) const;
        void calculateQuadratureCoeffs();
        void normalize();
        double radialPofile(double r) const;
        double intervalInteg(unsigned intervalNumber,
                             double stepFraction=1.0) const;
        double intervalIntegLowD(unsigned intervalNumber,
                                 double stepFraction=1.0) const;
        double radialQuantile(double cdf) const;
 
        std::vector<double> profile_;
        std::vector<double> cdf_;
        double ndSphereVol_;
        double step_;
        unsigned len_;
        unsigned deg_;
 
        // Coefficients for Gaussian quadratures
        double xg_[4], wg_[4];
 
#ifdef SWIG
    public:
        inline RadialProfileND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const double* data1, unsigned dataLen1,
                               unsigned interpolationDegree)
            : AbsScalableDistributionND(location1, scale2, dim2),
              ndSphereVol_(ndUnitSphereVolume(dim2)),
              len_(dataLen1), deg_(interpolationDegree)
        {
            if (dim1 != dim2) throw std::invalid_argument(
                "In npstat::RadialProfileND constructor: "
                "sizes of locations and scales are incompatible");
            this->initialize(data1, dataLen1);
        }
#endif // SWIG
    };
 
    /**
    // Distribution defined by an interpolation table inside
    // the unit box (which can be shifted and scaled). All grid
    // points are inside the box, in the bin centers. Currently,
    // interpolationDegree could be only 0 (faster, no interpolation)
    // or 1 (multilinear interpolation).
    */
    class BinnedDensityND : public AbsScalableDistributionND
    {
    public:
        template <typename Num1, unsigned Len1, unsigned Dim1>
        BinnedDensityND(const double* location, const double* scale,
                        unsigned locationAndScaleLength,
                        const ArrayND<Num1,Len1,Dim1>& histogram,
                        const unsigned interpolationDegree);
 
        inline virtual BinnedDensityND* clone() const
            {return new BinnedDensityND(*this);}
 
        inline virtual ~BinnedDensityND() {}
 
        inline bool mappedByQuantiles() const {return true;}
 
        inline const ArrayND<double>& gridData() const
            {return randomizer_.gridData();}
 
        inline unsigned interpolationDegree() const
            {return randomizer_.interpolationDegree();}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname() {return "npstat::BinnedDensityND";}
        static inline unsigned version() {return 1;}
        static BinnedDensityND* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        inline double unscaledDensity(const double* x) const
            {return randomizer_.density(x, dim_);}
 
        inline void unscaledUnitMap(const double* rnd,
                                    const unsigned bufLen, double* x) const
            {randomizer_.generate(rnd, bufLen, x);}
 
        GridRandomizer randomizer_;
 
#ifdef SWIG
    public:
        inline BinnedDensityND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const ArrayND<double>& histogram,
                               const unsigned interpolationDegree)
            :  AbsScalableDistributionND(location1, scale2, dim2),
               randomizer_(histogram, BoxND<double>::unitBox(histogram.rank()),
                           interpolationDegree)
        {
            if (dim1 != dim2 || dim1 != histogram.rank())
                throw std::invalid_argument(
                    "In npstat::BinnedDensityND constructor: "
                    "incompatible argument dimensions");
        }
        inline BinnedDensityND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const ArrayND<float>& histogram,
                               const unsigned interpolationDegree)
            :  AbsScalableDistributionND(location1, scale2, dim2),
               randomizer_(histogram, BoxND<double>::unitBox(histogram.rank()),
                           interpolationDegree)
        {
            if (dim1 != dim2 || dim1 != histogram.rank())
                throw std::invalid_argument(
                    "In npstat::BinnedDensityND constructor: "
                    "incompatible argument dimensions");
        }
        inline BinnedDensityND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const ArrayND<int>& histogram,
                               const unsigned interpolationDegree)
            :  AbsScalableDistributionND(location1, scale2, dim2),
               randomizer_(histogram, BoxND<double>::unitBox(histogram.rank()),
                           interpolationDegree)
        {
            if (dim1 != dim2 || dim1 != histogram.rank())
                throw std::invalid_argument(
                    "In npstat::BinnedDensityND constructor: "
                    "incompatible argument dimensions");
        }
        inline BinnedDensityND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const ArrayND<long>& histogram,
                               const unsigned interpolationDegree)
            :  AbsScalableDistributionND(location1, scale2, dim2),
               randomizer_(histogram, BoxND<double>::unitBox(histogram.rank()),
                           interpolationDegree)
        {
            if (dim1 != dim2 || dim1 != histogram.rank())
                throw std::invalid_argument(
                    "In npstat::BinnedDensityND constructor: "
                    "incompatible argument dimensions");
        }
        inline BinnedDensityND(const double* location1, unsigned dim1,
                               const double* scale2, unsigned dim2,
                               const ArrayND<unsigned char>& histogram,
                               const unsigned interpolationDegree)
            :  AbsScalableDistributionND(location1, scale2, dim2),
               randomizer_(histogram, BoxND<double>::unitBox(histogram.rank()),
                           interpolationDegree)
        {
            if (dim1 != dim2 || dim1 != histogram.rank())
                throw std::invalid_argument(
                    "In npstat::BinnedDensityND constructor: "
                    "incompatible argument dimensions");
        }
#endif // SWIG
    };
 
    /**
    // Linear transform of another multivariate distribution. If the other
    // distribution is a function of x, this distribution is obtained by
    // substituting S^(-1)*(x - shift) in place of x, where S is some
    // non-degenerate matrix.
    */
    class LinTransformedDistroND : public AbsDistributionND
    {
    public:
        inline LinTransformedDistroND(const AbsDistributionND& d,
                                      const double* shift, unsigned i_dim,
                                      const Matrix<double>& S)
            : AbsDistributionND(i_dim), distro_(d.clone()),
              shift_(shift, shift+i_dim), buf_(2U*i_dim), S_(S)
        {
            assert(shift);
            initialize();
        }
 
        LinTransformedDistroND(const LinTransformedDistroND&);
        LinTransformedDistroND& operator=(const LinTransformedDistroND&);
        inline virtual LinTransformedDistroND* clone() const
            {return new LinTransformedDistroND(*this);}
 
        inline virtual ~LinTransformedDistroND() {delete distro_;}
 
        double density(const double* x, unsigned dim) const;
        void unitMap(const double* rnd, unsigned bufLen, double* x) const;
        inline bool mappedByQuantiles() const {return false;}
 
        inline double getShift(const unsigned index) const
            {return shift_.at(index);}
        inline const Matrix<double>& getS() const {return S_;}
        inline const Matrix<double>& getInverseS() const {return invS_;}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname() {return "npstat::LinTransformedDistroND";}
        static inline unsigned version() {return 1;}
        static LinTransformedDistroND* read(const gs::ClassId& id, std::istream&);
 
    protected:
        virtual bool isEqual(const AbsDistributionND&) const;
 
    private:
        void initialize();
 
        AbsDistributionND* distro_;
        std::vector<double> shift_;
        mutable std::vector<double> buf_;
        Matrix<double> S_;
        Matrix<double> invS_;
        double norm_;
 
#ifdef SWIG
    public:
        inline LinTransformedDistroND(const AbsDistributionND& d,
                                      const double* shift, unsigned i_dim,
                                      const double* covmat, unsigned nrows, unsigned ncols)
            : AbsDistributionND(i_dim), distro_(d.clone()), shift_(shift, shift+i_dim),
              buf_(2U*i_dim), S_(nrows, ncols, covmat)
        {
            assert(shift);
            assert(covmat);
            initialize();
        }
#endif // SWIG
    };
}
 
#include "npstat/stat/DistributionsND.icc"
 
#endif // NPSTAT_DISTRIBUTIONSND_HH_