Index: trunk/npstat/stat/AbsCGF1D.hh =================================================================== --- trunk/npstat/stat/AbsCGF1D.hh (revision 755) +++ trunk/npstat/stat/AbsCGF1D.hh (revision 756) @@ -1,92 +1,92 @@ #ifndef NPSTAT_ABSCGF1D_HH_ #define NPSTAT_ABSCGF1D_HH_ /*! // \file AbsCGF1D.hh // // \brief Interface definition for univariate cumulant generating functions // // Author: I. Volobouev // // December 2019 */ #include namespace npstat { struct AbsCGF1D { inline virtual ~AbsCGF1D() {} /** Virtual copy constructor */ virtual AbsCGF1D* clone() const = 0; /** // CGF for the variable (X - mu)/sigma. Created on the // heap and must be later deleted by the user. */ virtual AbsCGF1D* shiftAndScale(double mu, double sigma) const = 0; /** Infimum of the support */ virtual double smin() const = 0; /** Supremum of the support */ virtual double smax() const = 0; /** CGF value */ inline double operator()(double s) const {return derivative(0U, s);} - /** CGF derivatives */ + /** CGF derivatives. Get cum[order] by setting s to 0. */ virtual double derivative(unsigned order, double s) const = 0; /** Solution of the saddlepoint equation CGF'(s) == x */ virtual double saddlepoint(double x) const; /** // Search empirically for the maximum value of s // for which this CGF remains convex. This function // should search for the closest to 0 root of the // CGF second derivative on the interval [0, slimit]. // It should return 0.0 if no such root was found. */ virtual double convexLimit(double slimit, double stepSize) const; /** Standardized cumulant for s-tilted density */ double standardizedCumulant(unsigned order, double s) const; /** // Saddlepoint density approximation. Parameter "order" is the // approximation order (e.g., order = 2 means O(N^-2)). This // parameter must be positive. The raw, unnormalized approximation // is returned. */ double saddlepointDensityApprox(unsigned order, double x) const; /** // Lugannani-Rice approximation for the survival function. // The approximation employs a series which looks like // c0/N^{1/2} + c1/N^{3/2} + c2/N^{5/2} + ... // The parameter "nTerms" specifies the number of terms in // the series to use and, thereby, the approximation order. */ double lugannaniRiceSFApprox(unsigned nTerms, double x) const; /** // Derived classes should not implement "operator==", implement // "isEqual" instead */ inline bool operator==(const AbsCGF1D& r) const {return (typeid(*this) == typeid(r)) && this->isEqual(r);} /** Logical negation of operator== */ inline bool operator!=(const AbsCGF1D& r) const {return !(*this == r);} protected: /** Comparison for equality. To be implemented by derived classes. */ virtual bool isEqual(const AbsCGF1D&) const = 0; }; } #endif // NPSTAT_ABSCGF1D_HH_ Index: trunk/npstat/stat/SeriesCGF1D.cc =================================================================== --- trunk/npstat/stat/SeriesCGF1D.cc (revision 755) +++ trunk/npstat/stat/SeriesCGF1D.cc (revision 756) @@ -1,108 +1,233 @@ #include #include #include +#include +#include +#include #include "npstat/stat/SeriesCGF1D.hh" #include "npstat/rng/permutation.hh" #include "npstat/nm/Poly1D.hh" +#include "npstat/nm/findRootUsingBisections.hh" + +namespace { + inline double guessExtraCumulant(const npstat::SeriesCGF1D& cgf, + const double a, const unsigned newOrder) + { + assert(newOrder % 2U == 0U); + assert(newOrder > cgf.maxOrder()); + const double sderi = cgf.derivative(2U, a); + if (sderi < 0.0) + { + const unsigned nm2 = newOrder-2U; + return -sderi * npstat::ldfactorial(nm2)/powl(a, nm2); + } + else + return 0.0; + } + + class ConvexificationFcn : public npstat::Functor1 + { + public: + inline ConvexificationFcn(npstat::SeriesCGF1D& cgf, + const double a, const double b, + const unsigned tune) + : cgf_(cgf), a_(a), b_(b), tune_(tune) {} + + inline double operator()(const double& cumGuess) const + { + cgf_.setCumulant(tune_, cumGuess); + if (cgf_.isConvex(a_, b_)) + return 1.0; + else + return -1.0; + } + + private: + npstat::SeriesCGF1D& cgf_; + double a_; + double b_; + unsigned tune_; + }; +} namespace npstat { SeriesCGF1D::SeriesCGF1D(const double* cumulants, const unsigned nCumulants, const double i_smin, const double i_smax) : cumulants_(cumulants, cumulants+nCumulants), smin_(i_smin), smax_(i_smax) { assert(cumulants); if (nCumulants < 2U) throw std::invalid_argument( "In npstat::SeriesCGF1D constructor: not enough cumulants"); if (cumulants[1] <= 0.0) throw std::invalid_argument( "In npstat::SeriesCGF1D constructor: variance must be positive"); if (!(smin_ < 0.0)) throw std::invalid_argument( "In npstat::SeriesCGF1D constructor: invalid support infimum"); if (!(smax_ > 0.0)) throw std::invalid_argument( "In npstat::SeriesCGF1D constructor: invalid support supremum"); + truncateLeadingZeros(); + } + + void SeriesCGF1D::setCumulant(const unsigned order, const double value) + { + if (order) + { + const unsigned maxO = cumulants_.size(); + if (order <= maxO) + { + if (order == 2U && value <= 0.0) throw std::invalid_argument( + "In npstat::SeriesCGF1D setCumulant: variance must be positive"); + cumulants_[order-1U] = value; + if (order == maxO && value == 0.0) + truncateLeadingZeros(); + } + else if (value) + { + cumulants_.resize(order); + cumulants_[order-1U] = value; + } + } + else if (value) + { + throw std::invalid_argument( + "In npstat::SeriesCGF1D setCumulant: can't change 0th order cumulant"); + } + } + + void SeriesCGF1D::truncateLeadingZeros() + { + const unsigned sz = cumulants_.size(); + unsigned firstNon0 = sz; + for (unsigned i=firstNon0-1U; i>0U; --i) + if (cumulants_[i]) + { + firstNon0 = i; + break; + } + assert(firstNon0 < sz); + if (firstNon0 < sz - 1U) + cumulants_.resize(firstNon0 + 1U); } bool SeriesCGF1D::isEqual(const AbsCGF1D& other) const { const SeriesCGF1D& r = static_cast(other); return cumulants_ == r.cumulants_ && smin_ == r.smin_ && smax_ == r.smax_; } double SeriesCGF1D::derivative(const unsigned order, const double s) const { if (s == 0.0) { if (order) { const unsigned idx = order - 1U; if (idx < cumulants_.size()) return cumulants_[idx]; } } else if (!(smin_ < s && s < smax_)) throw std::invalid_argument( "In npstat::SeriesCGF1D::derivative: argument out of range"); else { const unsigned degree = cumulants_.size(); if (order <= degree) { const double* a = &cumulants_[0]; long double sum = 0.0L; if (order) { const int shift = order-1U; for (int deg=degree-order; deg>=0; --deg) { sum *= s; sum += a[deg+shift]/ldfactorial(deg); } } else { for (unsigned deg=degree; deg>0U; --deg) { sum += a[deg-1]/ldfactorial(deg); sum *= s; } } return sum; } } return 0.0; } SeriesCGF1D* SeriesCGF1D::shiftAndScale(const double mu, const double sigma) const { if (!(sigma > 0.0)) throw std::invalid_argument( "In npstat::SeriesCGF1D::shiftAndScale: scale must be positive"); std::vector cum(cumulants_); const unsigned sz = cum.size(); cum[0] -= mu; if (sigma != 1.0) for (unsigned i=0; i= b) throw std::invalid_argument( "In npstat::SeriesCGF1D::sisConvex: invalid interval"); const unsigned degree = cumulants_.size(); std::vector coeffs(degree + 1U); coeffs[0] = 0.0L; for (unsigned deg=0; deg 0.0L && !cgf2.nRoots(a, b); } + + SeriesCGF1D SeriesCGF1D::convexify(const double a, const double b) const + { + if (isConvex(a, b)) + return *this; + + const unsigned maxO = cumulants_.size(); + const unsigned tune = maxO % 2U ? maxO + 1U : maxO + 2U; + const double cumGuessA = guessExtraCumulant(*this, a, tune); + const double cumGuessB = guessExtraCumulant(*this, b, tune); + double cumGuess = std::max(cumGuessA, cumGuessB); + SeriesCGF1D copy(*this); + if (cumGuess > 0.0) + { + copy.setCumulant(tune, cumGuess); + const double step = (b - a)*sqrt(DBL_EPSILON); + if (copy.isConvex(a + step, b - step)) + return copy; + } + + cumGuess = std::max(cumGuess, 1.0); + copy.setCumulant(tune, cumGuess); + while (!copy.isConvex(a, b)) + { + assert(cumGuess < DBL_MAX/2.0); + cumGuess *= 2.0; + copy.setCumulant(tune, cumGuess); + } + + ConvexificationFcn fcn(copy, a, b, tune); + double cum = 0.0; + const double tol = cumGuess*DBL_EPSILON*2.0; + const bool status = findRootUsingBisections( + fcn, 0.0, 0.0, cumGuess, tol, &cum); + assert(status); + copy.setCumulant(tune, cum); + return copy; + } } Index: trunk/npstat/stat/SeriesCGF1D.hh =================================================================== --- trunk/npstat/stat/SeriesCGF1D.hh (revision 755) +++ trunk/npstat/stat/SeriesCGF1D.hh (revision 756) @@ -1,65 +1,78 @@ #ifndef NPSTAT_SERIESCGF1D_HH_ #define NPSTAT_SERIESCGF1D_HH_ /*! // \file SeriesCGF1D.hh // // \brief Univariate cumulant generating function defined by cumulants // // Author: I. Volobouev // // December 2019 */ #include #include "npstat/stat/AbsCGF1D.hh" namespace npstat { class SeriesCGF1D : public AbsCGF1D { public: /** // At least the first two cumulants (mean and variance) // must be provided in the "cumulants" array. // cumulants[0] must contain the mean, cumulants[1] the // variance, etc. "nCumulants" specifies the length of // the array (and also the highest cumulant order). // The variance must be positive (i.e., the degenerate // distribution of a single point mass is excluded from // consideration). The parameters "smin" and "smax" // should be chosen in such a way that the CGF remains // convex on the (smin, smax) interval. */ SeriesCGF1D(const double* cumulants, unsigned nCumulants, double smin, double smax); inline virtual SeriesCGF1D* clone() const {return new SeriesCGF1D(*this);} inline virtual ~SeriesCGF1D() {} virtual SeriesCGF1D* shiftAndScale(double mu, double sigma) const; inline virtual double smin() const {return smin_;} inline virtual double smax() const {return smax_;} virtual double derivative(unsigned order, double s) const; + void setCumulant(unsigned order, double value); + + inline unsigned maxOrder() const {return cumulants_.size();} + /** // Check whether the CGF is convex on the interval (a, b]. // This method ignores smin and smax. */ - virtual bool isConvex(double a, double b) const; + bool isConvex(double a, double b) const; + + /** + // Convexify this CFG on the given interval by adding the + // smallest possible higher order cumulant of the lowest + // even order not specified for this CGF. + */ + SeriesCGF1D convexify(double a, double b) const; protected: virtual bool isEqual(const AbsCGF1D& r) const; private: + void truncateLeadingZeros(); + std::vector cumulants_; double smin_; double smax_; }; } #endif // NPSTAT_SERIESCGF1D_HH_