OSDE1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_OSDE1D_HH_
#define NPSTAT_OSDE1D_HH_
 
/*!
// \file OSDE1D.hh
//
// \brief Orthogonal Series Density Estimation (OSDE) in one dimension
//
// Author: I. Volobouev
//
// May 2017
*/
 
#include <utility>
#include <vector>
 
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
#include "npstat/stat/AbsDistribution1D.hh"
 
#ifdef SWIG
#include <cassert>
#include <stdexcept>
 
#include "npstat/wrap/arrayNDToNumpy.hh"
#include "npstat/wrap/numpyArrayUtils.hh"
#endif // SWIG
 
namespace npstat {
    class OSDE1D
    {
    public:
        /**
        // Main constructor. The arguments are as follows:
        //
        //  poly1d -- classical orthogonal polynomial system to use
        //
        //  xmin   -- lower bound of the density support region
        //
        //  xmax   -- upper bound of the density support region
        //
        // Note that the density estimate will be automatically
        // normalized for Legendre polynomials only.
        */
        OSDE1D(const AbsClassicalOrthoPoly1D& poly1d,
               double xmin, double xmax);
 
        OSDE1D(const OSDE1D&);
        OSDE1D& operator=(const OSDE1D&);
 
        ~OSDE1D();
 
        //@{
        /** Basic inspection of object properties */
        inline const AbsClassicalOrthoPoly1D& clpoly() const {return *poly_;}
        inline double xmin() const {return xmin_;}
        inline double xmax() const {return xmax_;}
        double weight(double x) const;
        //@}
 
        /** Polynomial values */
        inline double poly(const unsigned deg, const double x) const
            {return poly_->poly(deg, (x - shift_)/scale_);}
 
        /**
        // A pair of polynomial values. Faster than
        // calling "poly" two times.
        */
        std::pair<double,double> twopoly(
            unsigned deg1, unsigned deg2, double x) const;
 
        /** Polynomial series representing the fitted density */
        double series(const double* coeffs, unsigned maxdeg, double x) const;
 
        /** Integral of the series */
        inline double integrateSeries(const double* coeffs,
                                      const unsigned maxdeg) const
            {return scale_*poly_->integrateSeries(coeffs, maxdeg);}
 
        /** Unweighted polynomial integral */
        inline double integratePoly(const unsigned degree,
                                    const unsigned power) const
            {return scale_*poly_->integratePoly(degree, power);}
 
        /** Unweighted integral of a product of multiple polynomials */
        inline double jointIntegral(
            const unsigned* degrees, unsigned nDegrees) const
            {return scale_*poly_->jointIntegral(degrees, nDegrees);}
 
        /**
        // Matrix of unweighted pairwise integrals. The returned
        // matrix will be dimensioned (maxdeg + 1) x (maxdeg + 1).
        */
        Matrix<double> tMatrix(unsigned maxdeg) const;
 
        /**
        // Build the coefficients of the orthogonal polynomial series
        // for the given function. The length of the array "coeffs"
        // should be at least maxdeg + 1. Note that the coefficients
        // are returned in the order of increasing degree (same order
        // is used by the "series" function).
        */
        template <class Functor, class Quadrature>
        void calculateCoeffs(const Functor& fcn, const Quadrature& quad,
                             double* coeffs, unsigned maxdeg) const;
 
        /**
        // Build the coefficients of the orthogonal polynomial series
        // for the given density using Gauss-Legendre quadratures.
        // The length of the array "coeffs" should be at least maxdeg + 1.
        */
        void densityCoeffs(const AbsDistribution1D& distro,
                           unsigned nIntegrationPoints,
                           double* coeffs, unsigned maxdeg) const;
 
        /**
        // Expected variances of the coefficients obtained previously
        // with the "densityCoeffs" method. The length of the array
        // "variances" should be at least maxdeg + 1.
        */
        void densityCoeffsVariance(const AbsDistribution1D& distro,
                                   unsigned nIntegrationPoints,
                                   const double* coeffs, unsigned maxdeg,
                                   double expectedSampleSize,
                                   double* variances) const;
 
        /**
        // Expected covariances of the coefficients obtained previously
        // with the "densityCoeffs" method.
        */
        npstat::Matrix<double>
        densityCoeffsCovariance(const AbsDistribution1D& distro,
                                unsigned nIntegrationPoints,
                                const double* coeffs, unsigned maxdeg,
                                double expectedSampleSize) const;
 
        /**
        // Build the coefficients of the orthogonal polynomial series
        // for the given sample of points (empirical density function).
        // The length of the array "coeffs" should be at least maxdeg + 1.
        // Note that the coefficients are returned in the order of
        // increasing degree (same order is used by the "series" function).
        */
        template <class Numeric>
        void sampleCoeffs(const Numeric* coords, unsigned long lenCoords,
                          double* coeffs, unsigned maxdeg) const;
 
        /**
        // Estimate variances of the coefficients obtained previously
        // with the "sampleCoeffs" method.
        */
        template <class Numeric>
        void sampleCoeffsVariance(const Numeric* coords, unsigned long lenCoords,
                                  const double* coeffs, unsigned maxdeg,
                                  double* variances) const;
 
        /**
        // Estimate covariances of the coefficients obtained previously
        // with the "sampleCoeffs" method.
        */
        template <class Numeric>
        npstat::Matrix<double>
        sampleCoeffsCovariance(const Numeric* coords, unsigned long lenCoords,
                               const double* coeffs, unsigned maxdeg) const;
 
        /**
        // Build the coefficients of the orthogonal polynomial series for
        // the given sample of weighted points (empirical density function).
        // The first element of the pair will be treated as the coordinate
        // and the second element will be treated as weight. Weights must
        // not be negative. The length of the array "coeffs" should be
        // at least maxdeg + 1. Note that the coefficients are returned
        // in the order of increasing degree (same order is used by the
        // "series" function).
        */
        template <class Pair>
        void weightedSampleCoeffs(const Pair* points, unsigned long numPoints,
                                  double* coeffs, unsigned maxdeg) const;
 
        /**
        // Estimate variances of the coefficients obtained previously
        // with the "weightedSampleCoeffs" method. This code assumes that
        // weights and coordinates are statistically independent from
        // each other.
        */
        template <class Pair>
        void weightedSampleCoeffsVariance(const Pair* points, unsigned long nPoints,
                                          const double* coeffs, unsigned maxdeg,
                                          double* variances) const;
 
        /**
        // Estimate covariances of the coefficients obtained previously
        // with the "weightedSampleCoeffs" method. This code assumes that
        // weights and coordinates are statistically independent from
        // each other.
        */
        template <class Pair>
        npstat::Matrix<double>
        weightedSampleCoeffsCovariance(const Pair* points, unsigned long nPoints,
                                       const double* coeffs, unsigned maxdeg) const;
 
        /**
        // Integrated squared error between the given function and the
        // polynomial series
        */
        template <class Functor, class Quadrature>
        double integratedSquaredError(const Functor& fcn,
                                      const Quadrature& quad,
                                      const double* coeffs,
                                      unsigned maxdeg) const;
 
        /**
        // Integrated squared error between the given function and the
        // polynomial series on an arbitrary interval. If the given
        // interval has regions beyond xmin() and xmax(), the integration
        // will assume that on those regions the value of the series is 0,
        // and the quadrature will be carried out there separately.
        */
        template <class Functor, class Quadrature>
        double integratedSquaredError(const Functor& fcn,
                                      const Quadrature& quad,
                                      const double* coeffs,
                                      unsigned maxdeg,
                                      double xmin, double xmax) const;
 
        /**
        // A helper function for choosing the density support interval if
        // it is unknown and has to be estimated from the sample. If the
        // limit is not from the sample (as indicated by the corresponding
        // boolean argument), it is left unchanged. The first element of
        // the returned pair will correspond to the lower limit of the
        // interval and the second element to the upper.
        */
        static std::pair<double,double> supportRegion(
            const AbsClassicalOrthoPoly1D& poly1d,
            double leftLimit, bool leftIsFromSample,
            double rightLimit, bool rightIsFromSample,
            unsigned long nPoints);
 
        /**
        // Optimal OSDE truncation degree according to the Hart scheme.
        // Minimizes Sum_{i=0}^M (k*variances[i] - coeffs[i]^2) over M.
        // The default value of k = 2 corresponds to the original scheme.
        // Argument arrays "coeffs" and "variances" should have at least
        // maxdeg + 1 elements and should be calculated in advance
        // with "sampleCoeffs"/"sampleCoeffsVariance" or with
        // "weightedSampleCoeffs"/"weightedSampleCoeffsVariance" methods.
        */
        static unsigned optimalDegreeHart(const double* coeffs,
                                          const double* variances,
                                          unsigned maxdeg, double k = 2.0);
 
    private:
        const AbsClassicalOrthoPoly1D* poly_;
        double xmin_;
        double xmax_;
        double shift_;
        double scale_;
 
#ifdef SWIG
    public:
        inline double integrateSeries_2(
            const double* coeffs, unsigned nCoeffs) const
        {
            assert(nCoeffs);
            return scale_*poly_->integrateSeries(coeffs, nCoeffs-1U);
        }
 
        inline double series_2(const double* coeffs, unsigned nCoeffs,
                               const double x) const
        {
            assert(nCoeffs);
            return this->series(coeffs, nCoeffs-1U, x);
        }
 
        inline PyObject* seriesScan(const double* coeffs, unsigned nCoeffs,
                                    const unsigned nOut) const
        {
            assert(nCoeffs);
            const unsigned maxdeg = nCoeffs - 1U;
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = nOut;
            PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
            PyObject* yarr = PyArray_SimpleNew(1, &sh, typenum);
            if (xarr && yarr)
            {
                if (nOut)
                {
                    try
                    {
                        double* x = (double*)PyArray_DATA((PyArrayObject*)xarr);
                        double* y = (double*)PyArray_DATA((PyArrayObject*)yarr);
                        const double dx = (xmax_ - xmin_)/nOut;
                        for (unsigned i=0; i<nOut; ++i)
                        {
                            x[i] = xmin_ + (i + 0.5)*dx;
                            y[i] = this->series(coeffs, maxdeg, x[i]);
                        }
                    }
                    catch (const std::exception& e)
                    {
                        Py_DECREF(xarr);
                        Py_DECREF(yarr);
                        throw e;
                    }
                }
                return Py_BuildValue("(OO)", xarr, yarr);
            }
            else
            {
                Py_XDECREF(xarr);
                Py_XDECREF(yarr);
                return 0;
            }
        }
 
        inline void densityCoeffs_2(const AbsDistribution1D& distro,
                                    unsigned nIntegrationPoints,
                                    double* coeffsOut, unsigned nCoeffsOut) const
        {
            assert(nCoeffsOut);
            this->densityCoeffs(distro, nIntegrationPoints, coeffsOut, nCoeffsOut-1);
        }
 
        inline void densityCoeffsVariance_2(const AbsDistribution1D& distro,
                                            unsigned nIntegrationPoints,
                                            const double* coeffs, unsigned nCoeffs,
                                            double expectedSampleSize,
                                            double* coeffsOut, unsigned nCoeffsOut) const
        {
            assert(nCoeffs);
            if (!(nCoeffs == nCoeffsOut)) throw std::invalid_argument(
                "In npstat::OSDE1D::densityCoeffsVariance_2: incompatible output array");
            this->densityCoeffsVariance(distro, nIntegrationPoints, coeffs,
                                        nCoeffs-1U, expectedSampleSize, coeffsOut);
        }
 
        inline npstat::Matrix<double>
        densityCoeffsCovariance_2(const AbsDistribution1D& distro,
                                  unsigned nIntegrationPoints,
                                  const double* coeffs, unsigned nCoeffs,
                                  double expectedSampleSize) const
        {
            assert(nCoeffs);
            return this->densityCoeffsCovariance(distro, nIntegrationPoints,
                                                 coeffs, nCoeffs-1, expectedSampleSize);
        }
 
        template <class Vec>
        inline void sampleCoeffs_2(const Vec& sample,
                                   double* coeffsOut, unsigned nCoeffsOut) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffsOut);
            this->sampleCoeffs(in, sz, coeffsOut, nCoeffsOut-1U);
        }
 
        template <class Vec>
        inline void sampleCoeffsVariance_2(const Vec& sample,
                                           const double *coeffs, unsigned nCoeffs,
                                           double* coeffsOut, unsigned nCoeffsOut) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffs);
            if (!(nCoeffs == nCoeffsOut)) throw std::invalid_argument(
                "In npstat::OSDE1D::sampleCoeffsVariance_2: incompatible output array");
            this->sampleCoeffsVariance(in, sz, coeffs, nCoeffs-1U, coeffsOut);
        }
 
        template <class Vec>
        inline npstat::Matrix<double> sampleCoeffsCovariance_2(
            const Vec& sample, const double *coeffs, unsigned nCoeffs) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffs);
            return this->sampleCoeffsCovariance(in, sz, coeffs, nCoeffs-1U);
        }
 
        template <class Quadrature>
        inline void calculateCoeffs_2(const AbsDistribution1D& d,
                                      const Quadrature& quad,
                                      double* coeffsOut, unsigned nCoeffsOut) const
        {
            DensityFunctor1D fcn(d);
            assert(nCoeffsOut);
            return this->calculateCoeffs(fcn, quad, coeffsOut, nCoeffsOut-1U);
        }
 
        template <class Quadrature>
        inline double integratedSquaredError_2(const AbsDistribution1D& d,
                                               const Quadrature& quad,
                                               const double* coeffs, unsigned nCoeffs) const
        {
            DensityFunctor1D fcn(d);
            assert(nCoeffs);
            return this->integratedSquaredError(fcn, quad, coeffs, nCoeffs-1U);
        }
 
        template <class Vec>
        inline void weightedSampleCoeffs_2(const Vec& sample,
                                           double* coeffsOut, unsigned nCoeffsOut) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffsOut);
            this->weightedSampleCoeffs(in, sz, coeffsOut, nCoeffsOut-1U);
        }
 
        template <class Vec>
        inline npstat::Matrix<double> weightedSampleCoeffsCovariance_2(
            const Vec& sample, const double *coeffs, unsigned nCoeffs) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffs);
            return this->weightedSampleCoeffsCovariance(in, sz, coeffs, nCoeffs-1U);
        }
 
        template <class Vec>
        inline void
        weightedSampleCoeffsVariance_2(const Vec& sample,
                                       const double *coeffs, unsigned nCoeffs,
                                       double* coeffsOut, unsigned nCoeffsOut) const
        {
            const unsigned long sz = sample.size();
            const typename Vec::value_type *in = 0;
            if (sz)
                in = &sample[0];
            assert(nCoeffs);
            if (!(nCoeffs == nCoeffsOut)) throw std::invalid_argument(
                "In npstat::OSDE1D::weightedSampleCoeffsVariance_2: incompatible output array");
            this->weightedSampleCoeffsVariance(in, sz, coeffs, nCoeffs-1U, coeffsOut);
        }
 
        inline static unsigned optimalDegreeHart_2(const double* coeffs, unsigned nCoeffs,
                                                   const double* variances, unsigned nVariances,
                                                   double k = 2.0)
        {
            assert(nCoeffs);
            if (!(nCoeffs == nVariances)) throw std::invalid_argument(
                "In npstat::OSDE1D::optimalDegreeHart_2: incompatible input arrays");
            return optimalDegreeHart(coeffs, variances, nCoeffs-1U, k);
        }
#endif // SWIG
    };
}
 
#include "npstat/stat/OSDE1D.icc"
 
#endif // NPSTAT_OSDE1D_HH_