Index: trunk/npstat/stat/saddlepointDistribution1D.cc =================================================================== --- trunk/npstat/stat/saddlepointDistribution1D.cc (revision 917) +++ trunk/npstat/stat/saddlepointDistribution1D.cc (revision 918) @@ -1,27 +1,27 @@ #include #include "npstat/stat/saddlepointDistribution1D.hh" namespace npstat { CPP11_auto_ptr saddlepointDistribution1D( const AbsCGF1D& cgf, const unsigned approximationOrder, const unsigned interpolationDegree, unsigned nCoords, const double xmin, const double xmax) { if (!(xmin < xmax)) throw std::invalid_argument( "In npstat::saddlepointDistribution1D: invalid support interval"); if (nCoords < 2U) nCoords = 2U; const double range = xmax - xmin; - const double step = range/(nCoords - 1U); const unsigned ncM1 = nCoords - 1U; + const double step = range/ncM1; std::vector density(nCoords); for (unsigned i=0; i(new Tabulated1D( xmin, range, density, interpolationDegree)); } } Index: trunk/npstat/stat/arrayStats.icc =================================================================== --- trunk/npstat/stat/arrayStats.icc (revision 917) +++ trunk/npstat/stat/arrayStats.icc (revision 918) @@ -1,704 +1,704 @@ #include #include #include #include #include #include "npstat/nm/BoxNDScanner.hh" #include "npstat/rng/permutation.hh" namespace npstat { template void arrayCentralMoments(const Numeric* arr, const unsigned long sz, const unsigned maxOrder, long double* moments, long double* momentUncertainties) { if (!sz) throw std::invalid_argument("In npstat::arrayCentralMoments: " "can not process arrays with zero elements"); assert(arr); assert(moments); moments[0] = 1.0L; if (momentUncertainties) momentUncertainties[0] = 0.0L; if (maxOrder) { long double mean = 0.0L; for (unsigned long i=0; i(arr[i]); mean /= sz; if (momentUncertainties) { std::vector mu(2*maxOrder+1); arrayMoments(arr, sz, mean, 2*maxOrder, &mu[0]); for (unsigned k=2; k<=maxOrder; ++k) { const long double mm1 = k > 2U ? mu[k-1] : 0.0L; const long double tmp = k*mm1; const long double v = mu[2*k] - 2*k*mm1*mu[k+1] - mu[k]*mu[k] + tmp*tmp*mu[2]; assert(v >= 0.0L); moments[k] = mu[k]; momentUncertainties[k] = sqrtl(v/sz); } momentUncertainties[1] = sqrtl(mu[2]/sz); } else if (maxOrder > 1U) arrayMoments(arr, sz, mean, maxOrder, moments); moments[1] = mean; } } template void arrayMoments(const Numeric* arr, const unsigned long sz, const long double center, const unsigned maxOrder, long double* moments) { if (!sz) throw std::invalid_argument("In npstat::arrayMoments: " "can not process arrays with zero elements"); assert(arr); assert(moments); moments[0] = 1.0L; for (unsigned k=1; k<=maxOrder; ++k) moments[k] = 0.0L; if (maxOrder) { for (unsigned long i=0; i(arr[i]) - center; long double dk = 1.0L; for (unsigned k=1; k<=maxOrder; ++k) { dk *= delta; moments[k] += dk; } } } for (unsigned k=1; k<=maxOrder; ++k) moments[k] /= sz; } template void arrayCumulants(const Numeric* arr, const unsigned long sz, const unsigned maxOrder, long double* cumulants) { const unsigned orderSupported = 6U; long double s[orderSupported+1U] = {0.0L}; if (maxOrder > orderSupported) throw std::invalid_argument("In npstat::arrayCumulants: order argument " "outside of supported range"); if (!sz || sz < maxOrder) throw std::invalid_argument("In npstat::arrayCumulants: " "insufficient array length"); assert(arr); assert(cumulants); - cumulants[0] = 0.0L; + cumulants[0] = 1.0L; if (maxOrder) { long double mean = 0.0L; for (unsigned long i=0; i(arr[i]); mean /= sz; cumulants[1] = mean; if (maxOrder > 1U) { for (unsigned long i=0; i(arr[i]) - mean; long double dk = delta; for (unsigned k=2U; k<=maxOrder; ++k) { dk *= delta; s[k] += dk; } } cumulants[2] = s[2]/(sz - 1U); if (maxOrder > 2U) cumulants[3] = sz*s[3]/(sz-1U)/(sz-2U); if (maxOrder > 3U) { const long double n = sz; const long double n2 = n*n; cumulants[4] = ((n2 + n)*s[4] - s[2]*s[2]*3*(sz - 1U))/(sz-1U)/(sz-2U)/(sz-3U); if (maxOrder > 4U) { const long double n3 = n2*n; cumulants[5] = ((n3 + 5*n2)*s[5] - s[3]*s[2]*10*(n2 - n))/(sz-1U)/(sz-2U)/(sz-3U)/(sz-4U); if (maxOrder > 5U) { const long double n4 = n2*n2; cumulants[6] = ((n4 + 16*n3 + 11*n2 - 4*n)*s[6] - 15*(n-1)*(n-1)*(n+4)*s[4]*s[2] - 10*(n3 - 2*n2 + 5*n - 4)*s[3]*s[3] + 30*(n2 - 3*n + 2)*s[2]*s[2]*s[2])/(sz-1U)/(sz-2U)/(sz-3U)/(sz-4U)/(sz-5U); } } } } } } template long double arrayMoment(const Numeric* arr, const unsigned long sz, const long double center, const unsigned order) { if (!sz) throw std::invalid_argument("In npstat::arrayMoment: " "can not process arrays with zero elements"); assert(arr); long double sum = 0.0L; switch (order) { case 0U: return 1.0L; case 1U: for (unsigned long i=0; i(arr[i]) - center); break; case 2U: for (unsigned long i=0; i(arr[i]) - center; sum += delta*delta; } break; case 3U: for (unsigned long i=0; i(arr[i]) - center; sum += delta*delta*delta; } break; case 4U: for (unsigned long i=0; i(arr[i]) - center; const long double tmp = delta*delta; sum += tmp*tmp; } break; default: for (unsigned long i=0; i(arr[i]) - center; sum += powl(delta, order); } } return sum/sz; } template inline void arrayStats(const Data* arr, const unsigned long sz, double* pmean, double* pstdev, double* pskew, double* pkurt) { if (pmean || pstdev || pskew || pkurt) { if (!sz) throw std::invalid_argument("In npstat::arrayStats: " "can not process arrays with zero elements"); assert(arr); if (pstdev) *pstdev = 0.0; if (pskew) *pskew = 0.0; if (pkurt) *pkurt = 3.0; if (sz == 1UL) { if (pmean) *pmean = static_cast(arr[0]); } else { long double sm = 0.0L; for (unsigned long i=0; i(mean); const bool b = (pskew && sz > 2UL) || (pkurt && sz > 3UL); if (pstdev || b) { long double sumsq = 0.0L, skewsum = 0.0L, kursum = 0.0L; for (unsigned long i=0; i 0.0L) { if (pstdev) *pstdev = static_cast(sqrtl(K2)); if (pskew && sz > 2UL) { const long double K3 = sz/NM1/(sz - 2UL)*skewsum; *pskew = static_cast(K3/powl(K2,1.5L)); } if (pkurt && sz > 3UL) { const long double g2 = kursum*sz/sumsq/sumsq-3.0L; *pkurt = static_cast(3.0L + NM1/(sz-2UL)/(sz-3UL)*((sz+1UL)*g2 + 6.0L)); } } } } } } template inline void arrayCoordMean(const Array& a, const BoxND& limits, double* mean, const unsigned storageLength) { const unsigned lenMean = a.rank(); if (lenMean) { long double meanAcc[CHAR_BIT*sizeof(unsigned long)] = {0.0L,}; if (lenMean > CHAR_BIT*sizeof(unsigned long)) throw std::out_of_range("In npstat::arrayCoordMean: " "array dimensionality is too large"); assert(mean); if (storageLength < lenMean) throw std::invalid_argument("In npstat::arrayCoordMean: " "result buffer is too small"); long double weightSum = 0.0L; for (BoxNDScanner scanner(limits, a.shape()); scanner.isValid(); ++scanner) { scanner.getCoords(mean, lenMean); const long double v = a.linearValue(scanner.state()); weightSum += v; for (unsigned i=0; i(meanAcc[i]/weightSum); } } template inline void arrayCoordCovariance(const Array& a, const BoxND& limits, Matrix* covmat) { const unsigned aRank = a.rank(); if (aRank) { assert(covmat); if (covmat->nRows() != aRank) throw std::invalid_argument("In npstat::arrayCoordCovariance: " "incompatible number of rows in the result matrix"); if (covmat->nColumns() != aRank) throw std::invalid_argument("In npstat::arrayCoordCovariance: " "incompatible number of columns in the result matrix"); double mean[CHAR_BIT*sizeof(unsigned long)]; arrayCoordMean(a, limits, mean, aRank); long double weightSum = 0.0L, weightSumSq = 0.0L; const unsigned matLen = aRank*(aRank+1U)/2; std::vector ldbuf(matLen, 0.0L); long double* acc = &ldbuf[0]; double coords[CHAR_BIT*sizeof(unsigned long)]; for (BoxNDScanner scanner(limits, a.shape()); scanner.isValid(); ++scanner) { scanner.getCoords(coords, aRank); const long double v = a.linearValue(scanner.state()); weightSum += v; weightSumSq += v*v; unsigned ilin = 0; for (unsigned i=0; i(covmat->data()); const unsigned len = covmat->length(); for (unsigned i=0; i( acc[ilin]/weightSum*(effM1 + 1.0L)/effM1); (*covmat)[i][j] = cov; (*covmat)[j][i] = cov; } } } } // The formulae used below assume that the effective number // of events is given by the ratio of squared sum of weights // to sum of weights squared. Once the effective number of events // is known, it is used with the standard population estimates // consistent with SAS (see, for example, the article "Comparing // measures of sample skewness and kurtosis", by D.N. Joanes // and C.A. Gill, The Statistician, Vol 47, pp 183-189 (1998). template inline void arrayShape1D(const Array& a, const double xmin, const double xmax, double* pmean, double* pstdev, double* pskewness, double* pkurtosis) { typedef typename Array::value_type Numeric; if (a.rank() != 1U) throw std::invalid_argument( "In npstat::arrayShape1D: array dimensionality must be 1"); if (xmin >= xmax) throw std::invalid_argument( "In npstat::arrayShape1D: invalid interval specification"); if (pmean || pstdev || pskewness || pkurtosis) { const Numeric* data = a.data(); const Numeric zero = static_cast(0); const unsigned long nbins = a.length(); const double binWidth = (xmax - xmin)/nbins; long double wsum = 0.0L, xsum = 0.0L; for (unsigned long i=0; i(xsum/wsum); if (pmean) *pmean = mean; if (pstdev || pskewness || pkurtosis) { long double wsumsq = 0.0L, varsum = 0.0L; const double shiftedMin = xmin - mean; for (unsigned long i=0; i(sqrtl(K2)); if (pskewness || pkurtosis) { long double skewsum = 0.0L, kursum = 0.0L; for (unsigned long i=0; i 0.0L) { const long double m3 = skewsum/wsum; const long double K3 = effN*effN/effM1/effM2*m3; *pskewness = static_cast(K3/powl(K2,1.5L)); } else *pskewness = 0.0; } if (pkurtosis) { const long double effM3 = effM2 - 1.0L; if (effM3 > 0.0L) { const long double m4 = kursum/wsum; const long double g2 = m4/m2/m2 - 3.0L; *pkurtosis = static_cast(3.0L + effM1/effM2/effM3*((effN+1.0L)*g2 + 6.0L)); } else *pkurtosis = 3.0; } } } } } template inline void arrayQuantiles1D(const Numeric* data, const unsigned long len, const double xmin, const double xmax, const double* qvalues, double* quantiles, const unsigned nqvalues) { if (nqvalues == 0U) return; if (!len) throw std::invalid_argument("In npstat::arrayQuantiles1D: " "can not process arrays with zero elements"); assert(data); if (xmin > xmax) throw std::invalid_argument("In npstat::arrayQuantiles1D: " "invalid interval specification"); assert(qvalues); assert(quantiles); for (unsigned i=0; i= 0.0 && qvalues[i] <= 1.0)) throw std::domain_error("In npstat::arrayQuantiles1D: " "cdf argument outside of [0, 1] interval"); if (xmin == xmax) { for (unsigned i=0; i zero) { if (!positiveElementFound) { positiveElementFound = true; idlo = id; } sum += value; idhi = id; } } if (!positiveElementFound) throw std::invalid_argument( "In npstat::arrayQuantiles1D: all weights are zero"); const double xlo = xmin + idlo*step; const double xhi = idhi == len - 1UL ? xmax : xmin + (idhi+1UL)*step; const unsigned long lenm1 = len - 1UL; unsigned long scanbin = idlo; long double sumbelow = 0.0L; long double sumabove = data[idlo]; for (unsigned i=0; i= sumabove && scanbin < lenm1) { sumbelow = sumabove; sumabove += data[++scanbin]; } const double dbin = (sumneeded - sumbelow)/data[scanbin]; double qval = xmin + (scanbin + dbin)*step; if (qval > xmax) qval = xmax; quantiles[i] = qval; } } } template inline double arrayEntropy(const Numeric* p, const unsigned long len, const bool normalize) { const Numeric zero = Numeric(); bool havePositive = false; long double entropy = 0.0L; long double sum = 1.0L; if (normalize) { sum = std::accumulate(p, p+len, 0.0L); if (sum <= 0.0L) throw std::invalid_argument("In npstat::arrayEntropy: sum of" " array elements is not positive"); } for (unsigned long i=0; i zero) { havePositive = true; const long double prob = p[i]; entropy -= prob*std::log(prob/sum); } } if (!havePositive) throw std::invalid_argument("In npstat::arrayEntropy: there are no" " positive array elements"); entropy /= sum; return entropy/len; } // C++ does not support partial template specialization for // functions, only for classes namespace Private { template struct PoissonLogLikelihood { inline static double calculate(const Real* means, const Numeric* counts, const unsigned long len) { long double sum = 0.0L; for (unsigned long i=0; i(counts[i])); if (dc < 0.0) throw std::invalid_argument( "In npstat::poissonLogLikelihood: " "this function can not handle negative counts"); if (dc > static_cast(ULONG_MAX)) throw std::invalid_argument( "In npstat::poissonLogLikelihood: " "count is too large"); const unsigned long cnt = static_cast(dc); const double m = means[i]; if (cnt) { if (m <= 0.0) throw std::invalid_argument( "In npstat::poissonLogLikelihood: " "non-positive mean encountered with counts present"); sum += (cnt*log(m) - logfactorial(cnt) - m); } else sum -= m; } return sum; } }; template struct PoissonLogLikelihood { inline static double calculate(const Real* means, const unsigned* counts, const unsigned long len) { long double sum = 0.0L; for (unsigned long i=0; i struct PoissonLogLikelihood { inline static double calculate(const Real* means, const unsigned long* counts, const unsigned long len) { long double sum = 0.0L; for (unsigned long i=0; i inline double poissonLogLikelihood(const Real* means, const Numeric* counts, const unsigned long len) { double logli = 0.0; if (len) { assert(means); assert(counts); logli = Private::PoissonLogLikelihood::calculate( means, counts, len); } return logli; } } Index: trunk/npstat/stat/JohnsonCurves.hh =================================================================== --- trunk/npstat/stat/JohnsonCurves.hh (revision 917) +++ trunk/npstat/stat/JohnsonCurves.hh (revision 918) @@ -1,216 +1,217 @@ #ifndef NPSTAT_JOHNSONCURVES_HH_ #define NPSTAT_JOHNSONCURVES_HH_ /*! // \file JohnsonCurves.hh // // \brief Johnson frequency curves // // The S_u distribution is rather easy to fit by moments, and the // fitting procedure works well. The S_b is much harder to fit, and // the implementation may be rough at the moment. It is, however, // better than any previously published algorithm. // // Author: I. Volobouev // // April 2010 */ #include "npstat/stat/Distribution1DFactory.hh" namespace npstat { /** Johnson S_u (unbounded) curve */ class JohnsonSu : public AbsScalableDistribution1D { public: JohnsonSu(double location, double scale, double skewness, double kurtosis); inline virtual JohnsonSu* clone() const {return new JohnsonSu(*this);} inline virtual ~JohnsonSu() {} inline double skewness() const {return skew_;} inline double kurtosis() const {return kurt_;} inline bool isValid() const {return isValid_;} inline double getDelta() const {return delta_;} inline double getLambda() const {return lambda_;} inline double getGamma() const {return gamma_;} inline double getXi() const {return xi_;} // Methods needed for I/O inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::JohnsonSu";} static inline unsigned version() {return 1;} static JohnsonSu* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; JohnsonSu(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} // The I/O helper constructor inline JohnsonSu(double location, double scale) : AbsScalableDistribution1D(location, scale) {} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; void initialize(); double skew_; double kurt_; double delta_; double lambda_; double gamma_; double xi_; bool isValid_; }; /** Johnson S_b (bounded) curve */ class JohnsonSb : public AbsScalableDistribution1D { public: JohnsonSb(double location, double scale, double skewness, double kurtosis); inline virtual JohnsonSb* clone() const {return new JohnsonSb(*this);} inline virtual ~JohnsonSb() {} inline double skewness() const {return skew_;} inline double kurtosis() const {return kurt_;} inline bool isValid() const {return isValid_;} inline double getDelta() const {return delta_;} inline double getLambda() const {return lambda_;} inline double getGamma() const {return gamma_;} inline double getXi() const {return xi_;} static bool fitParameters(double skewness, double kurtosis, double *gamma, double *delta, double *lambda, double *xi); // Methods needed for I/O inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::JohnsonSb";} static inline unsigned version() {return 1;} static JohnsonSb* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; JohnsonSb(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} // The I/O helper constructor inline JohnsonSb(double location, double scale) : AbsScalableDistribution1D(location, scale) {} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; double skew_; double kurt_; double delta_; double lambda_; double gamma_; double xi_; bool isValid_; }; /** This class selects an appropriate Johnson curve automatically */ class JohnsonSystem : public AbsScalableDistribution1D { public: enum CurveType { GAUSSIAN = 0, LOGNORMAL, SU, SB, INVALID }; // std::invalid_argument will be thrown if the combination // of skewness and kurtosis arguments is impossible JohnsonSystem(double location, double scale, double skewness, double kurtosis); JohnsonSystem(const JohnsonSystem&); JohnsonSystem& operator=(const JohnsonSystem&); virtual ~JohnsonSystem(); inline virtual JohnsonSystem* clone() const {return new JohnsonSystem(*this);} inline double skewness() const {return skew_;} inline double kurtosis() const {return kurt_;} inline CurveType curveType() const {return curveType_;} + inline bool isValid() const {return !(curveType_ == INVALID);} static CurveType select(double skewness, double kurtosis); // Methods needed for I/O inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::JohnsonSystem";} static inline unsigned version() {return 1;} static JohnsonSystem* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; JohnsonSystem(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} // The I/O helper constructor inline JohnsonSystem(double location, double scale) : AbsScalableDistribution1D(location, scale), fcn_(0) {} void initialize(); inline double unscaledDensity(double x) const {return fcn_->density(x);} inline double unscaledCdf(double x) const {return fcn_->cdf(x);} inline double unscaledExceedance(double x) const {return fcn_->exceedance(x);} inline double unscaledQuantile(double x) const {return fcn_->quantile(x);} AbsScalableDistribution1D* fcn_; double skew_; double kurt_; CurveType curveType_; }; } #endif // NPSTAT_JOHNSONCURVES_HH_ Index: trunk/npstat/stat/EdgeworthSeries1D.cc =================================================================== --- trunk/npstat/stat/EdgeworthSeries1D.cc (revision 917) +++ trunk/npstat/stat/EdgeworthSeries1D.cc (revision 918) @@ -1,644 +1,646 @@ #include #include "npstat/nm/MathUtils.hh" #include "npstat/nm/SpecialFunctions.hh" #include "npstat/nm/findRootUsingBisections.hh" #include "npstat/nm/GaussHermiteQuadrature.hh" #include "npstat/stat/EdgeworthSeries1D.hh" #include "npstat/stat/distributionReadError.hh" #include "geners/binaryIO.hh" #include "geners/vectorIO.hh" #define EDGEWORTH_MAXORDER 14U #define MAXCOEFFS (EDGEWORTH_MAXORDER*3U) #define LSQR2 1.41421356237309504880169L #define SQRTWOPIL 2.50662827463100050241577L static long double ldgexceedance(const long double x) { return erfcl(x/LSQR2)/2.0L; } static long double ldgcdf(const long double x) { return erfcl(-x/LSQR2)/2.0L; } namespace npstat { EdgeworthSeries1D::EdgeworthSeries1D( const std::vector& cumulants, const EdgeworthSeriesMethod m, const unsigned i_order, const bool i_slrMode) : cumulants_(cumulants), m_(m), order_(i_order), slrMode_(i_slrMode) { unsigned maxOrder = 4U; if (m == EDGEWORTH_CLASSICAL && !slrMode_) maxOrder = EDGEWORTH_MAXORDER; if (order_ > maxOrder) throw std::invalid_argument( "In npstat::EdgeworthSeries1D constructor: " "order parameter out of range"); if (order_) { const unsigned nMin = minNCumulants(); if (cumulants_.size() < nMin) throw std::invalid_argument( "In npstat::EdgeworthSeries1D constructor: " "not enough cumulants provided"); if (nMin > 1U) if (cumulants_[1] <= 0.0) throw std::invalid_argument( "In npstat::EdgeworthSeries1D constructor: " "variance must be positive"); } } double EdgeworthSeries1D::edgMean() const { if (minNCumulants() && m_ == EDGEWORTH_CLASSICAL) return cumulants_.at(0); else return 0.0; } double EdgeworthSeries1D::edgStdev() const { if (minNCumulants() > 1) return normalizedSigma(); else return 1.0; } unsigned EdgeworthSeries1D::minNCumulants() const { if (slrMode_) return order_; else return order_ ? order_ + 2U : 0U; } bool EdgeworthSeries1D::isEqual(const AbsDistribution1D& o) const { const EdgeworthSeries1D& r = static_cast(o); return cumulants_ == r.cumulants_ && m_ == r.m_ && order_ == r.order_ && slrMode_ == r.slrMode_; } bool EdgeworthSeries1D::write(std::ostream& os) const { const unsigned m = static_cast(m_); const unsigned char slr = slrMode_; gs::write_pod(os, m); gs::write_pod(os, order_); gs::write_pod(os, slr); return gs::write_item(os, cumulants_, false) && !os.fail(); } EdgeworthSeries1D* EdgeworthSeries1D::read( const gs::ClassId& id, std::istream& in) { static const gs::ClassId myClassId( gs::ClassId::makeId()); myClassId.ensureSameId(id); unsigned m, order; unsigned char slr; gs::read_pod(in, &m); gs::read_pod(in, &order); gs::read_pod(in, &slr); std::vector cumulants; gs::restore_item(in, &cumulants, false); if (in.fail()) { distributionReadError(in, classname()); return 0; } else return new EdgeworthSeries1D( cumulants, static_cast(m), order, slr); } double EdgeworthSeries1D::densityFactor(const double xin) const { long double coeffs[MAXCOEFFS+1]; unsigned maxdeg = 0; if (order_) { if (slrMode_) hermiteCoeffsSLR(coeffs+1, &maxdeg); else hermiteCoeffsNormal(coeffs+1, &maxdeg); ++maxdeg; } coeffs[0] = 0.0L; + + // CHECK THIS!!! Do we really need to divide the series by normalizedSigma()? return 1.0 + hermiteSeriesSumProb(coeffs, maxdeg, normalizedCoord(xin))/normalizedSigma(); } double EdgeworthSeries1D::normalizedSigma() const { double sigma = 1.0; if (order_ > 1U && m_ == EDGEWORTH_CLASSICAL) sigma = sqrt(cumulants_.at(1)); return sigma; } double EdgeworthSeries1D::density(const double xin) const { const double x = normalizedCoord(xin); return exp(-x*x/2)/SQRTWOPIL/normalizedSigma()*densityFactor(xin); } double EdgeworthSeries1D::cdf(const double xin) const { const double z = normalizedCoord(xin); return ldgcdf(z) - expl(-z*z/2.0)/SQRTWOPIL/normalizedSigma()*cdfFactor(xin); } double EdgeworthSeries1D::exceedance(const double xin) const { const double z = normalizedCoord(xin); return ldgexceedance(z) + expl(-z*z/2.0)/SQRTWOPIL/normalizedSigma()*cdfFactor(xin); } double EdgeworthSeries1D::quantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::EdgeworthSeries1D::quantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0 || r1 == 1.0) return inverseGaussCdf(r1); const double x0 = inverseGaussCdf(0.0); const double x1 = inverseGaussCdf(1.0); CdfFunctor1D fcn(*this); double q = 0.0; if (!findRootUsingBisections(fcn, r1, x0, x1, 1.0e-15, &q)) throw std::runtime_error("In npstat::EdgeworthSeries1D::quantile: " "root finding failed"); return q; } double EdgeworthSeries1D::inverseExceedance(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::EdgeworthSeries1D::inverseExceedance: " "exceedance argument outside of [0, 1] interval"); if (r1 == 0.0 || r1 == 1.0) return inverseGaussCdf(1.0 - r1); const double x0 = inverseGaussCdf(0.0); const double x1 = inverseGaussCdf(1.0); ExceedanceFunctor1D fcn(*this); double q = 0.0; if (!findRootUsingBisections(fcn, r1, x0, x1, 1.0e-15, &q)) throw std::runtime_error("In npstat::EdgeworthSeries1D::inverseExceedance: " "root finding failed"); return q; } double EdgeworthSeries1D::normalizedCoord(const double xin) const { double x = xin; if (order_) { switch (m_) { case EDGEWORTH_SEVERINI: break; case EDGEWORTH_CLASSICAL: { double sigma = 1.0; if (order_ > 1U) sigma = sqrt(cumulants_.at(1)); x = (xin - cumulants_.at(0))/sigma; } break; default: assert(!"Incomplete swich statement. This is a bug. Please report."); } } return x; } double EdgeworthSeries1D::cdfFactor(const double xin) const { if (order_) { long double coeffs[MAXCOEFFS]; unsigned maxdeg; if (slrMode_) hermiteCoeffsSLR(coeffs, &maxdeg); else hermiteCoeffsNormal(coeffs, &maxdeg); return hermiteSeriesSumProb(coeffs, maxdeg, normalizedCoord(xin)); } else return 0.0; } double EdgeworthSeries1D::empiricalCentralMoment(const unsigned k) const { switch (k) { case 0U: return 1.0; case 1U: return 0.0; default: if (k > 200U) throw std::invalid_argument("In npstat::EdgeworthSeries1D::empiricalCentralMoment: " "moment order is too high"); GaussHermiteQuadrature quad(128); EdgeworthSeriesMomentFcn fcn(*this, k); return quad.integrateProb(edgMean(), edgStdev(), fcn); } } void EdgeworthSeries1D::hermiteCoeffsSLR(long double* coeffs, unsigned* maxdeg) const { for (unsigned i=0; i= order_); *maxdeg = order_ - 1U; assert(MAXCOEFFS > *maxdeg); const double k1 = m_ == EDGEWORTH_SEVERINI ? cumulants_[0] : 0.0; coeffs[0] = k1; if (order_ > 1U) { const double k2 = m_ == EDGEWORTH_SEVERINI ? cumulants_[1] : 1.0; const double k2m1 = k2 - 1.0; coeffs[1] = (k1*k1 + k2m1)/2.0; if (order_ > 2U) { const double k3 = cumulants_[2]; coeffs[2] = (k1*(k1*k1 + 3*k2m1) + k3)/6; if (order_ > 3U) { const double k4 = cumulants_[3]; const double k1sq = k1*k1; coeffs[3] = (k1sq*k1sq + 6*k1sq*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24; if (order_ > 4U) // We should not be here assert(0); } if (m_ == EDGEWORTH_CLASSICAL) { const double sigma = sqrt(cumulants_[1]); for (unsigned i=2; i<=*maxdeg; ++i) coeffs[i] /= powl(sigma, i); } } } } void EdgeworthSeries1D::hermiteCoeffsClassical(long double* coeffs, unsigned* maxdeg) const { for (unsigned i=0; i 2U); const double kappa3 = cumulants_[2]; *maxdeg = 2U; coeffs[2] = kappa3/6.0; if (order_ >= 2U) { *maxdeg = 5U; const double kappa4 = cumulants_[3]; const double kappa3_2 = kappa3*kappa3; coeffs[3] += kappa4/24; coeffs[5] += kappa3_2/72; if (order_ >= 3U) { *maxdeg = 8U; const double kappa5 = cumulants_[4]; const double kappa3_3 = kappa3_2*kappa3; coeffs[4] += kappa5/120; coeffs[6] += (kappa3*kappa4)/144; coeffs[8] += kappa3_3/1296; if (order_ >= 4U) { *maxdeg = 11U; const double kappa6 = cumulants_[5]; const double kappa3_4 = kappa3_3*kappa3; const double kappa4_2 = kappa4*kappa4; coeffs[5] += kappa6/720; coeffs[7] += kappa4_2/1152 + (kappa3*kappa5)/720; coeffs[9] += (kappa3_2*kappa4)/1728; coeffs[11] += kappa3_4/31104; if (order_ >= 5U) { *maxdeg = 14U; const double kappa7 = cumulants_[6]; const double kappa3_5 = kappa3_4*kappa3; coeffs[6] += kappa7/5040; coeffs[8] += (kappa4*kappa5)/2880 + (kappa3*kappa6)/4320; coeffs[10] += (kappa3*kappa4_2)/6912 + (kappa3_2*kappa5)/8640; coeffs[12] += (kappa3_3*kappa4)/31104; coeffs[14] += kappa3_5/933120; if (order_ >= 6U) { *maxdeg = 17U; const double kappa8 = cumulants_[7]; const double kappa5_2 = kappa5*kappa5; const double kappa4_3 = kappa4_2*kappa4; const double kappa3_6 = kappa3_5*kappa3; coeffs[7] += kappa8/40320; coeffs[9] += kappa5_2/28800 + (kappa4*kappa6)/17280 + (kappa3*kappa7)/30240; coeffs[11] += kappa4_3/82944 + (kappa3*kappa4*kappa5)/17280 + (kappa3_2*kappa6)/51840; coeffs[13] += (kappa3_2*kappa4_2)/82944 + (kappa3_3*kappa5)/155520; coeffs[15] += (kappa3_4*kappa4)/746496; coeffs[17] += kappa3_6/33592320; if (order_ >= 7U) { *maxdeg = 20U; const double kappa9 = cumulants_[8]; const double kappa3_7 = kappa3_6*kappa3; coeffs[8] += kappa9/362880.0; coeffs[10] += (kappa5*kappa6)/86400 + (kappa4*kappa7)/120960 + (kappa3*kappa8)/241920.0; coeffs[12] += (kappa4_2*kappa5)/138240 + (kappa3*kappa5_2)/172800 + (kappa3*kappa4*kappa6)/103680 + (kappa3_2*kappa7)/362880.0; coeffs[14] += (kappa3*kappa4_3)/497664 + (kappa3_2*kappa4*kappa5)/207360 + (kappa3_3*kappa6)/933120.0; coeffs[16] += (kappa3_3*kappa4_2)/1492992 + (kappa3_4*kappa5)/3732480.0; coeffs[18] += (kappa3_5*kappa4)/22394880.0; coeffs[20] += kappa3_7/1410877440.0; if (order_ >= 8U) { *maxdeg = 23U; const double kappa10 = cumulants_[9]; const double kappa6_2 = kappa6*kappa6; const double kappa3_8 = kappa3_7*kappa3; const double kappa4_4 = kappa4_3*kappa4; coeffs[9] += kappa10/3628800.0; coeffs[11] += kappa6_2/1036800 + (kappa5*kappa7)/604800 + (kappa4*kappa8)/967680 + (kappa3*kappa9)/2177280.0; coeffs[13] += (kappa4*kappa5_2)/691200 + (kappa4_2*kappa6)/829440 + (kappa3*kappa5*kappa6)/518400 + (kappa3*kappa4*kappa7)/725760 + (kappa3_2*kappa8)/2903040.0; coeffs[15] += kappa4_4/7962624 + (kappa3*kappa4_2*kappa5)/829440 + (kappa3_2*kappa5_2)/2073600 + (kappa3_2*kappa4*kappa6)/1244160 + (kappa3_3*kappa7)/6531840.0; coeffs[17] += (kappa3_2*kappa4_3)/5971968 + (kappa3_3*kappa4*kappa5)/3732480 + (kappa3_4*kappa6)/22394880.0; coeffs[19] += (kappa3_4*kappa4_2)/35831808 + (kappa3_5*kappa5)/111974400.0; coeffs[21] += (kappa3_6*kappa4)/806215680.0; coeffs[23] += kappa3_8/67722117120.0; if (order_ >= 9U) { *maxdeg = 26U; const double kappa11 = cumulants_[10]; const double kappa5_3 = kappa5_2*kappa5; const double kappa3_9 = kappa3_8*kappa3; coeffs[10] += kappa11/39916800.0; coeffs[12] += (kappa6*kappa7)/3628800.0 + (kappa5*kappa8)/4838400.0 + (kappa4*kappa9)/8709120.0 + (kappa3*kappa10)/21772800.0; coeffs[14] += kappa5_3/10368000.0 + (kappa4*kappa5*kappa6)/2073600.0 + (kappa3*kappa6_2)/6220800.0 + (kappa4_2*kappa7)/5806080.0 + (kappa3*kappa5*kappa7)/3628800.0 + (kappa3*kappa4*kappa8)/5806080.0 + (kappa3_2*kappa9)/26127360.0; coeffs[16] += (kappa4_3*kappa5)/9953280.0 + (kappa3*kappa4*kappa5_2)/4147200.0 + (kappa3*kappa4_2*kappa6)/4976640.0 + (kappa3_2*kappa5*kappa6)/6220800.0 + (kappa3_2*kappa4*kappa7)/8709120.0 + (kappa3_3*kappa8)/52254720.0; coeffs[18] += (kappa3*kappa4_4)/47775744 + (kappa3_2*kappa4_2*kappa5)/9953280.0 + (kappa3_3*kappa5_2)/37324800.0 + (kappa3_3*kappa4*kappa6)/22394880.0 + (kappa3_4*kappa7)/156764160.0; coeffs[20] += (kappa3_3*kappa4_3)/107495424 + (kappa3_4*kappa4*kappa5)/89579520.0 + (kappa3_5*kappa6)/671846400.0; coeffs[22] += (kappa3_5*kappa4_2)/1074954240.0 + (kappa3_6*kappa5)/4031078400.0; coeffs[24] += (kappa3_7*kappa4)/33861058560.0; coeffs[26] += kappa3_9/3656994324480.0; if (order_ >= 10U) { *maxdeg = 29U; const double kappa12 = cumulants_[11]; const double kappa7_2 = kappa7*kappa7; const double kappa3_10 = kappa3_9*kappa3; const double kappa4_5 = kappa4_4*kappa4; coeffs[11] += kappa12/479001600.0; coeffs[13] += kappa7_2/50803200.0 + (kappa6*kappa8)/29030400.0 + (kappa5*kappa9)/43545600.0 + (kappa4*kappa10)/87091200.0 + (kappa3*kappa11)/239500800.0; coeffs[15] += (kappa5_2*kappa6)/20736000.0 + (kappa4*kappa6_2)/24883200.0 + (kappa4*kappa5*kappa7)/14515200.0 + (kappa3*kappa6*kappa7)/21772800.0 + (kappa4_2*kappa8)/46448640.0 + (kappa3*kappa5*kappa8)/29030400.0 + (kappa3*kappa4*kappa9)/52254720.0 + (kappa3_2*kappa10)/261273600.0; coeffs[17] += (kappa4_2*kappa5_2)/33177600.0 + (kappa3*kappa5_3)/62208000.0 + (kappa4_3*kappa6)/59719680.0 + (kappa3*kappa4*kappa5*kappa6)/12441600.0 + (kappa3_2*kappa6_2)/74649600.0 + (kappa3*kappa4_2*kappa7)/34836480.0 + (kappa3_2*kappa5*kappa7)/43545600.0 + (kappa3_2*kappa4*kappa8)/69672960.0 + (kappa3_3*kappa9)/470292480.0; coeffs[19] += kappa4_5/955514880.0 + (kappa3*kappa4_3*kappa5)/59719680.0 + (kappa3_2*kappa4*kappa5_2)/49766400.0 + (kappa3_2*kappa4_2*kappa6)/59719680.0 + (kappa3_3*kappa5*kappa6)/111974400.0 + (kappa3_3*kappa4*kappa7)/156764160.0 + (kappa3_4*kappa8)/1254113280.0; coeffs[21] += (kappa3_2*kappa4_4)/573308928 + (kappa3_3*kappa4_2*kappa5)/179159040.0 + (kappa3_4*kappa5_2)/895795200.0 + (kappa3_4*kappa4*kappa6)/537477120.0 + (kappa3_5*kappa7)/4702924800.0; coeffs[23] += (kappa3_4*kappa4_3)/2579890176.0 + (kappa3_5*kappa4*kappa5)/2687385600.0 + (kappa3_6*kappa6)/24186470400.0; coeffs[25] += (kappa3_6*kappa4_2)/38698352640.0 + (kappa3_7*kappa5)/169305292800.0; coeffs[27] += (kappa3_8*kappa4)/1625330810880.0; coeffs[29] += kappa3_10/219419659468800.0; if (order_ >= 11U) { *maxdeg = 32U; const double kappa13 = cumulants_[12]; const double kappa3_11 = kappa3_10*kappa3; coeffs[12] += kappa13/6227020800.0L; coeffs[14] += (kappa7*kappa8)/203212800.0L + (kappa6*kappa9)/261273600.0L + (kappa5*kappa10)/435456000.0L + (kappa4*kappa11)/958003200.0L + (kappa3*kappa12)/2874009600.0L; coeffs[16] += (kappa5*kappa6_2)/124416000.0L + (kappa5_2*kappa7)/145152000.0L + (kappa4*kappa6*kappa7)/87091200.0L + (kappa3*kappa7_2)/304819200.0L + (kappa4*kappa5*kappa8)/116121600.0L + (kappa3*kappa6*kappa8)/174182400.0L + (kappa4_2*kappa9)/418037760.0L + (kappa3*kappa5*kappa9)/261273600.0L + (kappa3*kappa4*kappa10)/522547200.0L + (kappa3_2*kappa11)/2874009600.0L; coeffs[18] += (kappa4*kappa5_3)/248832000.0L + (kappa4_2*kappa5*kappa6)/99532800.0L + (kappa3*kappa5_2*kappa6)/124416000.0L + (kappa3*kappa4*kappa6_2)/149299200.0L + (kappa4_3*kappa7)/418037760.0L + (kappa3*kappa4*kappa5*kappa7)/87091200.0L + (kappa3_2*kappa6*kappa7)/261273600.0L + (kappa3*kappa4_2*kappa8)/278691840.0L + (kappa3_2*kappa5*kappa8)/348364800.0L + (kappa3_2*kappa4*kappa9)/627056640.0L + (kappa3_3*kappa10)/4702924800.0L; coeffs[20] += (kappa4_4*kappa5)/955514880.0L + (kappa3*kappa4_2*kappa5_2)/199065600.0L + (kappa3_2*kappa5_3)/746496000.0L + (kappa3*kappa4_3*kappa6)/358318080.0L + (kappa3_2*kappa4*kappa5*kappa6)/149299200.0L + (kappa3_3*kappa6_2)/1343692800.0L + (kappa3_2*kappa4_2*kappa7)/418037760.0L + (kappa3_3*kappa5*kappa7)/783820800.0L + (kappa3_3*kappa4*kappa8)/1254113280.0L + (kappa3_4*kappa9)/11287019520.0L; coeffs[22] += (kappa3*kappa4_5)/5733089280.0L + (kappa3_2*kappa4_3*kappa5)/716636160.0L + (kappa3_3*kappa4*kappa5_2)/895795200.0L + (kappa3_3*kappa4_2*kappa6)/1074954240.0L + (kappa3_4*kappa5*kappa6)/2687385600.0L + (kappa3_4*kappa4*kappa7)/3762339840.0L + (kappa3_5*kappa8)/37623398400.0L; coeffs[24] += (kappa3_3*kappa4_4)/10319560704 + (kappa3_4*kappa4_2*kappa5)/4299816960.0L + (kappa3_5*kappa5_2)/26873856000.0L + (kappa3_5*kappa4*kappa6)/16124313600.0L + (kappa3_6*kappa7)/169305292800.0L; coeffs[26] += (kappa3_5*kappa4_3)/77396705280.0L + (kappa3_6*kappa4*kappa5)/96745881600.0L + (kappa3_7*kappa6)/1015831756800.0L; coeffs[28] += (kappa3_7*kappa4_2)/1625330810880.0L + (kappa3_8*kappa5)/8126654054400.0L; coeffs[30] += (kappa3_9*kappa4)/87767863787520.0L; coeffs[32] += kappa3_11/14481697524940800.0L; if (order_ >= 12U) { *maxdeg = 35U; const double kappa14 = cumulants_[13]; const double kappa8_2 = kappa8*kappa8; const double kappa6_3 = kappa6_2*kappa6; const double kappa3_12 = kappa3_11*kappa3; const double kappa5_4 = kappa5_3*kappa5; const double kappa4_6 = kappa4_5*kappa4; coeffs[13] += kappa14/87178291200.0L; coeffs[15] += kappa8_2/3251404800.0L + (kappa7*kappa9)/1828915200.0L + (kappa6*kappa10)/2612736000.0L + (kappa5*kappa11)/4790016000.0L + (kappa4*kappa12)/11496038400.0L + (kappa3*kappa13)/37362124800.0L; coeffs[17] += kappa6_3/2239488000.0L + (kappa5*kappa6*kappa7)/435456000.0L + (kappa4*kappa7_2)/1219276800.0L + (kappa5_2*kappa8)/1161216000.0L + (kappa4*kappa6*kappa8)/696729600.0L + (kappa3*kappa7*kappa8)/1219276800.0L + (kappa4*kappa5*kappa9)/1045094400.0L + (kappa3*kappa6*kappa9)/1567641600.0L + (kappa4_2*kappa10)/4180377600.0L + (kappa3*kappa5*kappa10)/2612736000.0L + (kappa3*kappa4*kappa11)/5748019200.0L + (kappa3_2*kappa12)/34488115200.0L; coeffs[19] += kappa5_4/4976640000.0L + (kappa4*kappa5_2*kappa6)/497664000.0L + (kappa4_2*kappa6_2)/1194393600.0L + (kappa3*kappa5*kappa6_2)/746496000.0L + (kappa4_2*kappa5*kappa7)/696729600.0L + (kappa3*kappa5_2*kappa7)/870912000.0L + (kappa3*kappa4*kappa6*kappa7)/522547200.0L + (kappa3_2*kappa7_2)/3657830400.0L + (kappa4_3*kappa8)/3344302080.0L + (kappa3*kappa4*kappa5*kappa8)/696729600.0L + (kappa3_2*kappa6*kappa8)/2090188800.0L + (kappa3*kappa4_2*kappa9)/2508226560.0L + (kappa3_2*kappa5*kappa9)/3135283200.0L + (kappa3_2*kappa4*kappa10)/6270566400.0L + (kappa3_3*kappa11)/51732172800.0L; coeffs[21] += (kappa4_3*kappa5_2)/2388787200.0L + (kappa3*kappa4*kappa5_3)/1492992000.0L + (kappa4_4*kappa6)/5733089280.0L + (kappa3*kappa4_2*kappa5*kappa6)/597196800.0L + (kappa3_2*kappa5_2*kappa6)/1492992000.0L + (kappa3_2*kappa4*kappa6_2)/1791590400.0L + (kappa3*kappa4_3*kappa7)/2508226560.0L + (kappa3_2*kappa4*kappa5*kappa7)/1045094400.0L + (kappa3_3*kappa6*kappa7)/4702924800.0L + (kappa3_2*kappa4_2*kappa8)/3344302080.0L + (kappa3_3*kappa5*kappa8)/6270566400.0L + (kappa3_3*kappa4*kappa9)/11287019520.0L + (kappa3_4*kappa10)/112870195200.0L; coeffs[23] += kappa4_6/137594142720.0L + (kappa3*kappa4_4*kappa5)/5733089280.0L + (kappa3_2*kappa4_2*kappa5_2)/2388787200.0L + (kappa3_3*kappa5_3)/13436928000.0L + (kappa3_2*kappa4_3*kappa6)/4299816960.0L + (kappa3_3*kappa4*kappa5*kappa6)/2687385600.0L + (kappa3_4*kappa6_2)/32248627200.0L + (kappa3_3*kappa4_2*kappa7)/7524679680.0L + (kappa3_4*kappa5*kappa7)/18811699200.0L + (kappa3_4*kappa4*kappa8)/30098718720.0L + (kappa3_5*kappa9)/338610585600.0L; coeffs[25] += (kappa3_2*kappa4_5)/68797071360.0L + (kappa3_3*kappa4_3*kappa5)/12899450880.0L + (kappa3_4*kappa4*kappa5_2)/21499084800.0L + (kappa3_4*kappa4_2*kappa6)/25798901760.0L + (kappa3_5*kappa5*kappa6)/80621568000.0L + (kappa3_5*kappa4*kappa7)/112870195200.0L + (kappa3_6*kappa8)/1354442342400.0L; coeffs[27] += (kappa3_4*kappa4_4)/247669456896 + (kappa3_5*kappa4_2*kappa5)/128994508800.0L + (kappa3_6*kappa5_2)/967458816000.0L + (kappa3_6*kappa4*kappa6)/580475289600.0L + (kappa3_7*kappa7)/7110822297600.0L; coeffs[29] += (kappa3_6*kappa4_3)/2786281390080.0L + (kappa3_7*kappa4*kappa5)/4063327027200.0L + (kappa3_8*kappa6)/48759924326400.0L; coeffs[31] += (kappa3_8*kappa4_2)/78015878922240.0L + (kappa3_9*kappa5)/438839318937600.0L; coeffs[33] += (kappa3_10*kappa4)/5266071827251200.0L; coeffs[35] += kappa3_12/1042682221795737600.0L; if (order_ >= 13U) { *maxdeg = 38U; const double kappa15 = cumulants_[14]; const double kappa4_7 = kappa4_3*kappa4_4; const double kappa3_13 = kappa3_6*kappa3_7; const double kappa3_14 = kappa3_7*kappa3_7; const double kappa9_2 = kappa9*kappa9; coeffs[14] += kappa15/1307674368000.0L; coeffs[16] += (kappa8*kappa9)/14631321600.0L + (kappa7*kappa10)/18289152000.0L + (kappa6*kappa11)/28740096000.0L + (kappa5*kappa12)/57480192000.0L + (kappa4*kappa13)/149448499200.0L + (kappa3*kappa14)/523069747200.0L; coeffs[18] += (kappa6_2*kappa7)/5225472000.0L + (kappa5*kappa7_2)/6096384000.0L + (kappa5*kappa6*kappa8)/3483648000.0L + (kappa4*kappa7*kappa8)/4877107200.0L + (kappa3*kappa8_2)/19508428800.0L + (kappa5_2*kappa9)/10450944000.0L + (kappa4*kappa6*kappa9)/6270566400.0L + (kappa3*kappa7*kappa9)/10973491200.0L + (kappa4*kappa5*kappa10)/10450944000.0L + (kappa3*kappa6*kappa10)/15676416000.0L + (kappa4_2*kappa11)/45984153600.0L + (kappa3*kappa5*kappa11)/28740096000.0L + (kappa3*kappa4*kappa12)/68976230400.0L + (kappa3_2*kappa13)/448345497600.0L; coeffs[20] += (kappa5_3*kappa6)/7464960000.0L + (kappa4*kappa5*kappa6_2)/2985984000.0L + (kappa3*kappa6_3)/13436928000.0L + (kappa4*kappa5_2*kappa7)/3483648000.0L + (kappa4_2*kappa6*kappa7)/4180377600.0L + (kappa3*kappa5*kappa6*kappa7)/2612736000.0L + (kappa3*kappa4*kappa7_2)/7315660800.0L + (kappa4_2*kappa5*kappa8)/5573836800.0L + (kappa3*kappa5_2*kappa8)/6967296000.0L + (kappa3*kappa4*kappa6*kappa8)/4180377600.0L + (kappa3_2*kappa7*kappa8)/14631321600.0L + (kappa4_3*kappa9)/30098718720.0L + (kappa3*kappa4*kappa5*kappa9)/6270566400.0L + (kappa3_2*kappa6*kappa9)/18811699200.0L + (kappa3*kappa4_2*kappa10)/25082265600.0L + (kappa3_2*kappa5*kappa10)/31352832000.0L + (kappa3_2*kappa4*kappa11)/68976230400.0L + (kappa3_3*kappa12)/620786073600.0L; coeffs[22] += (kappa4_2*kappa5_3)/11943936000.0L + (kappa3*kappa5_4)/29859840000.0L + (kappa4_3*kappa5*kappa6)/7166361600.0L + (kappa3*kappa4*kappa5_2*kappa6)/2985984000.0L + (kappa3*kappa4_2*kappa6_2)/7166361600.0L + (kappa3_2*kappa5*kappa6_2)/8957952000.0L + (kappa4_4*kappa7)/40131624960.0L + (kappa3*kappa4_2*kappa5*kappa7)/4180377600.0L + (kappa3_2*kappa5_2*kappa7)/10450944000.0L + (kappa3_2*kappa4*kappa6*kappa7)/6270566400.0L + (kappa3_3*kappa7_2)/65840947200.0L + (kappa3*kappa4_3*kappa8)/20065812480.0L + (kappa3_2*kappa4*kappa5*kappa8)/8360755200.0L + (kappa3_3*kappa6*kappa8)/37623398400.0L + (kappa3_2*kappa4_2*kappa9)/30098718720.0L + (kappa3_3*kappa5*kappa9)/56435097600.0L + (kappa3_3*kappa4*kappa10)/112870195200.0L + (kappa3_4*kappa11)/1241572147200.0L; coeffs[24] += (kappa4_5*kappa5)/114661785600.0L + (kappa3*kappa4_3*kappa5_2)/14332723200.0L + (kappa3_2*kappa4*kappa5_3)/17915904000.0L + (kappa3*kappa4_4*kappa6)/34398535680.0L + (kappa3_2*kappa4_2*kappa5*kappa6)/7166361600.0L + (kappa3_3*kappa5_2*kappa6)/26873856000.0L + (kappa3_3*kappa4*kappa6_2)/32248627200.0L + (kappa3_2*kappa4_3*kappa7)/30098718720.0L + (kappa3_3*kappa4*kappa5*kappa7)/18811699200.0L + (kappa3_4*kappa6*kappa7)/112870195200.0L + (kappa3_3*kappa4_2*kappa8)/60197437440.0L + (kappa3_4*kappa5*kappa8)/150493593600.0L + (kappa3_4*kappa4*kappa9)/270888468480.0L + (kappa3_5*kappa10)/3386105856000.0L; coeffs[26] += (kappa3*kappa4_6)/825564856320.0L + (kappa3_2*kappa4_4*kappa5)/68797071360.0L + (kappa3_3*kappa4_2*kappa5_2)/42998169600.0L + (kappa3_4*kappa5_3)/322486272000.0L + (kappa3_3*kappa4_3*kappa6)/77396705280.0L + (kappa3_4*kappa4*kappa5*kappa6)/64497254400.0L + (kappa3_5*kappa6_2)/967458816000.0L + (kappa3_4*kappa4_2*kappa7)/180592312320.0L + (kappa3_5*kappa5*kappa7)/564350976000.0L + (kappa3_5*kappa4*kappa8)/902961561600.0L + (kappa3_6*kappa9)/12189981081600.0L; coeffs[28] += (kappa3_3*kappa4_5)/1238347284480.0L + (kappa3_4*kappa4_3*kappa5)/309586821120.0L + (kappa3_5*kappa4*kappa5_2)/644972544000.0L + (kappa3_5*kappa4_2*kappa6)/773967052800.0L + (kappa3_6*kappa5*kappa6)/2902376448000.0L + (kappa3_6*kappa4*kappa7)/4063327027200.0L + (kappa3_7*kappa8)/56886578380800.0L; coeffs[30] += (kappa3_5*kappa4_4)/7430083706880.0L + (kappa3_6*kappa4_2*kappa5)/4643802316800.0L + (kappa3_7*kappa5_2)/40633270272000.0L + (kappa3_7*kappa4*kappa6)/24379962163200.0L + (kappa3_8*kappa7)/341319470284800.0L; coeffs[32] += (kappa3_7*kappa4_3)/117023818383360.0L + (kappa3_8*kappa4*kappa5)/195039697305600.0L + (kappa3_9*kappa6)/2633035913625600.0L; coeffs[34] += (kappa3_9*kappa4_2)/4212857461800960.0L + (kappa3_10*kappa5)/26330359136256000.0L; coeffs[36] += (kappa3_11*kappa4)/347560740598579200.0L; coeffs[38] += kappa3_13/81329213300067532800.0L; if (order_ >= 14U) { *maxdeg = 41U; const double kappa16 = cumulants_[15]; coeffs[15] += kappa16/20922789888000.0L; coeffs[17] += kappa9_2/263363788800.0L + (kappa8*kappa10)/146313216000.0L + (kappa7*kappa11)/201180672000.0L + (kappa6*kappa12)/344881152000.0L + (kappa5*kappa13)/747242496000.0L + (kappa4*kappa14)/2092278988800.0L + (kappa3*kappa15)/7846046208000.0L; coeffs[19] += (kappa6*kappa7_2)/36578304000.0L + (kappa6_2*kappa8)/41803776000.0L + (kappa5*kappa7*kappa8)/24385536000.0L + (kappa4*kappa8_2)/78033715200.0L + (kappa5*kappa6*kappa9)/31352832000.0L + (kappa4*kappa7*kappa9)/43893964800.0L + (kappa3*kappa8*kappa9)/87787929600.0L + (kappa5_2*kappa10)/104509440000.0L + (kappa4*kappa6*kappa10)/62705664000.0L + (kappa3*kappa7*kappa10)/109734912000.0L + (kappa4*kappa5*kappa11)/114960384000.0L + (kappa3*kappa6*kappa11)/172440576000.0L + (kappa4_2*kappa12)/551809843200.0L + (kappa3*kappa5*kappa12)/344881152000.0L + (kappa3*kappa4*kappa13)/896690995200.0L + (kappa3_2*kappa14)/6276836966400.0L; coeffs[21] += (kappa5_2*kappa6_2)/29859840000.0L + (kappa4*kappa6_3)/53747712000.0L + (kappa5_3*kappa7)/52254720000.0L + (kappa4*kappa5*kappa6*kappa7)/10450944000.0L + (kappa3*kappa6_2*kappa7)/31352832000.0L + (kappa4_2*kappa7_2)/58525286400.0L + (kappa3*kappa5*kappa7_2)/36578304000.0L + (kappa4*kappa5_2*kappa8)/27869184000.0L + (kappa4_2*kappa6*kappa8)/33443020800.0L + (kappa3*kappa5*kappa6*kappa8)/20901888000.0L + (kappa3*kappa4*kappa7*kappa8)/29262643200.0L + (kappa3_2*kappa8_2)/234101145600.0L + (kappa4_2*kappa5*kappa9)/50164531200.0L + (kappa3*kappa5_2*kappa9)/62705664000.0L + (kappa3*kappa4*kappa6*kappa9)/37623398400.0L + (kappa3_2*kappa7*kappa9)/131681894400.0L + (kappa4_3*kappa10)/300987187200.0L + (kappa3*kappa4*kappa5*kappa10)/62705664000.0L + (kappa3_2*kappa6*kappa10)/188116992000.0L + (kappa3*kappa4_2*kappa11)/275904921600.0L + (kappa3_2*kappa5*kappa11)/344881152000.0L + (kappa3_2*kappa4*kappa12)/827714764800.0L + (kappa3_3*kappa13)/8070218956800.0L; coeffs[23] += (kappa4*kappa5_4)/119439360000.0L + (kappa4_2*kappa5_2*kappa6)/23887872000.0L + (kappa3*kappa5_3*kappa6)/44789760000.0L + (kappa4_3*kappa6_2)/85996339200.0L + (kappa3*kappa4*kappa5*kappa6_2)/17915904000.0L + (kappa3_2*kappa6_3)/161243136000.0L + (kappa4_3*kappa5*kappa7)/50164531200.0L + (kappa3*kappa4*kappa5_2*kappa7)/20901888000.0L + (kappa3*kappa4_2*kappa6*kappa7)/25082265600.0L + (kappa3_2*kappa5*kappa6*kappa7)/31352832000.0L + (kappa3_2*kappa4*kappa7_2)/87787929600.0L + (kappa4_4*kappa8)/321052999680.0L + (kappa3*kappa4_2*kappa5*kappa8)/33443020800.0L + (kappa3_2*kappa5_2*kappa8)/83607552000.0L + (kappa3_2*kappa4*kappa6*kappa8)/50164531200.0L + (kappa3_3*kappa7*kappa8)/263363788800.0L + (kappa3*kappa4_3*kappa9)/180592312320.0L + (kappa3_2*kappa4*kappa5*kappa9)/75246796800.0L + (kappa3_3*kappa6*kappa9)/338610585600.0L + (kappa3_2*kappa4_2*kappa10)/300987187200.0L + (kappa3_3*kappa5*kappa10)/564350976000.0L + (kappa3_3*kappa4*kappa11)/1241572147200.0L + (kappa3_4*kappa12)/14898865766400.0L; coeffs[25] += (kappa4_4*kappa5_2)/229323571200.0L + (kappa3*kappa4_2*kappa5_3)/71663616000.0L + (kappa3_2*kappa5_4)/358318080000.0L + (kappa4_5*kappa6)/687970713600.0L + (kappa3*kappa4_3*kappa5*kappa6)/42998169600.0L + (kappa3_2*kappa4*kappa5_2*kappa6)/35831808000.0L + (kappa3_2*kappa4_2*kappa6_2)/85996339200.0L + (kappa3_3*kappa5*kappa6_2)/161243136000.0L + (kappa3*kappa4_4*kappa7)/240789749760.0L + (kappa3_2*kappa4_2*kappa5*kappa7)/50164531200.0L + (kappa3_3*kappa5_2*kappa7)/188116992000.0L + (kappa3_3*kappa4*kappa6*kappa7)/112870195200.0L + (kappa3_4*kappa7_2)/1580182732800.0L + (kappa3_2*kappa4_3*kappa8)/240789749760.0L + (kappa3_3*kappa4*kappa5*kappa8)/150493593600.0L + (kappa3_4*kappa6*kappa8)/902961561600.0L + (kappa3_3*kappa4_2*kappa9)/541776936960.0L + (kappa3_4*kappa5*kappa9)/1354442342400.0L + (kappa3_4*kappa4*kappa10)/2708884684800.0L + (kappa3_5*kappa11)/37247164416000.0L; coeffs[27] += kappa4_7/23115815976960.0L + (kappa3*kappa4_5*kappa5)/687970713600.0L + (kappa3_2*kappa4_3*kappa5_2)/171992678400.0L + (kappa3_3*kappa4*kappa5_3)/322486272000.0L + (kappa3_2*kappa4_4*kappa6)/412782428160.0L + (kappa3_3*kappa4_2*kappa5*kappa6)/128994508800.0L + (kappa3_4*kappa5_2*kappa6)/644972544000.0L + (kappa3_4*kappa4*kappa6_2)/773967052800.0L + (kappa3_3*kappa4_3*kappa7)/541776936960.0L + (kappa3_4*kappa4*kappa5*kappa7)/451480780800.0L + (kappa3_5*kappa6*kappa7)/3386105856000.0L + (kappa3_4*kappa4_2*kappa8)/1444738498560.0L + (kappa3_5*kappa5*kappa8)/4514807808000.0L + (kappa3_5*kappa4*kappa9)/8126654054400.0L + (kappa3_6*kappa10)/121899810816000.0L; coeffs[29] += (kappa3_2*kappa4_6)/9906778275840.0L + (kappa3_3*kappa4_4*kappa5)/1238347284480.0L + (kappa3_4*kappa4_2*kappa5_2)/1031956070400.0L + (kappa3_5*kappa5_3)/9674588160000.0L + (kappa3_4*kappa4_3*kappa6)/1857520926720.0L + (kappa3_5*kappa4*kappa5*kappa6)/1934917632000.0L + (kappa3_6*kappa6_2)/34828517376000.0L + (kappa3_5*kappa4_2*kappa7)/5417769369600.0L + (kappa3_6*kappa5*kappa7)/20316635136000.0L + (kappa3_6*kappa4*kappa8)/32506616217600.0L + (kappa3_7*kappa9)/511979205427200.0L; coeffs[31] += (kappa3_4*kappa4_5)/29720334827520.0L + (kappa3_5*kappa4_3*kappa5)/9287604633600.0L + (kappa3_6*kappa4*kappa5_2)/23219011584000.0L + (kappa3_6*kappa4_2*kappa6)/27862813900800.0L + (kappa3_7*kappa5*kappa6)/121899810816000.0L + (kappa3_7*kappa4*kappa7)/170659735142400.0L + (kappa3_8*kappa8)/2730555762278400.0L; coeffs[33] += (kappa3_6*kappa4_4)/267483013447680.0L + (kappa3_7*kappa4_2*kappa5)/195039697305600.0L + (kappa3_8*kappa5_2)/1950396973056000.0L + (kappa3_8*kappa4*kappa6)/1170238183833600.0L + (kappa3_9*kappa7)/18431251395379200.0L; coeffs[35] += (kappa3_8*kappa4_3)/5617143282401280.0L + (kappa3_9*kappa4*kappa5)/10532143654502400.0L + (kappa3_10*kappa6)/157982154817536000.0L; coeffs[37] += (kappa3_10*kappa4_2)/252771447708057600.0L + (kappa3_11*kappa5)/1737803702992896000.0L; coeffs[39] += (kappa3_12*kappa4)/25024373323097702400.0L; coeffs[41] += kappa3_14/6831653917205672755200.0L; if (order_ > EDGEWORTH_MAXORDER) assert(0); } } } } } } } } } } } } } if (order_ > 1U) { const double sigma = sqrt(cumulants_[1]); if (sigma != 1.0) for (unsigned i=2; i<=*maxdeg; ++i) coeffs[i] /= powl(sigma, i); } } void EdgeworthSeries1D::hermiteCoeffsNormal(long double* coeffs, unsigned* maxdeg) const { if (m_ == EDGEWORTH_CLASSICAL) { hermiteCoeffsClassical(coeffs, maxdeg); return; } for (unsigned i=0; i 2U); const double k1 = cumulants_[0]; const double k2 = cumulants_[1]; const double k3 = cumulants_[2]; *maxdeg = 2U; assert(MAXCOEFFS > *maxdeg); coeffs[0] = k1; coeffs[2] = k3/6.0; if (order_ > 1U) { assert(nk > 3U); const double k4 = cumulants_[3]; const double k2m1 = k2 - 1.0; *maxdeg = 5; assert(MAXCOEFFS > *maxdeg); coeffs[1] += (k1*k1 + k2m1)/2.0; coeffs[3] += (4*k1*k3 + k4)/24.0; coeffs[5] += k3*k3/72.0; if (order_ > 2U) { assert(nk > 4U); const double k5 = cumulants_[4]; *maxdeg = 8; assert(MAXCOEFFS > *maxdeg); coeffs[2] += k1*(k1*k1 + 3*k2m1)/6; coeffs[4] += (10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120; coeffs[6] += (k3*(2*k1*k3 + k4))/144; coeffs[8] += k3*k3*k3/1296; if (order_ > 3U) { assert(nk > 5U); const double k6 = cumulants_[5]; const double k1sq = k1*k1; const double k3sq = k3*k3; *maxdeg = 11; assert(MAXCOEFFS > *maxdeg); coeffs[3] += (k1sq*k1sq + 6*k1sq*k2m1 + 3*k2m1*k2m1)/24; coeffs[5] += (20*k1*k1sq*k3 + 60*k1*k2m1*k3 + 15*k1sq*k4 + 15*k2m1*k4 + 6*k1*k5 + k6)/720; coeffs[7] += (40*k1sq*k3sq + 40*k2m1*k3sq + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760; coeffs[9] += (k3sq*(4*k1*k3 + 3*k4))/5184; coeffs[11] += k3sq*k3sq/31104; if (order_ > 4U) assert(0); } } } } double EdgeworthSeriesMomentFcn::operator()(const double& x) const { const double fact = series_.densityFactor(x); switch (degree_) { case 0: return fact; case 1: return x*fact; default: return pow(x - mean_, degree_)*fact; } } } Index: trunk/examples/C++/CmdLine.hh =================================================================== --- trunk/examples/C++/CmdLine.hh (revision 917) +++ trunk/examples/C++/CmdLine.hh (revision 918) @@ -1,391 +1,403 @@ //========================================================================= // CmdLine.hh // // Simple command line parser for the C++ "main" program. It provides // functionality of "getopt" and "getopt_long" with a convenient interface. // Typical usage is as follows: // // #include // // #include "CmdLine.hh" // // int main(int argc, char *argv[]) // { // CmdLine cmdline(argc, argv); // // int i = 0; // double d = 0; // bool b = false; // std::string requiredOption; // std::vector positionalArgs; // // try { // /* Arguments are short and long versions of the option name. */ // /* Long version can be omitted. If you want to use the long */ // /* version only, call the two-argument method with the short */ -// /* version set to NULL. */ +// /* version set to 0. */ // cmdline.option("-i", "--integer") >> i; // cmdline.option("-d") >> d; // // /* Options that must be present on the command line */ // cmdline.require("-r", "--required") >> requiredOption; // -// /* Switches that do not require additional arguments */ +// /* Switches that do not require subsequent arguments */ // b = cmdline.has("-b"); // // /* Declare the end of option processing. Unconsumed options */ // /* will cause an exception to be thrown. */ // cmdline.optend(); // // /* Process all remaining arguments */ // while (cmdline) // { // std::string s; // cmdline >> s; // positionalArgs.push_back(s); // } // } // catch (const CmdLineError& e) { // std::cerr << "Error in " << cmdline.progname() << ": " // << e.str() << std::endl; // return 1; // } // // /* ..... Do main processing here ..... */ // // return 0; // } // // Short version options must use a single character. It is possible // to combine several short options together on the command line, -// for exampe, "-xzvf" is equivalent to "-x -z -v -f". +// for example, "-xzvf" is equivalent to "-x -z -v -f". // // Use standalone "-" to indicate that the next argument is either // an option value or the program argument (but not a switch), even // if it starts with "-". This is useful if the option value has to // be set to a negative number. // // Use standalone "--" (not preceded by standalone "-") to indicate // the end of options. All remaining arguments will be treated as // program arguments, even if they start with "-". // // Note that each of the parsing methods of the CmdLine class ("option", // "has", and "require") is greedy. These methods will consume all // corresponding options or switches and will set the result to the last // option value seen. It is therefore impossible to provide a collection // of values by using an option more than once. This is done in order to // avoid difficulties with deciding what to do when multiple option values // were consumed by the user code only partially. // // After the "optend()" call, the "argc()" method of the CmdLine // class can be used to determine the number of remaining program // arguments. If the expected number of arguments is known in advance, // the simplest way to get the arguments out is like this (should be // inside the "try" block): // // if (cmdline.argc() != n_expected) // throw CmdLineError("wrong number of command line arguments"); // cmdline >> arg0 >> arg1 >> ...; // -// // I. Volobouev -// March 2013 +// June 2023 //========================================================================= #ifndef CMDLINE_HH_ #define CMDLINE_HH_ #include #include +#include #include #include #include #ifdef __GNUC__ #include #include #include #endif -#include "geners/CPP11_shared_ptr.hh" +#ifdef __GXX_EXPERIMENTAL_CXX0X__ +#include +#define CmdLine_shared_ptr std::shared_ptr +#else +#include +#define CmdLine_shared_ptr std::tr1::shared_ptr +#endif // Subsequent classes will throw exceptions of the following class class CmdLineError { public: - inline CmdLineError(const char* msg = 0) + inline CmdLineError() + : os_(new std::ostringstream()) {} + + inline CmdLineError(const char* msg) : os_(new std::ostringstream()) {if (msg) *os_ << msg;} + inline CmdLineError(const std::string& msg) + : os_(new std::ostringstream()) {if (!msg.empty()) *os_ << msg;} + template inline CmdLineError& operator<<(const T& obj) { *os_ << obj; return *this; } inline std::string str() const {return os_->str();} private: - CPP11_shared_ptr os_; + CmdLine_shared_ptr os_; }; template inline void OneShotExtract(std::istringstream& is, T& obj) {is >> obj;} template<> inline void OneShotExtract(std::istringstream& is, std::string& obj) {obj = is.str(); is.seekg(0, std::ios_base::end);} class OneShotIStream { public: inline OneShotIStream() : valid_(false), readout_(false) {} inline OneShotIStream(const std::string& s) : str_(s), valid_(true), readout_(false) {} inline operator void*() const { return valid_ && !readout_ ? (void*)this : (void*)0; } template inline bool operator>>(T& obj) { if (readout_) throw CmdLineError() << "can't reuse command line argument \"" << str_ << '"'; readout_ = true; if (valid_) { std::istringstream is(str_); OneShotExtract(is, obj); if (is.bad() || is.fail()) throw CmdLineError() << "failed to parse command line argument \"" << str_ << '"' #ifdef __GNUC__ << ", " << demangle(obj) << " expected" #endif ; if (is.peek() != EOF) throw CmdLineError() << "extra characters in command line argument \"" << str_ << '"' #ifdef __GNUC__ << ", " << demangle(obj) << " expected" #endif ; } return valid_; } inline bool isValid() const {return valid_;} private: std::string str_; bool valid_; bool readout_; #ifdef __GNUC__ template inline std::string demangle(T& obj) const { int status; const std::type_info &ti = typeid(obj); char *realname = abi::__cxa_demangle(ti.name(), 0, 0, &status); std::string s(realname); free(realname); return s; } #endif }; class CmdLine { // Argument codes (second member of the pair): // 0 -- possible option value (or program argument, not yet known which) // 1 -- short option switch // 2 -- long option switch // 3 -- program argument typedef std::pair Pair; typedef std::list Optlist; inline Optlist::iterator find(const char* shortOpt, const char* longOpt) { Optlist::iterator iend = args_.end(); for (Optlist::iterator it = args_.begin(); it != iend; ++it) { if (shortOpt && it->second == 1 && it->first == shortOpt) return it; if (longOpt && it->second == 2 && it->first == longOpt) return it; } return iend; } public: inline CmdLine(const unsigned argc, const char* const argv[]) : nprogargs_(0) { // Parse the program name const char *progname = std::strrchr(argv[0], '/'); if (progname) ++progname; else progname = argv[0]; // Take into account program name mangling by GNU autotools if (strncmp(progname, "lt-", 3) == 0) progname += 3; progname_ = progname; // Make a list of arguments noting on the way if this is // a short option, long option, or possible option argument bool previousIsOpt = false; bool nextIsArg = false; for (unsigned i=1; isecond == 0) it->second = 3; args_.erase(it0); } return found; } inline OneShotIStream option(const char* shortOpt, const char* longOpt=0) { OneShotIStream result; for (Optlist::iterator it = find(shortOpt, longOpt); it != args_.end(); it = find(shortOpt, longOpt)) { Optlist::iterator it0(it); if (++it != args_.end()) if (it->second == 0) { result = OneShotIStream(it->first); args_.erase(it0, ++it); --nprogargs_; continue; } throw CmdLineError() << "missing command line argument for option \"" << it0->first << '"'; } return result; } inline OneShotIStream require(const char* shortOpt, const char* longOpt=0) { const OneShotIStream& is(option(shortOpt, longOpt)); if (!is.isValid()) { const char empty[] = ""; const char* s = shortOpt ? shortOpt : (longOpt ? longOpt : empty); throw CmdLineError() << "required command line option \"" << s << "\" is missing"; } return is; } inline void optend() const { for (Optlist::const_iterator it = args_.begin(); it != args_.end(); ++it) if (it->second == 1 || it->second == 2) throw CmdLineError("invalid command line option \"") << it->first << '"'; } inline operator void*() const { return (void*)(static_cast(nprogargs_)); } inline unsigned argc() const { return nprogargs_; } template inline CmdLine& operator>>(T& obj) { if (!nprogargs_) throw CmdLineError("no more input available on the command line"); Optlist::iterator it = args_.begin(); for (; it != args_.end(); ++it) if (it->second == 0 || it->second == 3) break; OneShotIStream is(it->first); args_.erase(it); --nprogargs_; is >> obj; return *this; } private: CmdLine(); std::string progname_; Optlist args_; unsigned nprogargs_; }; #endif // CMDLINE_HH_ Index: trunk/tests/test_StatAccumulator.cc =================================================================== --- trunk/tests/test_StatAccumulator.cc (revision 917) +++ trunk/tests/test_StatAccumulator.cc (revision 918) @@ -1,347 +1,347 @@ #include #include #include "UnitTest++.h" #include "test_utils.hh" #include "TestAccumulator.hh" #include "npstat/stat/StatAccumulator.hh" #include "npstat/stat/SampleAccumulator.hh" #include "npstat/stat/BinSummary.hh" #include "npstat/stat/arrayStats.hh" using namespace std; using namespace npstat; namespace { TEST(arrayCumulants) { const long double eps = 1.0e-15L; const unsigned maxpoints = 33; SampleAccumulator acc; for (unsigned i=0; i acc; for (unsigned i=0; i acc; for (unsigned i=0; i acc; for (unsigned i=0; i acc; for (unsigned i=0; i acc; for (unsigned i=0; i