Index: trunk/npstat/stat/AbsCGF1D.cc =================================================================== --- trunk/npstat/stat/AbsCGF1D.cc (revision 710) +++ trunk/npstat/stat/AbsCGF1D.cc (revision 711) @@ -1,224 +1,263 @@ #include #include #include #include +#include #include "npstat/nm/findRootUsingBisections.hh" #include "npstat/stat/AbsCGF1D.hh" #include "npstat/stat/Distributions1D.hh" namespace { class AbsCGF1DDeriv : public npstat::Functor1 { public: - inline AbsCGF1DDeriv(const npstat::AbsCGF1D& cgf) - : cgf_(cgf) {} + inline AbsCGF1DDeriv(const npstat::AbsCGF1D& cgf, + const unsigned degree) + : cgf_(cgf), deg_(degree) {} + inline double operator()(const double& s) const - {return cgf_.derivative(1U, s);} + {return cgf_.derivative(deg_, s);} private: const npstat::AbsCGF1D& cgf_; + unsigned deg_; }; } namespace npstat { + double AbsCGF1D::convexLimit(const double slimit, + const double stepSize) const + { + if (derivative(2U, 0.0) <= 0.0) throw std::runtime_error( + "In npstat::AbsCGF1D::convexLimit: " + "the CGF second derivative is not positive at s = 0"); + if (slimit == 0.0) + return 0.0; + if (stepSize == 0.0) throw std::invalid_argument( + "In npstat::AbsCGF1D::convexLimit: " + "step size must not be 0"); + const double step = slimit > 0.0 ? std::abs(stepSize) : -std::abs(stepSize); + const double slo = smin(); + const double shi = smax(); + const double lim = std::abs(slimit); + const double eps = std::numeric_limits::epsilon(); + double olds = 0.0; + for (double s=step; std::abs(s) <= lim && slo < s && s < shi; s += step) + { + if (s == olds) throw std::invalid_argument( + "In npstat::AbsCGF1D::convexLimit: the step size is too small"); + if (derivative(2U, s) <= 0.0) + { + double root; + AbsCGF1DDeriv deriv(*this, 2U); + const bool status = findRootUsingBisections(deriv, 0.0, olds, s, + 2.0*eps, &root); + assert(status); + return root; + } + olds = s; + } + return 0.0; + } + double AbsCGF1D::saddlepoint(const double x) const { static const double factor = sqrt(2.0); const double infimum = smin(); assert(infimum < 0.0); const double supremum = smax(); assert(supremum > 0.0); // The first derivative is monotonously increasing const double d0 = derivative(1U, 0.0); if (x == d0) return 0.0; double xmin = 0.0, xmax = 0.0; double dmin = d0, dmax = d0; const double eps = std::numeric_limits::epsilon(); if (x > d0) { // Need to find useful xmax xmax = 1.0; while (xmax < supremum) { dmax = derivative(1U, xmax); if (x <= dmax) break; if (xmax >= std::numeric_limits::max()/1.001/factor) break; xmax *= factor; } if (x > dmax) { xmax = supremum/(1.0 + eps); dmax = derivative(1U, xmax); } if (!(x <= dmax)) throw std::invalid_argument( "Failed to solve the saddlepoint equation " "within the CGF support"); } else { // Need to find useful xmin xmin = -1.0; while (xmin > infimum) { dmin = derivative(1U, xmin); if (x >= dmin) break; if (xmin < -std::numeric_limits::max()/1.001/factor) break; xmin *= factor; } if (x < dmin) { xmin = infimum/(1.0 + eps); dmin = derivative(1U, xmin); } if (!(x >= dmin)) throw std::invalid_argument( "Failed to solve the saddlepoint equation " "within the CGF support"); } assert(x >= dmin); assert(x <= dmax); if (x == dmin) return xmin; if (x == dmax) return xmax; double root; - AbsCGF1DDeriv deriv(*this); + AbsCGF1DDeriv deriv(*this, 1U); const bool status = findRootUsingBisections(deriv, x, xmin, xmax, 2.0*eps, &root); assert(status); return root; } double AbsCGF1D::standardizedCumulant(const unsigned order, const double s) const { if (order) { if (!(smin() < s && s < smax())) throw std::invalid_argument( "In npstat::AbsCGF1D::standardizedCumulant: " "argument outside of CGF support"); const double var = derivative(2U, s); assert(var > 0.0); return derivative(order, s)/pow(var, order/2.0); } return 0.0; } double AbsCGF1D::saddlepointDensityApprox( const unsigned order, const double x) const { if (!order) throw std::invalid_argument( "In npstat::AbsCGF1D::saddlepointDensityApproximation: " "approximation order must be positive"); if (order > 4U) throw std::invalid_argument( "In npstat::AbsCGF1D::saddlepointDensityApproximation: " "unsupported approximation order (too high)"); const double shat = saddlepoint(x); const double Khat = derivative(0U, shat); const double varhat = derivative(2U, shat); assert(varhat > 0.0); long double factor = 1.0L; if (order > 1U) { const double k3 = standardizedCumulant(3U, shat); const double k4 = standardizedCumulant(4U, shat); const double k3_2 = k3*k3; factor += (k4/8 - 5*k3_2/24); if (order > 2U) { const double k5 = standardizedCumulant(5U, shat); const double k6 = standardizedCumulant(6U, shat); const double k3_4 = k3_2*k3_2; const double k4_2 = k4*k4; factor += ((385*k3_4)/1152 - (35*k3_2*k4)/64 + (35*k4_2)/384 + (7*k3*k5)/48 - k6/48); if (order > 3U) { const double k7 = standardizedCumulant(7U, shat); const double k8 = standardizedCumulant(8U, shat); const double k3_3 = k3*k3_2; const double k3_6 = k3_3*k3_3; const double k4_3 = k4*k4_2; const double k5_2 = k5*k5; factor += ((-85085*k3_6)/82944 + (25025*k3_4*k4)/9216 - (5005*k3_2*k4_2)/3072 + (385*k4_3)/3072 - (1001*k3_3*k5)/1152 + (77*k3*k4*k5)/128 - (21*k5_2)/640 + (77*k3_2*k6)/384 - (7*k4*k6)/128 - (k3*k7)/32 + k8/384); } } } if (factor < 0.0L) { std::ostringstream os; os << "In npstat::AbsCGF1D::saddlepointDensityApproximation: " << "negative density correction factor encountered for " << "approximation order " << order << " at x = " << x << ", shat = " << shat << ". Consider reducing approximation order or density support range."; throw std::invalid_argument(os.str()); } return static_cast(factor)/sqrt(2.0*M_PI*varhat)*exp(Khat - shat*x); } double AbsCGF1D::lugannaniRiceSFApprox(const unsigned nTerms, const double x) const { const double shat = saddlepoint(x); const double Khat = derivative(0U, shat); const double tmp = 2.0*(shat*x - Khat); assert(tmp >= 0.0); const double what = (shat >= 0.0 ? 1.0 : -1.0)*sqrt(tmp); const Gauss1D g(0.0, 1.0); const double mainTerm = g.exceedance(what); const double factor = g.density(what); long double sum = 0.0L; if (nTerms > 0U) { assert(what); const double varhat = derivative(2U, shat); assert(varhat > 0.0); const double uhat = shat*sqrt(varhat); assert(uhat); const double c0 = 1.0/uhat - 1.0/what; sum += c0; if (nTerms > 1U) { const double k3 = standardizedCumulant(3U, shat); const double k4 = standardizedCumulant(4U, shat); const double k3_2 = k3*k3; const double uhat_2 = uhat*uhat; const double uhat_3 = uhat_2*uhat; const double what_3 = what*what*what; const double c1 = (k4/8 - 5*k3_2/24)/uhat - k3/2/uhat_2 - 1.0/uhat_3 + 1.0/what_3; sum += c1; if (nTerms > 2U) { throw std::invalid_argument( "In npstat::AbsCGF1D::luugannaniRiceSFApprox: " "unsupported approximation order (too high)"); } } } return mainTerm + factor*static_cast(sum); } } Index: trunk/npstat/stat/AbsCGF1D.hh =================================================================== --- trunk/npstat/stat/AbsCGF1D.hh (revision 710) +++ trunk/npstat/stat/AbsCGF1D.hh (revision 711) @@ -1,83 +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 */ 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_