Index: trunk/npstat/stat/EdgeworthSeries1D.cc =================================================================== --- trunk/npstat/stat/EdgeworthSeries1D.cc (revision 592) +++ trunk/npstat/stat/EdgeworthSeries1D.cc (revision 593) @@ -1,356 +1,373 @@ #include #include "npstat/nm/MathUtils.hh" #include "npstat/nm/SpecialFunctions.hh" #include "npstat/nm/findRootUsingBisections.hh" #include "npstat/stat/EdgeworthSeries1D.hh" #include "npstat/stat/distributionReadError.hh" #include "geners/binaryIO.hh" #include "geners/vectorIO.hh" #define MAXCOEFFS 12 #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) { if (order_ > 4) 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"); } else cumulants_.clear(); } double EdgeworthSeries1D::mean() const { if (minNCumulants()) return cumulants_[0]; else return 0.0; } double EdgeworthSeries1D::stdev() 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; 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; } 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::hermiteCoeffsNormal(long double* coeffs, unsigned* maxdeg) const { for (unsigned i=0; i 2U); const double k1 = m_ == EDGEWORTH_SEVERINI ? cumulants_[0] : 0.0; const double k2 = m_ == EDGEWORTH_SEVERINI ? cumulants_[1] : 1.0; 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); } } if (m_ == EDGEWORTH_CLASSICAL) { const double sigma = sqrt(cumulants_[1]); for (unsigned i=2; i<=*maxdeg; ++i) coeffs[i] /= powl(sigma, i); } } } 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 - series_.mean(), degree_)*fact; } } } Index: trunk/npstat/stat/EdgeworthSeries1D.hh =================================================================== --- trunk/npstat/stat/EdgeworthSeries1D.hh (revision 592) +++ trunk/npstat/stat/EdgeworthSeries1D.hh (revision 593) @@ -1,118 +1,124 @@ #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 EdgeworthSeriesMethod.hh header 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 4. // // 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_;} double mean() const; double stdev() const; //@} double densityFactor(double x) const; double cdfFactor(double x) 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 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), degree_(degree) {} inline virtual ~EdgeworthSeriesMomentFcn() {} virtual double operator()(const double&) const; private: const EdgeworthSeries1D& series_; unsigned degree_; }; } #endif // NPSTAT_EDGEWORTHSERIES1D_HH_ Index: trunk/tests/test_Distributions1D.cc =================================================================== --- trunk/tests/test_Distributions1D.cc (revision 592) +++ trunk/tests/test_Distributions1D.cc (revision 593) @@ -1,1731 +1,1755 @@ #include #include #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/stat/Distributions1D.hh" #include "npstat/stat/GaussianMixture1D.hh" #include "npstat/stat/JohnsonCurves.hh" #include "npstat/stat/TruncatedDistribution1D.hh" #include "npstat/stat/LeftCensoredDistribution.hh" #include "npstat/stat/RightCensoredDistribution.hh" #include "npstat/stat/QuantileTable1D.hh" #include "npstat/stat/DistributionMix1D.hh" #include "npstat/stat/VerticallyInterpolatedDistribution1D.hh" #include "npstat/stat/UGaussConvolution1D.hh" #include "npstat/stat/RatioOfNormals.hh" #include "npstat/stat/InterpolatedDistro1D1P.hh" #include "npstat/stat/TransformedDistribution1D.hh" #include "npstat/stat/EdgeworthSeries1D.hh" #include "npstat/stat/distro1DStats.hh" #include "npstat/stat/IdentityTransform1D.hh" #include "npstat/stat/AsinhTransform1D.hh" #include "npstat/stat/LogRatioTransform1D.hh" #include "npstat/stat/LocationScaleTransform1.hh" #include "npstat/stat/Distribution1DReader.hh" #include "npstat/stat/DistributionTransform1DReader.hh" #include "npstat/stat/arrayStats.hh" #include "npstat/nm/EquidistantSequence.hh" #include "npstat/nm/LinearMapper1d.hh" #include "npstat/nm/ClassicalOrthoPolys1D.hh" #include "npstat/nm/MathUtils.hh" #include "npstat/nm/GaussHermiteQuadrature.hh" #include "LogLogQuadratic1D.hh" #include "geners/pseudoIO.hh" #define LSQR2 1.41421356237309504880169L #define SQRTWOPIL 2.50662827463100050241577L using namespace npstat; using namespace std; struct Two : public npstat::ConstValue1 { inline Two() : npstat::ConstValue1(2.0) {} }; gs_enable_pseudo_io(Two) struct Three : public npstat::ConstValue1 { inline Three() : npstat::ConstValue1(3.0) {} }; gs_enable_pseudo_io(Three) struct Xsq { inline double operator()(const double x) const {return x*x;} inline bool operator==(const Xsq&) const {return true;} }; gs_enable_pseudo_io(Xsq) struct Xp1 { inline double operator()(const double x) const {return x + 1.0;} inline bool operator==(const Xp1&) const {return true;} }; gs_enable_pseudo_io(Xp1) 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 { double he0(const double x) { return 1.0; } double he1(const double x) { return x; } double he2(const double x) { return x*x - 1.0; } double he3(const double x) { return x*(x*x - 3.0); } double he4(const long double x) { const long double xsq = x*x; return 3 - 6*xsq + xsq*xsq; } double he5(const long double x) { const long double xsq = x*x; return x*(xsq*xsq - 10*xsq + 15); } double he6(const long double x) { const long double xsq = x*x; return -15 + 45*xsq - 15*xsq*xsq + xsq*xsq*xsq; } double he7(const long double x) { const long double xsq = x*x; return x*(-105 + 105*xsq - 21*xsq*xsq + xsq*xsq*xsq); } double he8(const long double x) { const long double xsq = x*x; const long double x4 = xsq*xsq; return 105 - 420*xsq + 210*x4 - 28*x4*xsq + x4*x4; } double he9(const long double x) { const long double xsq = x*x; const long double x4 = xsq*xsq; return x*(945 - 1260*xsq + 378*x4 - 36*xsq*x4 + x4*x4); } double he10(const long double x) { const long double xsq = x*x; const long double x4 = xsq*xsq; const long double x6 = x4*xsq; const long double x8 = x6*xsq; const long double x10 = x8*xsq; return -945 + 4725*xsq - 3150*x4 + 630*x6 - 45*x8 + x10; } double he11(const long double x) { const long double xsq = x*x; const long double x4 = xsq*xsq; const long double x6 = x4*xsq; const long double x8 = x6*xsq; const long double x10 = x8*xsq; return x*(-10395 + 17325*xsq - 6930*x4 + 990*x6 - 55*x8 + x10); } double he12(const long double x) { const long double xsq = x*x; const long double x4 = xsq*xsq; const long double x6 = x4*xsq; const long double x8 = x6*xsq; const long double x10 = x8*xsq; const long double x12 = x10*xsq; return 10395 - 62370*xsq + 51975*x4 - 13860*x6 + 1485*x8 - 66*x10 + x12; } void get_cumulants_severini(const EdgeworthSeries1D& s, double cumulants[8]) { assert(s.method() == EDGEWORTH_SEVERINI); double m[8]; GaussHermiteQuadrature quad(128); for (unsigned i=0; i<8; ++i) { EdgeworthSeriesMomentFcn fcn(s, i); m[i] = quad.integrateProb(0.0, 1.0, fcn); } for (unsigned i=0; i<4; ++i) cumulants[i] = m[i]; cumulants[4] = m[4] - 3*m[2]*m[2]; cumulants[5] = m[5] - 10*m[2]*m[3]; cumulants[6] = m[6] - 15*m[2]*m[4] - 10*m[3]*m[3] + 30*m[2]*m[2]*m[2]; cumulants[7] = m[7] + 210*m[2]*m[2]*m[3] - 35*m[3]*m[4] - 21*m[2]*m[5]; } void get_cumulants_classical(const EdgeworthSeries1D& s, double cumulants[8]) { assert(s.method() == EDGEWORTH_CLASSICAL); double m[8]; GaussHermiteQuadrature quad(128); for (unsigned i=0; i<8; ++i) { EdgeworthSeriesMomentFcn fcn(s, i); m[i] = quad.integrateProb(s.mean(), s.stdev(), fcn); } for (unsigned i=0; i<4; ++i) cumulants[i] = m[i]; cumulants[4] = m[4] - 3*m[2]*m[2]; cumulants[5] = m[5] - 10*m[2]*m[3]; cumulants[6] = m[6] - 15*m[2]*m[4] - 10*m[3]*m[3] + 30*m[2]*m[2]*m[2]; cumulants[7] = m[7] + 210*m[2]*m[2]*m[3] - 35*m[3]*m[4] - 21*m[2]*m[5]; } void io_test(const AbsDistribution1D& d) { std::ostringstream os; CHECK(d.classId().write(os)); CHECK(d.write(os)); std::istringstream is(os.str()); gs::ClassId id(is, 1); AbsDistribution1D* phoenix = AbsDistribution1D::read(id, is); CHECK(*phoenix == d); delete phoenix; } void standard_test(const AbsDistribution1D& d, const double range, const double eps, const bool do_io = true) { const double cdf0 = d.cdf(0.0); for (unsigned i=0; i<100; ++i) { const double r = 2*range*(test_rng() - 0.5); const double cdfr = d.cdf(r); if (cdfr > 0.0 && cdfr < 1.0) { CHECK_CLOSE(cdfr, 1.0 - d.exceedance(r), eps); CHECK_CLOSE(r, d.quantile(cdfr), eps); const double integ = simpson_integral( DensityFunctor1D(d), 0.0, r); CHECK_CLOSE(integ, cdfr-cdf0, eps); } } if (do_io) io_test(d); } void invcdf_test(const AbsDistribution1D& d, const double eps) { for (unsigned i=0; i<1000; ++i) { const double r = test_rng(); const double x = d.quantile(r); const double cdf = d.cdf(x); CHECK_CLOSE(cdf, r, eps); CHECK_CLOSE(cdf, 1.0 - d.exceedance(x), eps); } } + void invexceedance_test(const EdgeworthSeries1D& d, const double eps) + { + for (unsigned i=0; i<1000; ++i) + { + const double r = test_rng(); + const double x = d.inverseExceedance(r); + const double exc = d.exceedance(x); + CHECK_CLOSE(exc, r, eps); + CHECK_CLOSE(exc, 1.0 - d.cdf(x), eps); + } + } + TEST(Distributions1D_Gauss) { Gauss1D g(0., 1.); double value = simpson_integral(DensityFunctor1D(g), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, 6, 1.e-7); CHECK_CLOSE(7.61985302416e-24, g.cdf(-10.0), 1.e-34); CHECK_CLOSE(7.61985302416e-24, g.exceedance(10.0), 1.e-34); CHECK_CLOSE(2.75362411861e-89, g.cdf(-20.0), 1.e-99); CHECK_CLOSE(2.75362411861e-89, g.exceedance(20.0), 1.e-99); } TEST(Distributions1D_Bifgauss_0) { const double lfrac = 0.3; const double lsigma = lfrac*2.0; const double rsigma = 2.0 - lsigma; BifurcatedGauss1D g(0.0, 1.0, lfrac, 2.0, 1.5); double value = simpson_integral(DensityFunctor1D(g), -lsigma*g.nSigmasLeft(), rsigma*g.nSigmasRight()); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, lsigma*g.nSigmasLeft(), 1.e-7); } TEST(Distributions1D_Bifgauss_1) { const double eps = 1.0e-12; BifurcatedGauss1D g1(0.0, 2.0, 0.5, 2.0, 2.0); TruncatedGauss1D g2(0.0, 2.0, 2.0); for (unsigned i=0; i<100; ++i) { const double r = test_rng(); const double q1 = g1.quantile(r); const double q2 = g2.quantile(r); CHECK_CLOSE(q1, q2, eps); CHECK_CLOSE(g1.density(q1), g2.density(q1), eps); CHECK_CLOSE(g1.cdf(q1), g2.cdf(q1), eps); CHECK_CLOSE(g1.exceedance(q1), g2.exceedance(q1), eps); } } TEST(Distributions1D_Bifgauss_2) { BifurcatedGauss1D g(0.0, 1.0, 0.0, 2.0, 2.0); double value = simpson_integral(DensityFunctor1D(g), 0.0, 4.0); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, 4.0, 1.e-7); } TEST(Distributions1D_Bifgauss_3) { BifurcatedGauss1D g(0.0, 1.0, 1.0, 2.0, 2.0); double value = simpson_integral(DensityFunctor1D(g), -4.0, 0.0); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, 4.0, 1.e-7); IdentityTransform1D it(5.0); TransformedDistribution1D td(it, g); value = simpson_integral(DensityFunctor1D(td), -4.0, 0.0); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(td, 4.0, 1.e-7); } TEST(Distributions1D_uniform) { Uniform1D g(0., 1.); double value = simpson_integral(DensityFunctor1D(g), -0.1, 1.1); CHECK_CLOSE(1.0, value, 1.e-3); standard_test(g, 1, 1.e-7); } TEST(Distributions1D_exp) { Exponential1D g(0., 1.); standard_test(g, 20.0, 1.e-7); g.setScale(3.); standard_test(g, 60.0, 1.e-7); } TEST(Distributions1D_Pareto) { Pareto1D p1(0., 1., 1.5); double value = simpson_integral(DensityFunctor1D(p1), 1, 6); CHECK_CLOSE(0.93195861825602283, value, 1.e-8); invcdf_test(p1, 1.e-10); io_test(p1); } TEST(Distributions1D_MirroredGauss) { { MirroredGauss1D g(0., 1., 0.3, 0.1); const double value = simpson_integral(DensityFunctor1D(g), 0, 1); CHECK_CLOSE(1.0, value, 1.e-12); invcdf_test(g, 1.e-10); io_test(g); for (unsigned i=0; i<100; ++i) { const double y = test_rng(); const double cdfy = simpson_integral(DensityFunctor1D(g), 0, y); CHECK_CLOSE(cdfy, g.cdf(y), 1.0e-10); } } { MirroredGauss1D g(0., 2., 0.7, 0.1); const double value = simpson_integral(DensityFunctor1D(g), 0, 2); CHECK_CLOSE(1.0, value, 1.e-12); invcdf_test(g, 1.e-10); io_test(g); for (unsigned i=0; i<100; ++i) { const double y = 2.0*test_rng(); const double cdfy = simpson_integral(DensityFunctor1D(g), 0, y); CHECK_CLOSE(cdfy, g.cdf(y), 1.0e-10); } } } TEST(Distributions1D_Huber) { Huber1D h1(0., 1., 0.0); double value = simpson_integral(DensityFunctor1D(h1), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(h1, 6, 1.e-7); double tailW = 0.2; double sigma = 0.8; Huber1D h2(0., sigma, tailW); CHECK_EQUAL(tailW, h2.tailWeight()); value = simpson_integral(DensityFunctor1D(h2), -20, 20, 10000); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(h2, 6, 1.e-7); const double a = sigma*h2.tailStart(); value = simpson_integral(DensityFunctor1D(h2), -a, a, 10000); CHECK_CLOSE(1.0-tailW, value, 1.e-8); for (unsigned i=0; i<1000; ++i) { const double rndx = 6*a*(test_rng() - 0.5); const double cdf = h2.cdf(rndx); value = simpson_integral(DensityFunctor1D(h2), -20, rndx, 10000); CHECK_CLOSE(cdf, value, 1.e-8); CHECK_CLOSE(h2.quantile(cdf), rndx, 1.e-8); } } TEST(Distributions1D_Symbeta) { SymmetricBeta1D sb(0.0, 1.0, 3); double value = simpson_integral(DensityFunctor1D(sb), -1, 1); CHECK_CLOSE(1.0, value, 1.e-8); for (unsigned i=0; i<10; ++i) { SymmetricBeta1D sb(0.0, 2.0, i); standard_test(sb, 1.8, 1.e-8); } } TEST(Distributions1D_Moyal) { Moyal1D mo(0.0, 1.0); double value = simpson_integral(DensityFunctor1D(mo), -7, 1000, 100000); CHECK_CLOSE(1.0, value, 1.e-9); CHECK_CLOSE(0.24197072451914335, mo.density(0.0), 1.0e-14); CHECK_CLOSE(0.20131624406488799, mo.density(1.0), 1.0e-14); CHECK_CLOSE(0.31731050786291410, mo.cdf(0.0), 1.0e-12); CHECK_CLOSE(0.54416242936230305, mo.cdf(1.0), 1.0e-12); invcdf_test(mo, 1.e-8); io_test(mo); } TEST(Distributions1D_StudentsT1D) { for (unsigned i=1; i<10; ++i) { StudentsT1D t(0.0, 1.0, i); invcdf_test(t, 1.e-8); } } TEST(Distributions1D_beta) { Beta1D b(0.0, 1.0, 3.0, 4.0); double value = simpson_integral(DensityFunctor1D(b), 0, 1); CHECK_CLOSE(1.0, value, 1.e-8); for (unsigned i=1; i<5; ++i) for (unsigned j=1; j<5; ++j) { Beta1D sb(0.5, 2.0, i, j); standard_test(sb, 0.9, 1.e-3); } } TEST(Distributions1D_gamma) { Gamma1D g1(0., 0.8, 1.0); standard_test(g1, 20.0, 5.e-3); Gamma1D g2(0., 0.9, 2.0); standard_test(g2, 20.0, 2.e-5); Gamma1D g3(0., 1.1, 3.0); standard_test(g3, 20.0, 2.e-5); } TEST(UGaussConvolution1D) { UGaussConvolution1D ug(0, 1, -1, 1); standard_test(ug, 6.0, 1.0e-6); } TEST(Distributions1D_Transformed) { StaticDistributionTransform1DReader::registerClass< LocationScaleTransform1 >(); const double teps = 1.0e-12; Two two; Three three; LocationScaleTransform1 tr(two, three); Gauss1D g(0.0, 1.0); TransformedDistribution1D d1(tr, g); io_test(d1); Gauss1D g1(2.0, 3.0); for (unsigned i=0; i<100; ++i) { const double x = 10.0*(test_rng() - 0.5); CHECK_CLOSE(g1.density(x), d1.density(x), teps); CHECK_CLOSE(g1.cdf(x), d1.cdf(x), teps); CHECK_CLOSE(g1.exceedance(x), d1.exceedance(x), teps); const double y = test_rng(); CHECK_CLOSE(g1.quantile(y), d1.quantile(y), teps); } } TEST(Distributions1D_JohnsonSu) { const double teps = 1.0e-12; JohnsonSu jsu1(0.0, 1.0, 1.0, 10.0); standard_test(jsu1, 6, 1.e-8); AsinhTransform1D tr1(jsu1.getDelta(), jsu1.getLambda(), jsu1.getGamma(), jsu1.getXi()); Gauss1D g(0.0, 1.0); TransformedDistribution1D d1(tr1, g); standard_test(d1, 6, 1.e-8); for (unsigned i=0; i<100; ++i) { const double x = 10.0*(test_rng() - 0.5); CHECK_CLOSE(jsu1.density(x), d1.density(x), teps); CHECK_CLOSE(jsu1.cdf(x), d1.cdf(x), teps); CHECK_CLOSE(jsu1.exceedance(x), d1.exceedance(x), teps); const double y = test_rng(); CHECK_CLOSE(jsu1.quantile(y), d1.quantile(y), teps); } JohnsonSu jsu2(0.1, 1.0, -1.0, 10.0); standard_test(jsu2, 6, 1.e-8); AsinhTransform1D tr2(jsu2.getDelta(), jsu2.getLambda(), jsu2.getGamma(), jsu2.getXi()); TransformedDistribution1D d2(tr2, g); standard_test(d2, 6, 1.e-8); double mean, sigma, skew, kurt; distro1DStats(jsu2, -30.0, 20.0, 500000, &mean, &sigma, &skew, &kurt); CHECK_CLOSE(0.1, mean, 1.0e-5); CHECK_CLOSE(1.0, sigma, 1.0e-4); CHECK_CLOSE(-1.0, skew, 0.01); CHECK_CLOSE(10.0, kurt, 0.2); } TEST(Distributions1D_JohnsonSb) { const double teps = 1.0e-12; JohnsonSb jsb1(0.0, 1.0, 1.0, 2.9); LogRatioTransform1D tr1(jsb1.getDelta(), jsb1.getLambda(), jsb1.getGamma(), jsb1.getXi()); Gauss1D g(0.0, 1.0); TransformedDistribution1D d1(tr1, g); const double lolim = jsb1.getXi(); const double range = jsb1.getLambda(); for (unsigned i=0; i<100; ++i) { const double x = lolim + range*test_rng(); CHECK_CLOSE(jsb1.density(x), d1.density(x), teps); CHECK_CLOSE(jsb1.cdf(x), d1.cdf(x), teps); CHECK_CLOSE(jsb1.exceedance(x), d1.exceedance(x), teps); const double y = test_rng(); CHECK_CLOSE(jsb1.quantile(y), d1.quantile(y), teps); } } TEST(Distributions1D_LogNormal0) { LogNormal g(0., 1., 0.); double value = simpson_integral(DensityFunctor1D(g), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, 6, 1.e-6); } TEST(Distributions1D_LogNormal1) { LogNormal ln1(0., 1., 2.0); double value = simpson_integral(DensityFunctor1D(ln1), -6, 20); CHECK_CLOSE(1.0, value, 1.e-6); standard_test(ln1, 7, 2.e-4); } TEST(Distributions1D_LogNormal2) { LogNormal ln2(0., 1., -2.0); double value = simpson_integral(DensityFunctor1D(ln2), -20, 6); CHECK_CLOSE(1.0, value, 1.e-6); standard_test(ln2, 7, 4.e-4); } TEST(Distributions1D_LogNormal3) { LogNormal ln3(-0.5, 1., 2.0); standard_test(ln3, 7, 1.e-3); } TEST(Distributions1D_LogNormal4) { LogNormal ln4(-0.5, 1., -2.0); standard_test(ln4, 7, 1.e-3); } TEST(Distributions1D_LogNormal5) { LogNormal ln5(0.5, 2., 2.0); standard_test(ln5, 7, 1.e-3); } TEST(Distributions1D_LogNormal6) { LogNormal ln6(0.5, 2., -2.0); standard_test(ln6, 7, 1.e-3); } TEST(Distributions1D_LogNormal7) { LogNormal ln7(1., 2., 10.0); standard_test(ln7, 7, 1.e-3); } TEST(Distributions1D_LogNormal8) { LogNormal ln8(1., 2., -10.0); standard_test(ln8, 7, 1.e-3); } TEST(Distributions1D_misc) { Cauchy1D ch(0.0, 0.5); standard_test(ch, 5, 1.e-8); TruncatedGauss1D tg(0.0, 2.0, 3.0); standard_test(tg, 6.0, 1.e-8); Gauss1D gauss(0.0, 2.0); TruncatedDistribution1D td(gauss, -6.0, 6.0); standard_test(td, 6.0, 1.e-8); CHECK_CLOSE(td.density(0.5), tg.density(0.5), 1.0e-14); CHECK_CLOSE(td.cdf(0.5), tg.cdf(0.5), 1.0e-14); CHECK_CLOSE(td.quantile(0.3), tg.quantile(0.3), 1.0e-14); CHECK_CLOSE(td.exceedance(0.3), tg.exceedance(0.3), 1.0e-14); } TEST(EdgeworthSeries1D_io) { const unsigned order = 2; std::vector cumulants(4); for (unsigned i=0; i dummy; double empCum[8]; EdgeworthSeries1D es1(dummy, EDGEWORTH_SEVERINI, order, false); CHECK_CLOSE(0.241970724519143, es1.density(1.0), 1.e-14); CHECK_CLOSE(7.61985302416e-24, es1.cdf(-10.0), 1.e-34); CHECK_CLOSE(7.61985302416e-24, es1.exceedance(10.0), 1.e-34); standard_test(es1, 6, 1.e-8); + invexceedance_test(es1, 1.e-8); get_cumulants_severini(es1, empCum); CHECK_CLOSE(1.0, empCum[0], 1.e-14); CHECK_CLOSE(0.0, empCum[1], 1.e-14); CHECK_CLOSE(1.0, empCum[2], 1.e-14); CHECK_CLOSE(0.0, empCum[3], 1.e-14); CHECK_CLOSE(0.0, empCum[4], 1.e-14); EdgeworthSeries1D es2(dummy, EDGEWORTH_SEVERINI, order, true); CHECK_CLOSE(0.241970724519143, es2.density(1.0), 1.e-14); CHECK_CLOSE(7.61985302416e-24, es2.cdf(-10.0), 1.e-34); CHECK_CLOSE(7.61985302416e-24, es2.exceedance(10.0), 1.e-34); standard_test(es2, 6, 1.e-8); + invexceedance_test(es2, 1.e-8); get_cumulants_severini(es2, empCum); CHECK_CLOSE(1.0, empCum[0], 1.e-14); CHECK_CLOSE(0.0, empCum[1], 1.e-14); CHECK_CLOSE(1.0, empCum[2], 1.e-14); CHECK_CLOSE(0.0, empCum[3], 1.e-14); CHECK_CLOSE(0.0, empCum[4], 1.e-14); EdgeworthSeries1D es3(dummy, EDGEWORTH_CLASSICAL, order, false); CHECK_CLOSE(0.241970724519143, es3.density(1.0), 1.e-14); CHECK_CLOSE(7.61985302416e-24, es3.cdf(-10.0), 1.e-34); CHECK_CLOSE(7.61985302416e-24, es3.exceedance(10.0), 1.e-34); standard_test(es3, 6, 1.e-8); + invexceedance_test(es3, 1.e-8); get_cumulants_classical(es3, empCum); CHECK_CLOSE(1.0, empCum[0], 1.e-14); CHECK_CLOSE(0.0, empCum[1], 1.e-14); CHECK_CLOSE(1.0, empCum[2], 1.e-14); CHECK_CLOSE(0.0, empCum[3], 1.e-14); CHECK_CLOSE(0.0, empCum[4], 1.e-14); EdgeworthSeries1D es4(dummy, EDGEWORTH_CLASSICAL, order, true); CHECK_CLOSE(0.241970724519143, es4.density(1.0), 1.e-14); CHECK_CLOSE(7.61985302416e-24, es4.cdf(-10.0), 1.e-34); CHECK_CLOSE(7.61985302416e-24, es4.exceedance(10.0), 1.e-34); standard_test(es4, 6, 1.e-8); + invexceedance_test(es4, 1.e-8); get_cumulants_classical(es4, empCum); CHECK_CLOSE(1.0, empCum[0], 1.e-14); CHECK_CLOSE(0.0, empCum[1], 1.e-14); CHECK_CLOSE(1.0, empCum[2], 1.e-14); CHECK_CLOSE(0.0, empCum[3], 1.e-14); CHECK_CLOSE(0.0, empCum[4], 1.e-14); } TEST(EdgeworthSeries1D_1) { const double range = 5.0; const unsigned ntries = 10; const unsigned order = 1; const double k1 = 3, k2 = 2, k3 = 0.5; double empCum[8]; std::vector c_slr(1); c_slr[0] = k1; std::vector c_gen(3); c_gen[0] = k1; c_gen[1] = k2; c_gen[2] = k3; EdgeworthSeries1D es1(c_slr, EDGEWORTH_SEVERINI, order, true); standard_test(es1, 6, 1.e-8); + invexceedance_test(es1, 1.e-8); for (unsigned i=0; i c_slr(2); c_slr[0] = k1; c_slr[1] = k2; std::vector c_gen(4); c_gen[0] = k1; c_gen[1] = k2; c_gen[2] = k3; c_gen[3] = k4; EdgeworthSeries1D es1(c_slr, EDGEWORTH_SEVERINI, order, true); standard_test(es1, 6, 1.e-8); for (unsigned i=0; i c_slr(3); c_slr[0] = k1; c_slr[1] = k2; c_slr[2] = k3; std::vector c_gen(5); c_gen[0] = k1; c_gen[1] = k2; c_gen[2] = k3; c_gen[3] = k4; c_gen[4] = k5; EdgeworthSeries1D es1(c_slr, EDGEWORTH_SEVERINI, order, true); - standard_test(es1, 6, 1.e-8); + standard_test(es1, 3, 1.e-8); for (unsigned i=0; i c_slr(4); c_slr[0] = k1; c_slr[1] = k2; c_slr[2] = k3; c_slr[3] = k4; std::vector c_gen(6); c_gen[0] = k1; c_gen[1] = k2; c_gen[2] = k3; c_gen[3] = k4; c_gen[4] = k5; c_gen[5] = k6; EdgeworthSeries1D es1(c_slr, EDGEWORTH_SEVERINI, order, true); - standard_test(es1, 6, 1.e-8); + standard_test(es1, 3, 1.e-8); + invexceedance_test(es1, 1.e-8); for (unsigned i=0; i quantiles(nquant); arrayQuantiles1D(data, sizeof(data)/sizeof(data[0]), 0.0, 1.0, &qs[0], &quantiles[0], nquant); BinnedDensity1D t1(location, width, data, sizeof(data)/sizeof(data[0]), 0); for (unsigned i=0; i(); LogLogQuadratic1D q1(1.0, 10.0, 0.5, 1.0); standard_test(q1, 1.0, 1.e-9); invcdf_test(q1, 1.e-14); } TEST(GaussianMixture1D_1) { std::vector entries; entries.push_back(GaussianMixtureEntry(3, 0, 1)); GaussianMixture1D g(0, 1, &entries[0], entries.size()); double value = simpson_integral(DensityFunctor1D(g), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(g, 6, 1.e-7); Gauss1D g2(0, 1); DistributionMix1D mix; mix.add(g2, 3.0); value = simpson_integral(DensityFunctor1D(mix), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(mix, 6, 1.e-7); VerticallyInterpolatedDistribution1D shmix(1U); shmix.add(g2, 3.0); value = simpson_integral(DensityFunctor1D(shmix), -6, 6); CHECK_CLOSE(1.0, value, 1.e-8); standard_test(shmix, 6, 1.e-7, false); for (unsigned i=0; i<100; ++i) { const double r = test_rng(); CHECK_CLOSE(shmix.quantile(r), mix.quantile(r), 1.e-15); } } TEST(GaussianMixture1D_2) { for (unsigned i=0; i<100; ++i) { Gauss1D g(test_rng(), 0.1+test_rng()); GaussianMixture1D mix(g); for (unsigned j=0; j<100; ++j) { const double x = -3 + 6*test_rng(); CHECK_CLOSE(g.density(x), mix.density(x), 1.0e-15); CHECK_CLOSE(g.cdf(x), mix.cdf(x), 1.0e-15); } } } TEST(DistributionMix1D_2) { for (unsigned i=0; i<100; ++i) { Gauss1D g(test_rng(), 0.1+test_rng()); DistributionMix1D mix; VerticallyInterpolatedDistribution1D shmix; mix.add(g, 1); shmix.add(g, 1); for (unsigned j=0; j<100; ++j) { const double x = -3 + 6*test_rng(); CHECK_CLOSE(g.density(x), mix.density(x), 1.0e-15); CHECK_CLOSE(g.cdf(x), mix.cdf(x), 1.0e-15); CHECK_CLOSE(g.density(x), shmix.density(x), 1.0e-15); CHECK_CLOSE(g.cdf(x), shmix.cdf(x), 1.0e-15); } } } TEST(GaussianMixture1D_3) { std::vector entries; entries.push_back(GaussianMixtureEntry(3, 0, 1)); entries.push_back(GaussianMixtureEntry(2, -0.5, 1.5)); entries.push_back(GaussianMixtureEntry(1, 0.5, 0.5)); GaussianMixture1D g(0, 1, &entries[0], entries.size()); double value = simpson_integral(DensityFunctor1D(g), -8, 8); CHECK_CLOSE(1.0, value, 1.e-7); standard_test(g, 4, 1.e-7); } TEST(DistributionMix1D_3) { DistributionMix1D mix; VerticallyInterpolatedDistribution1D shmix(3); Gauss1D g1(0, 1), g2(-0.5, 1.5), g3(0.5, 0.5); mix.add(g1, 3); mix.add(g2, 2); mix.add(g3, 1); shmix.add(g1, 3); shmix.add(g2, 2); shmix.add(g3, 1); double value = simpson_integral(DensityFunctor1D(mix), -8, 8); CHECK_CLOSE(1.0, value, 1.e-7); value = simpson_integral(DensityFunctor1D(shmix), -8, 8); CHECK_CLOSE(1.0, value, 1.e-7); standard_test(mix, 4, 1.e-7); standard_test(shmix, 4, 1.e-7, false); } TEST(GaussianMixture1D_4) { const double m0 = 1.3; const double s0 = 3/2.0; std::vector entries; entries.push_back(GaussianMixtureEntry(3, 0, 1)); entries.push_back(GaussianMixtureEntry(2, -0.5, 1.5)); entries.push_back(GaussianMixtureEntry(1, 0.5, 0.5)); GaussianMixture1D g(m0, s0, &entries[0], entries.size()); // G[x_, m_, s_] := 1/Sqrt[2 Pi]/s Exp[-(x - m)^2/s^2/2] // w1 = 3/6; m1 = 0; s1 = 1 // w2 = 2/6; m2 = -1/2; s2 = 3/2 // w3 = 1/6; m3 = 1/2; s3 = 1/2 // m0 = 1 + 3/10 // s0 = 3/2 // f = (w1*G[(x-m0)/s0,m1,s1]+w2*G[(x-m0)/s0,m2,s2]+w3*G[(x-m0)/s0,m3,s3])/s0 // mean = Integrate[x f, {x, -Infinity, Infinity}] // stdev = Sqrt[Integrate[(x - mean)^2 f, {x, -Infinity, Infinity}]] // q0 = Integrate[f, {x, -Infinity, 5/2}] // N[q0, 20] const double expmean = 47.0/40; const double expstdev = sqrt(203.0)/8.0; const double expcdf = 0.78400933536932041748; CHECK_CLOSE(expmean, g.mean(), 1.0e-15); CHECK_CLOSE(expstdev, g.stdev(), 1.0e-15); CHECK_CLOSE(expcdf, g.cdf(2.5), 1.0e-15); CHECK_CLOSE(2.5, g.quantile(expcdf), 1.0e-12); } TEST(LeftCensoredDistribution) { const double cutoff = 0.1234; const double eps = 1.0e-9; const double infty = -100.0; Gauss1D gauss(0.0, 1.0); TruncatedDistribution1D td(gauss, cutoff, 100.0); const double frac = gauss.exceedance(cutoff); LeftCensoredDistribution lc(td, frac, infty); CHECK_CLOSE(lc.cdf(cutoff-1.0), 1.0-frac, eps); CHECK_CLOSE(lc.exceedance(cutoff-1.0), frac, eps); CHECK_EQUAL(lc.cdf(infty-10.0), 0.0); CHECK_EQUAL(lc.exceedance(infty-10.0), 1.0); for (unsigned i=0; i<1000; ++i) { const double r = cutoff + test_rng()*5.0; CHECK_CLOSE(gauss.density(r), lc.density(r), eps); CHECK_CLOSE(gauss.exceedance(r), lc.exceedance(r), eps); CHECK_CLOSE(gauss.cdf(r), lc.cdf(r), eps); CHECK_CLOSE(lc.quantile(lc.cdf(r)), r, eps); } io_test(lc); } TEST(RightCensoredDistribution) { const double cutoff = 0.1234; const double eps = 1.0e-10; const double infty = 100.0; Gauss1D gauss(0.0, 1.0); TruncatedDistribution1D td(gauss, -100.0, cutoff); const double frac = gauss.cdf(cutoff); RightCensoredDistribution rc(td, frac, infty); CHECK_CLOSE(rc.cdf(cutoff+1.0), frac, eps); CHECK_CLOSE(rc.exceedance(cutoff+1.0), 1.0-frac, eps); CHECK_EQUAL(rc.cdf(infty+10.0), 1.0); CHECK_EQUAL(rc.exceedance(infty+10.0), 0.0); for (unsigned i=0; i<1000; ++i) { const double r = cutoff - test_rng()*5.0; CHECK_CLOSE(gauss.density(r), rc.density(r), eps); CHECK_CLOSE(gauss.exceedance(r), rc.exceedance(r), eps); CHECK_CLOSE(gauss.cdf(r), rc.cdf(r), eps); CHECK_CLOSE(rc.quantile(rc.cdf(r)), r, eps); } io_test(rc); } TEST(RatioOfNormals) { RatioOfNormals r1(-2.0, 1.0, 0.5, 1.0, 0.7); standard_test(r1, 10, 1.e-8); RatioOfNormals r2(1.0, 3.0, 2.0, 1.0, 0.5); standard_test(r2, 10, 1.e-7); const double x = 1.0; const double cdf = r2.cdf(x); const double quantile = r2.quantile(cdf); CHECK_CLOSE(x, quantile, 1.0e-8); } TEST(InterpolatedDistro1D1P) { const double eps = 1.0e-8; const double gridPoints[] = {1.0, 2.0}; double param = 1.3; const unsigned nGridPoints = sizeof(gridPoints)/sizeof(gridPoints[0]); GridAxis ax(std::vector(gridPoints, gridPoints+nGridPoints)); std::vector > distros; distros.push_back(CPP11_shared_ptr(new Gauss1D(1.0, 2.0))); distros.push_back(CPP11_shared_ptr(new Gauss1D(-1.0, 1.0))); InterpolatedDistro1D1P idist(ax, distros, param); standard_test(idist, 10.0, 1.e-7); LinearMapper1d meanmap(gridPoints[0], 1.0, gridPoints[1], -1.0); LinearMapper1d sigmap(gridPoints[0], 2.0, gridPoints[1], 1.0); for (unsigned i=0; i<100; ++i) { param = 1.0 + test_rng(); idist.setParameter(param); Gauss1D g(meanmap(param), sigmap(param)); const double r = test_rng(); const double x = g.quantile(r); CHECK_CLOSE(x, idist.quantile(r), eps); CHECK_CLOSE(g.density(x), idist.density(x), eps); CHECK_CLOSE(g.cdf(x), idist.cdf(x), eps); } } TEST(distro1DStats_1) { const double eps = 1.0e-7; for (unsigned i=0; i<100; ++i) { const double mean = test_rng(); const double sigma = 0.5 + test_rng(); double m, s, skew, kurt; distro1DStats(Gauss1D(mean, sigma), -10., 11., 100000, &m, &s, &skew, &kurt); CHECK_CLOSE(mean, m, eps); CHECK_CLOSE(sigma, s, eps); CHECK_CLOSE(0.0, skew, eps); CHECK_CLOSE(3.0, kurt, eps); } } }