Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh =================================================================== --- trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh (revision 700) +++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh (revision 701) @@ -1,411 +1,414 @@ #ifndef NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ #define NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ /*! // \file AbsClassicalOrthoPoly1D.hh // // \brief Base class for classical continuous orthonormal polynomials // // Author: I. Volobouev // // May 2017 */ #include #include #include #include #include "geners/CPP11_auto_ptr.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" namespace npstat { class StorablePolySeries1D; class AbsClassicalOrthoPoly1D { public: inline virtual ~AbsClassicalOrthoPoly1D() {} /** Virtual copy constructor */ virtual AbsClassicalOrthoPoly1D* clone() const = 0; /** The weight function should be normalized */ virtual long double weight(long double x) const = 0; //@{ /** Support of the polynomial */ virtual double xmin() const = 0; virtual double xmax() const = 0; //@} /** Maximum polynomial degree supported */ inline virtual unsigned maxDegree() const {return UINT_MAX - 1U;} /** Polynomial values */ long double poly(unsigned deg, long double x) const; /** // Values of two orthonormal polynomials. // Faster than calling "poly" two times. */ std::pair twopoly( unsigned deg1, unsigned deg2, long double x) const; /** // Values of all orthonormal polynomials up to some degree. // Faster than calling "poly" multiple times. The size of // the "values" array should be at least maxdeg + 1. */ void allpoly(long double x, long double* values, unsigned maxdeg) const; /** Polynomial series */ double series(const double* coeffs, unsigned maxdeg, double x) const; + /** Integral of the series */ + double integrateSeries(const double* coeffs, 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 void calculateCoeffs(const Functor& fcn, const Quadrature& quad, double* coeffs, unsigned maxdeg) 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 void sampleCoeffs(const Numeric* coords, unsigned long lenCoords, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs". The "variances" array // should have at least maxdeg + 1 elements. It will contain // the respective variances upon return. */ template void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords, const double* coeffs, unsigned maxdeg, double* variances) 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). // Before calculating the coefficients, the coordinates are shifted // and scaled according to x_new = (x_original - location)/scale. // The resulting coefficients are also divided by scale. */ template void sampleCoeffs(const Numeric* coords, unsigned long lenCoords, Numeric location, Numeric scale, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs" with the same location // and scale. The "variances" array should have at least maxdeg + 1 // elements. It will contain the respective variances upon return. */ template void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords, Numeric location, Numeric scale, const double* coeffs, unsigned maxdeg, double* variances) 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 void weightedSampleCoeffs(const Pair* points, unsigned long numPoints, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs". The "variances" // array should have at least maxdeg + 1 elements. It will contain // the respective variances upon return. This code assumes that // weights and coordinates are statistically independent from // each other. */ template void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints, const double* coeffs, unsigned maxdeg, double* variances) 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). Before calculating the coefficients, the // coordinates are shifted and scaled according to // x_new = (x_original - location)/scale. The resulting coefficients // are also divided by scale. */ template void weightedSampleCoeffs(const Pair* points, unsigned long numPoints, typename Pair::first_type location, typename Pair::first_type scale, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs" with the same // location and scale. The "variances" array should have at least // maxdeg + 1 elements. It will contain the respective variances // upon return. This code assumes that weights and coordinates // are statistically independent from each other. */ template void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints, typename Pair::first_type location, typename Pair::first_type scale, const double* coeffs, unsigned maxdeg, double* variances) const; /** // This method is useful for testing the numerical precision // of the orthonormalization procedure. It returns the scalar // products between various polynomials. */ template double empiricalKroneckerDelta(const Quadrature& quad, unsigned deg1, unsigned deg2) const; /** // A measure-weighted average of a product of four orthonormal // poly values */ template double jointAverage(const Quadrature& q, unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; /** // A measure-weighted average of a product of multiple orthonormal // poly values. "checkedForZeros" argument should be set to "true" // if we are sure that there are no zero elements in the "degrees" // array. */ template double jointAverage(const Quadrature& q, const unsigned* degrees, unsigned nDegrees, bool checkedForZeros = false) const; /** // An unweighted integral of the orthonormal polynomial with the // given degree to some power. For this method, it is assumed // that the polynomials are supported on a closed interval (without // such an assumption unweighted integrals do not make much sense) // and that Gauss-Legendre quadratures can be used. */ double integratePoly(unsigned degree, unsigned power) const; /** // An unweighted integral of a product of multiple orthonormal // polynomials. For this method, it is assumed that the polynomials // are supported on a closed interval (without such an assumption // unweighted integrals do not make much sense) and that Gauss-Legendre // quadratures can be used. */ double jointIntegral(const unsigned* degrees, unsigned nDegrees) const; /** // Unweighted averages of the polynomial values for the given sample. // The length of array "averages" should be at least maxdeg + 1. */ template void sampleAverages(const Numeric* coords, unsigned long lenCoords, double* averages, unsigned maxdeg) const; /** // Unweighted averages of the pairwise products of the polynomial // values for the given sample. The returned matrix will be symmetric // and will have the dimensions (maxdeg + 1) x (maxdeg + 1). */ template Matrix sampleProductAverages(const Numeric* coords, unsigned long lenCoords, unsigned maxdeg) const; /** // Unweighted average of a product of polynomial values for the // given sample. "degrees" is the collection of polynomial degrees. // Polynomials of these degrees will be included in the product. */ template double jointSampleAverage(const Numeric* coords, unsigned long lenCoords, const unsigned* degrees, unsigned lenDegrees) const; // Similar to the previous method but calculates two averages // simultaneously template std::pair twoJointSampleAverages( const Numeric* coords, unsigned long lenCoords, const unsigned* degrees1, unsigned lenDegrees1, const unsigned* degrees2, unsigned lenDegrees2) const; CPP11_auto_ptr makeStorablePolySeries( unsigned maxPolyDeg, const double *coeffs=0, unsigned maxdeg=0) const; protected: inline static int kdelta(const unsigned i, const unsigned j) {return i == j ? 1 : 0;} // Recurrence relationship function should return alpha // as the first element of the pair and the square root // of beta as the second virtual std::pair recurrenceCoeffs(unsigned deg) const = 0; private: // Helper classes class ProdFcn : public Functor1 { public: inline ProdFcn(const AbsClassicalOrthoPoly1D& p, const unsigned d1, const unsigned d2) : poly(p), deg1(d1), deg2(d2) {} inline long double operator()(const long double& x) const { const std::pair& p = poly.twopoly(deg1, deg2, x); return p.first*p.second*poly.weight(x); } private: const AbsClassicalOrthoPoly1D& poly; unsigned deg1; unsigned deg2; }; class MultiProdFcn; friend class MultiProdFcn; class MultiProdFcn : public Functor1 { public: inline MultiProdFcn(const AbsClassicalOrthoPoly1D& p, const unsigned* degrees, unsigned lenDegrees, const bool includeWeight=true) : poly(p), degs(degrees, degrees+lenDegrees), maxdeg(0), weighted(includeWeight) { if (lenDegrees) maxdeg = *std::max_element(degrees, degrees+lenDegrees); } inline long double operator()(const long double& x) const { long double w = 1.0L, p = 1.0L; if (weighted) w = poly.weight(x); if (maxdeg) p = poly.normpolyprod(°s[0], degs.size(), maxdeg, x); return w*p; } private: const AbsClassicalOrthoPoly1D& poly; std::vector degs; unsigned maxdeg; bool weighted; }; // For calling the function below, maxdeg should be calculated as // *std::max_element(degrees, degrees+nDegrees) long double normpolyprod(const unsigned* degrees, unsigned nDegrees, unsigned maxdeg, long double x) const; #ifdef SWIG public: inline double weight2(const double x) const {return weight(x);} inline double poly2(const unsigned deg, const double x) const {return poly(deg, x);} inline std::vector allpoly2(const unsigned maxdeg, const double x) const { std::vector p(maxdeg+1); allpoly(x, &p[0], p.size()); return std::vector(p.begin(), p.end()); } inline StorablePolySeries1D* makeStorablePolySeries_2( unsigned maxPolyDeg, const std::vector& coeffs) const { CPP11_auto_ptr ptr; if (coeffs.empty()) ptr = this->makeStorablePolySeries(maxPolyDeg); else ptr = this->makeStorablePolySeries(maxPolyDeg, &coeffs[0], coeffs.size() - 1); return ptr.release(); } #endif // SWIG }; /** // A functor for the weight function of the given ortho poly system. // The poly system is not copied, only a reference is used. It is // a responsibility of the user to make sure that the lifetime of // the poly system object exceeds the lifetime of the functor. */ class OrthoPoly1DWeight : public Functor1 { public: inline explicit OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn, const long double normfactor=1.0L) : fcn_(fcn), norm_(normfactor) {} inline virtual ~OrthoPoly1DWeight() {} inline virtual long double operator()(const long double& a) const {return norm_*fcn_.weight(a);} private: OrthoPoly1DWeight(); const AbsClassicalOrthoPoly1D& fcn_; const long double norm_; }; } #include "npstat/nm/AbsClassicalOrthoPoly1D.icc" #endif // NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.cc =================================================================== --- trunk/npstat/nm/AbsClassicalOrthoPoly1D.cc (revision 700) +++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.cc (revision 701) @@ -1,245 +1,281 @@ #include #include #include #include #include "npstat/nm/StorablePolySeries1D.hh" #include "npstat/nm/GaussLegendreQuadrature.hh" namespace npstat { namespace Private { + class ClassicalPolySeries : public Functor1 + { + public: + inline ClassicalPolySeries(const AbsClassicalOrthoPoly1D& p, + const double *coeffs, const unsigned maxdeg) + : poly_(p), coeffs_(coeffs), maxdeg_(maxdeg) {} + + inline double operator()(const double& x) const + { + return poly_.series(coeffs_, maxdeg_, x); + } + + private: + const AbsClassicalOrthoPoly1D& poly_; + const double* coeffs_; + unsigned maxdeg_; + }; + class ClassicalPolyPower : public Functor1 { public: inline ClassicalPolyPower(const AbsClassicalOrthoPoly1D& p, const unsigned deg, const unsigned power, const bool i_useWeight) : poly(p), deg1(deg), polypow(power), useWeight(i_useWeight) {} inline long double operator()(const long double& x) const { long double p = 1.0L, w = 1.0L; if (useWeight) w = poly.weight(x); if (polypow) { p = poly.poly(deg1, x); switch (polypow) { case 1U: break; case 2U: p *= p; break; default: p = powl(p, polypow); break; } } return w*p; } private: const AbsClassicalOrthoPoly1D& poly; unsigned deg1; unsigned polypow; bool useWeight; }; } long double AbsClassicalOrthoPoly1D::poly( const unsigned degree, const long double x) const { long double polyk = 1.0L; if (degree) { if (degree > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::poly: degree argument is out of range"); long double polykm1 = 0.0L; std::pair rcurrent = recurrenceCoeffs(0); for (unsigned k=0; k& rnext = recurrenceCoeffs(k+1); const long double p = ((x - rcurrent.first)*polyk - rcurrent.second*polykm1)/rnext.second; polykm1 = polyk; polyk = p; rcurrent = rnext; } } return polyk; } long double AbsClassicalOrthoPoly1D::normpolyprod( const unsigned* degrees, const unsigned nDegrees, const unsigned maxdeg, const long double x) const { long double prod = 1.0L; if (nDegrees) { assert(degrees); if (maxdeg) { if (maxdeg > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::normpolyprod: " "maximum degree argument is out of range"); long double polyk = 1.0L, polykm1 = 0.0L; std::pair rcurrent = recurrenceCoeffs(0); for (unsigned k=0; k& rnext = recurrenceCoeffs(k+1); const long double p = ((x - rcurrent.first)*polyk - rcurrent.second*polykm1)/rnext.second; polykm1 = polyk; polyk = p; rcurrent = rnext; } for (unsigned j=0; j AbsClassicalOrthoPoly1D::twopoly( const unsigned deg1, const unsigned deg2, const long double x) const { long double p1 = 1.0L, p2 = 1.0L; const unsigned degree = std::max(deg1, deg2); if (degree) { if (degree > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::twopoly: " "degree argument is out of range"); long double polyk = 1.0L, polykm1 = 0.0L; std::pair rcurrent = recurrenceCoeffs(0); for (unsigned k=0; k& rnext = recurrenceCoeffs(k+1); const long double p = ((x - rcurrent.first)*polyk - rcurrent.second*polykm1)/rnext.second; polykm1 = polyk; polyk = p; rcurrent = rnext; } if (deg1 == degree) p1 = polyk; if (deg2 == degree) p2 = polyk; } return std::pair(p1, p2); } double AbsClassicalOrthoPoly1D::series(const double *coeffs, const unsigned degree, const double xIn) const { assert(coeffs); long double sum = coeffs[0]; if (degree) { if (degree > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::series: " "degree argument is out of range"); const long double x = xIn; long double polyk = 1.0L, polykm1 = 0.0L; std::pair rcurrent = recurrenceCoeffs(0); for (unsigned k=0; k& rnext = recurrenceCoeffs(k+1); const long double p = ((x - rcurrent.first)*polyk - rcurrent.second*polykm1)/rnext.second; sum += p*coeffs[k+1]; polykm1 = polyk; polyk = p; rcurrent = rnext; } } return sum; } + double AbsClassicalOrthoPoly1D::integrateSeries(const double *coeffs, + const unsigned degree) const + { + assert(coeffs); + if (degree) + { + const unsigned nPt = GaussLegendreQuadrature::minimalExactRule(degree); + if (!nPt) throw std::invalid_argument( + "In npstat::AbsClassicalOrthoPoly1D::integrateSeries: " + "the number of coefficients is too large"); + GaussLegendreQuadrature quad(nPt); + Private::ClassicalPolySeries fcn(*this, coeffs, degree); + return quad.integrate(fcn, xmin(), xmax()); + } + else + return coeffs[0]*(xmax() - xmin()); + } + void AbsClassicalOrthoPoly1D::allpoly(const long double x, long double* values, const unsigned degree) const { assert(values); values[0] = 1.0L; if (degree) { if (degree > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::allpoly: " "degree argument is out of range"); long double polyk = 1.0L, polykm1 = 0.0L; std::pair rcurrent = recurrenceCoeffs(0); for (unsigned k=0; k& rnext = recurrenceCoeffs(k+1); const long double p = ((x - rcurrent.first)*polyk - rcurrent.second*polykm1)/rnext.second; polykm1 = polyk; polyk = p; rcurrent = rnext; values[k+1] = p; } } } double AbsClassicalOrthoPoly1D::jointIntegral( const unsigned* degrees, const unsigned nDegrees) const { if (nDegrees) { assert(degrees); const unsigned degSum = std::accumulate(degrees, degrees+nDegrees, 0U); const unsigned nPt = GaussLegendreQuadrature::minimalExactRule(degSum); if (!nPt) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::jointIntegral: " "joint poly degree is too high"); GaussLegendreQuadrature quad(nPt); MultiProdFcn fcn(*this, degrees, nDegrees, false); return quad.integrate(fcn, xmin(), xmax()); } else return xmax() - xmin(); } double AbsClassicalOrthoPoly1D::integratePoly( const unsigned degree, const unsigned power) const { const unsigned degSum = degree*power; if (degSum) { const unsigned nPt = GaussLegendreQuadrature::minimalExactRule(degSum); if (!nPt) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::integratePoly: " "poly degree or power is too high"); GaussLegendreQuadrature quad(nPt); Private::ClassicalPolyPower fcn(*this, degree, power, false); return quad.integrate(fcn, xmin(), xmax()); } else return xmax() - xmin(); } CPP11_auto_ptr AbsClassicalOrthoPoly1D::makeStorablePolySeries( const unsigned maxPolyDeg, const double *coeffs, unsigned maxdeg) const { if (maxPolyDeg > maxDegree()) throw std::invalid_argument( "In npstat::AbsClassicalOrthoPoly1D::makeStorablePolySeries: " "poly degree argument is out of range"); std::vector > rc(maxPolyDeg+1U); for (unsigned i=0; i<=maxPolyDeg; ++i) rc[i] = this->recurrenceCoeffs(i); return CPP11_auto_ptr( new StorablePolySeries1D(rc, this->xmin(), this->xmax(), 0.0, coeffs, maxdeg)); } } Index: trunk/npstat/stat/OSDE1D.cc =================================================================== --- trunk/npstat/stat/OSDE1D.cc (revision 700) +++ trunk/npstat/stat/OSDE1D.cc (revision 701) @@ -1,124 +1,130 @@ #include #include #include #include #include "npstat/nm/LinearMapper1d.hh" #include "npstat/stat/OSDE1D.hh" namespace npstat { OSDE1D::OSDE1D(const AbsClassicalOrthoPoly1D& poly1d, const double xlo, const double xhi) : poly_(poly1d.clone()), xmin_(xlo), xmax_(xhi) { if (xmin_ == xmax_) throw std::invalid_argument( "In npstat::OSDE1D constructor: " "density estimation interval has zero length"); if (xmin_ > xmax_) std::swap(xmin_, xmax_); const double polyWidth = poly_->xmax() - poly_->xmin(); scale_ = (xmax_ - xmin_)/polyWidth; shift_ = (poly_->xmax()*xmin_ - poly_->xmin()*xmax_)/polyWidth; } OSDE1D::OSDE1D(const OSDE1D& r) : poly_(r.poly_->clone()), xmin_(r.xmin_), xmax_(r.xmax_), shift_(r.shift_), scale_(r.scale_) { } OSDE1D& OSDE1D::operator=(const OSDE1D& r) { if (&r != this) { delete poly_; poly_ = 0; poly_ = r.poly_->clone(); xmin_ = r.xmin_; xmax_ = r.xmax_; shift_ = r.shift_; scale_ = r.scale_; } return *this; } OSDE1D::~OSDE1D() { delete poly_; } double OSDE1D::series(const double *coeffs, const unsigned maxdeg, const double x) const { return poly_->series(coeffs, maxdeg, (x - shift_)/scale_); } + double OSDE1D::integrateSeries(const double *coeffs, + const unsigned maxdeg) const + { + return scale_*poly_->integrateSeries(coeffs, maxdeg); + } + double OSDE1D::weight(const double x) const { return poly_->weight((x - shift_)/scale_)/scale_; } unsigned OSDE1D::optimalDegreeHart(const double* coeffs, const double* variances, const unsigned maxdeg, const double k) { assert(coeffs); assert(variances); long double sum = 0.0L, bestSum = DBL_MAX; unsigned bestM = 0; for (unsigned i=1U; i<=maxdeg; ++i) { sum += k*variances[i] - coeffs[i]*coeffs[i]; if (sum < bestSum) { bestSum = sum; bestM = i; } } return bestM; } std::pair OSDE1D::supportRegion( const AbsClassicalOrthoPoly1D& poly1d, const double leftLimit, const bool leftIsFromSample, const double rightLimit, const bool rightIsFromSample, const unsigned long nPoints) { if (leftLimit >= rightLimit) throw std::invalid_argument( "In npstat::OSDE1D::supportRegion: " "the left limit value must be less than the right limit value"); if (!(leftIsFromSample || rightIsFromSample)) return std::pair(leftLimit, rightLimit); const unsigned nMin = leftIsFromSample && rightIsFromSample ? 2U : 1U; if (nPoints < nMin) throw std::invalid_argument( "In npstat::OSDE1D::supportRegion: insufficient sample size"); const double xmin = poly1d.xmin(); const double xmax = poly1d.xmax(); const double halfstep = (xmax - xmin)/2.0/nPoints; if (leftIsFromSample && rightIsFromSample) { LinearMapper1d mp(xmin + halfstep, leftLimit, xmax - halfstep, rightLimit); return std::pair(mp(xmin), mp(xmax)); } else if (leftIsFromSample) { LinearMapper1d mp(xmin + halfstep, leftLimit, xmax, rightLimit); return std::pair(mp(xmin), rightLimit); } else { LinearMapper1d mp(xmin, leftLimit, xmax - halfstep, rightLimit); return std::pair(leftLimit, mp(xmax)); } } } Index: trunk/npstat/stat/OSDE1D.hh =================================================================== --- trunk/npstat/stat/OSDE1D.hh (revision 700) +++ trunk/npstat/stat/OSDE1D.hh (revision 701) @@ -1,317 +1,320 @@ #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 #include #include "npstat/nm/AbsClassicalOrthoPoly1D.hh" #include "npstat/stat/AbsDistribution1D.hh" #ifdef SWIG #include #include #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 */ 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_);} /** Polynomial series representing the fitted density */ double series(const double* coeffs, unsigned maxdeg, double x) const; + /** Integral of the series */ + double integrateSeries(const double* coeffs, 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 void calculateCoeffs(const Functor& fcn, const Quadrature& quad, double* coeffs, unsigned maxdeg) 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 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 void sampleCoeffsVariance(const Numeric* coords, unsigned long lenCoords, const double* coeffs, unsigned maxdeg, double* variances) 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 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 void weightedSampleCoeffsVariance(const Pair* points, unsigned long nPoints, const double* coeffs, unsigned maxdeg, double* variances) const; /** // Integrated squared error between the given function and the // polynomial series */ template 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 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 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"/"sampleCoeffVars" or with // "weightedSampleCoeffs"/"weightedSampleCoeffVars" 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 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::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; iseries(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; } } template 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 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 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 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 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 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_