DeltaMixture1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_DELTAMIXTURE1D_HH_
#define NPSTAT_DELTAMIXTURE1D_HH_
 
/*!
// \file DeltaMixture1D.hh
//
// \brief Weighted mixture of 1-d Dirac delta functions
//
// Author: I. Volobouev
//
// July 2019
*/
 
#include <vector>
#include <stdexcept>
 
#include "npstat/stat/AbsDistribution1D.hh"
#include "npstat/stat/Distribution1DFactory.hh"
 
namespace npstat {
    class DeltaMixture1D : public AbsScalableDistribution1D
    {
    public:
        /**
        // For this constructor, coordinates must be sorted
        // in the increasing order
        */
        template <typename Real1, typename Real2>
        inline DeltaMixture1D(const double location, const double scale,
                              const Real1* unscaledCoords,
                              const Real2* unscaledWeights,
                              const unsigned arrayLength)
            : AbsScalableDistribution1D(location, scale)
        {initialize(unscaledCoords, unscaledWeights, arrayLength);}
 
        inline virtual DeltaMixture1D* clone() const
            {return new DeltaMixture1D(*this);}
 
        inline virtual ~DeltaMixture1D() {}
 
        inline unsigned nComponents() const {return locations_.size();}
        inline double getLocation(const unsigned i) const {return locations_.at(i);}
        inline double getWeight(const unsigned i) const {return weights_.at(i);}
 
        /** Density integral on the intervals (a, b), [a, b], etc */
        double integral(double a, bool includePointAtA,
                        double b, bool includePointAtB) const;
 
        /** Distribution moment */
        long double moment(long double center, unsigned order) const;
 
        /** 
        // A collection of distribution moments. The array "moments"
        // should have at least "maxOrder+1" elements.
        */
        void moments(long double center, unsigned maxOrder,
                     long double* moments) const;
 
        /** 
        // A collection of distribution central moments. The array "moments"
        // should have at least "maxOrder+1" elements. On exit, moments[1]
        // will be set to the mean while all other moments will be central.
        */
        void centralMoments(unsigned maxOrder, long double* moments) const;
 
        // Methods needed for I/O
        virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream&) const;
 
        static inline const char* classname() {return "npstat::DeltaMixture1D";}
        static inline unsigned version() {return 1;}
        static DeltaMixture1D* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        virtual bool isEqual(const AbsDistribution1D&) const;
 
    private:
        friend class ScalableDistribution1DFactory<DeltaMixture1D>;
 
        DeltaMixture1D(double location, double scale,
                       const std::vector<double>& params);
        inline static int nParameters() {return -1;}
 
        template <typename Real1, typename Real2>
        void initialize(const Real1* unscaledCoords,
                        const Real2* unscaledWeights,
                        const unsigned arrayLength);
        void buildCdfTable();
 
        inline double unscaledDensity(double /* x */) const {return 0.0;}
        double unscaledCdf(double x) const;
        double unscaledExceedance(double x) const;
        double unscaledQuantile(double x) const;
 
        double weightSum(unsigned imin, unsigned imax) const;
 
        // Find all locations in the range
        void determineRange(double a, bool includeLeftBoundary,
                            double b, bool includeRightBoundary,
                            unsigned* imin, unsigned* imax) const;
 
        std::vector<double> locations_;
        std::vector<double> weights_;
        std::vector<double> cdf_;
        std::vector<double> exceedance_;
 
#ifdef SWIG
    public:
        inline DeltaMixture1D(const double location, const double scale,
                              const double* coords, unsigned lenCoords,
                              const double* weights, unsigned lenWeights)
            : AbsScalableDistribution1D(location, scale)
        {
            if (lenCoords != lenWeights) throw std::invalid_argument(
                "In npstat::DeltaMixture1D constructor: "
                "lengths of arrays of coordinates and weights are incompatible");
            initialize(coords, weights, lenCoords);
        }
 
        inline DeltaMixture1D(const double location, const double scale,
                              const double* coords, unsigned lenCoords)
            : AbsScalableDistribution1D(location, scale)
        {
            if (!lenCoords) throw std::invalid_argument(
                "In npstat::DeltaMixture1D constructor: "
                "empty array of coordinates");
            std::vector<double> weights(lenCoords, 1.0);
            initialize(coords, &weights[0], lenCoords);
        }
 
        inline double moment_2(const double center, const unsigned order) const
        {
            return moment(center, order);
        }
 
        inline std::vector<double> moments_2(const double center, const unsigned maxOrder) const
        {
            std::vector<long double> tmp(maxOrder+1U);
            moments(center, maxOrder, &tmp[0]);
            return std::vector<double>(tmp.begin(), tmp.end());
        }
 
        inline std::vector<double> centralMoments_2(const unsigned maxOrder) const
        {
            std::vector<long double> tmp(maxOrder+1U);
            centralMoments(maxOrder, &tmp[0]);
            return std::vector<double>(tmp.begin(), tmp.end());
        }
#endif
    };
}
 
#include "npstat/stat/DeltaMixture1D.icc"
 
#endif // NPSTAT_DELTAMIXTURE1D_HH_