Index: trunk/npstat/stat/EdgeworthSeries1D.hh =================================================================== --- trunk/npstat/stat/EdgeworthSeries1D.hh (revision 684) +++ trunk/npstat/stat/EdgeworthSeries1D.hh (revision 685) @@ -1,151 +1,151 @@ #ifndef NPSTAT_EDGEWORTHSERIES1D_HH_ #define NPSTAT_EDGEWORTHSERIES1D_HH_ /*! // \file EdgeworthSeries1D.hh // // \brief Distribution defined by Edgeworth series w.r.t. normal // // Author: I. Volobouev // // June 2019 */ #include #include "npstat/nm/SimpleFunctors.hh" #include "npstat/stat/AbsDistribution1D.hh" #include "npstat/stat/EdgeworthSeriesMethod.hh" namespace npstat { class EdgeworthSeries1D : public AbsDistribution1D { public: /** // Constructor arguments are as follows: // // cumulants -- The vector of cumulants. The first element // of the vector (with index 0) is the first // cumulant, the second element is the second, etc. // // m -- The method used to construct the distribition. // See header EdgeworthSeriesMethod.hh for further // explanations. // // order -- The expansion order in powers of 1/sqrt(n). // O(1/sqrt(n)) is the 0th order, O(1/n) is the first // order, O(1/(n sqrt(n)) is the second order, etc. - // Currently, the highest supported order is 12 for + // Currently, the highest supported order is 14 for // EDGEWORTH_CLASSICAL method in a non-SLR mode, and // 4 in all other situations. // // slrMode -- If true, we will assume that the kth cumulant // vanishes (for k > 2) to O(n^{-k/2}) instead of // the usual O(n^{-(k-2)/2}). This is the case for // the SLR statistic. */ EdgeworthSeries1D(const std::vector& cumulants, EdgeworthSeriesMethod m, unsigned order, bool slrMode = false); inline virtual EdgeworthSeries1D* clone() const {return new EdgeworthSeries1D(*this);} inline virtual ~EdgeworthSeries1D() {} virtual double density(double x) const; virtual double cdf(double x) const; virtual double exceedance(double x) const; virtual double quantile(double x) const; /** // Inverse exceedance function. More precise than "quantile(1 - x)" // for small x. */ double inverseExceedance(double x) const; //@{ /** Inspector */ inline EdgeworthSeriesMethod method() const {return m_;} inline bool slrMode() const {return slrMode_;} inline unsigned order() const {return order_;} //@} /** // This function returns the constructor argument, not the // real cumulant which depends on the expansion order used. // Use the "empiricalCentralMoment" function and then convert // central moments to cumulants if you need real cumulants. // The real mean is 0 if order == 0 and cum(0) if order > 0. */ inline double cum(const unsigned k) const {return cumulants_.at(k);}; /** Mean of the function used to build the expansion */ double edgMean() const; /** Standard deviation of the function used to build the expansion */ double edgStdev() const; /** Edgeworth expansion factor for the density */ double densityFactor(double x) const; /** Edgeworth expansion factor for the CDF */ double cdfFactor(double x) const; /** // Actual central moment of the density // (evaluated numerically for degree > 1) */ double empiricalCentralMoment(unsigned degree) const; //@{ /** Prototype needed for I/O */ inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream&) const; //@} static inline const char* classname() {return "npstat::EdgeworthSeries1D";} static inline unsigned version() {return 1;} static EdgeworthSeries1D* read( const gs::ClassId& id, std::istream& is); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: void hermiteCoeffsClassical(long double* buf, unsigned* maxdeg) const; void hermiteCoeffsNormal(long double* buf, unsigned* maxdeg) const; void hermiteCoeffsSLR(long double* buf, unsigned* maxdeg) const; double normalizedCoord(double x) const; double normalizedSigma() const; unsigned minNCumulants() const; std::vector cumulants_; EdgeworthSeriesMethod m_; unsigned order_; bool slrMode_; }; /** This class is mostly useful for testing purposes */ class EdgeworthSeriesMomentFcn : public Functor1 { public: inline EdgeworthSeriesMomentFcn(const EdgeworthSeries1D& series, const unsigned degree) : series_(series), mean_(0.0), degree_(degree) {if (series_.order()) {mean_ = series_.cum(0);}} inline virtual ~EdgeworthSeriesMomentFcn() {} virtual double operator()(const double&) const; private: const EdgeworthSeries1D& series_; double mean_; unsigned degree_; }; } #endif // NPSTAT_EDGEWORTHSERIES1D_HH_ Index: trunk/npstat/stat/AbsCGF1D.cc =================================================================== --- trunk/npstat/stat/AbsCGF1D.cc (revision 684) +++ trunk/npstat/stat/AbsCGF1D.cc (revision 685) @@ -1,136 +1,162 @@ #include #include #include #include "npstat/stat/AbsCGF1D.hh" #include "npstat/nm/findRootUsingBisections.hh" namespace { class AbsCGF1DDeriv : public npstat::Functor1 { public: inline AbsCGF1DDeriv(const npstat::AbsCGF1D& cgf) : cgf_(cgf) {} inline double operator()(const double& s) const {return cgf_.derivative(1U, s);} private: const npstat::AbsCGF1D& cgf_; }; } namespace npstat { double AbsCGF1D::saddlepoint(const double x) const { 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()/2.01) break; xmax *= 2.0; } 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()/2.01) break; xmin *= 2.0; } 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); 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 > 2U) throw std::invalid_argument( + 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); double factor = 1.0; if (order > 1U) { const double k3 = standardizedCumulant(3U, shat); const double k4 = standardizedCumulant(4U, shat); - factor += (k4/8 - 5*k3*k3/24); + 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); + } + } } return factor/sqrt(2.0*M_PI*varhat)*exp(Khat - shat*x); } }