Page MenuHomeHEPForge
Files F10881747

No OneTemporary
Actions

View Options

Index: trunk/tests/test_Distributions1D.cc
===================================================================
--- trunk/tests/test_Distributions1D.cc (revision 719)
+++ trunk/tests/test_Distributions1D.cc (revision 720)
@@ -1,2018 +1,2061 @@
#include <complex>
#include <iostream>
#include "UnitTest++.h"
#include "test_utils.hh"
#include "EdgeworthSeries1DOld.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/DeltaMixture1D.hh"
#include "npstat/stat/CompositeDistribution1D.hh"
#include "npstat/stat/ComparisonDistribution1D.hh"
#include "npstat/stat/distro1DStats.hh"
#include "npstat/stat/cumulantConversion.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 "npstat/nm/GaussLegendreQuadrature.hh"
#include "npstat/nm/SpecialFunctions.hh"
#include "npstat/rng/MersenneTwister.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<double,double>
{
inline Two() : npstat::ConstValue1<double,double>(2.0) {}
};
gs_enable_pseudo_io(Two)
struct Three : public npstat::ConstValue1<double,double>
{
inline Three() : npstat::ConstValue1<double,double>(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 cumulants_from_moments(const double m[10], double cumulants[10])
{
cumulants[0] = 0.0;
for (unsigned i=1; i<4; ++i)
cumulants[i] = m[i];
const double m2sq = m[2]*m[2];
const double m3sq = m[3]*m[3];
cumulants[4] = m[4] - 3*m2sq;
cumulants[5] = m[5] - 10*m[2]*m[3];
cumulants[6] = m[6] - 15*m[2]*m[4] - 10*m[3]*m[3] + 30*m2sq*m[2];
cumulants[7] = m[7] + 210*m2sq*m[3] - 35*m[3]*m[4] - 21*m[2]*m[5];
cumulants[8] = -630*m2sq*m2sq + 420*m2sq*m[4] - 35*m[4]*m[4] - 56*m[3]*m[5] +
28*m[2]*(20*m3sq - m[6]) + m[8];
cumulants[9] = -7560*m2sq*m[2]*m[3] + 560*m[3]*m3sq + 756*m2sq*m[5] - 126*m[4]*m[5] -
84*m[3]*m[6] + 36*m[2]*(70*m[3]*m[4] - m[7]) + m[9];
}
*/
void get_edgeworth_cumulants(const EdgeworthSeries1D& s, double cumulants[10])
{
double m[10];
for (unsigned i=0; i<10; ++i)
m[i] = s.empiricalCentralMoment(i);
if (s.order())
m[1] = s.cum(0);
npstat::convertCentralMomentsToCumulants(m, 9, cumulants);
}
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 standard_test_01(const AbsDistribution1D& d,
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 = test_rng();
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);
}
}
+ class MirroredGauss1D_meanFunctor : public Functor1<double,double>
+ {
+ public:
+ inline MirroredGauss1D_meanFunctor(const MirroredGauss1D& mg,
+ const double y)
+ : location_(mg.location()), scale_(mg.scale()),
+ sigmaOn0_1_(mg.sigmaOn0_1()), y_(y),
+ mapper_(location_, 0.0, location_+scale_, 1.0) {}
+
+ inline double operator()(const double& mean) const
+ {
+ MirroredGauss1D mg(location_, scale_, mapper_(mean), sigmaOn0_1_);
+ return mg.density(y_);
+ }
+
+ private:
+ double location_;
+ double scale_;
+ double sigmaOn0_1_;
+ double y_;
+ LinearMapper1d mapper_;
+ };
+
TEST(expectation)
{
const double eps = 1.0e-14;
Gauss1D g(2.0, 1.0);
auto lamx = [](double x) {return x;};
auto lamx2 = [](double x) {return x*x;};
auto lamx3 = [](double x) {return x*x*x;};
CHECK_CLOSE(2.0, g.expectation(lamx), eps);
CHECK_CLOSE(5.0, g.expectation(lamx2), eps);
CHECK_CLOSE(14.0, g.expectation(lamx3), eps);
}
TEST(cumulantConversion)
{
const double eps1 = 1.0e-15;
const double eps2 = 1.0e-13;
const unsigned maxOrder = 20U;
const unsigned npoints = 100U;
long double moments[maxOrder+1U];
long double moments2[maxOrder+1U];
long double cumulants[maxOrder+1U];
std::vector<double> points(npoints);
for (unsigned i=0; i<npoints; ++i)
points[i] = test_rng();
arrayCentralMoments(&points[0], points.size(), maxOrder, moments);
convertCentralMomentsToCumulants(moments, maxOrder, cumulants);
convertCumulantsToCentralMoments(cumulants, maxOrder, moments2);
for (unsigned k=0; k<=maxOrder; ++k)
{
if (k <= 12)
CHECK_CLOSE(moments[k], moments2[k], eps1);
else
CHECK_CLOSE(moments[k], moments2[k], eps2);
}
}
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_MirroredGauss_doublestochastic)
+ {
+ const double loc = 1.0;
+ const double scale = 3.0;
+ const double eps = 1.0e-12;
+
+ MirroredGauss1D g(loc, scale, 0.7, 0.1);
+ GaussLegendreQuadrature glq(1024);
+ CHECK_CLOSE(1.0, glq.integrate(DensityFunctor1D(g), loc, loc+scale), eps);
+
+ const unsigned ntries = 50;
+ for (unsigned i=0; i<ntries; ++i)
+ {
+ const double y = loc + test_rng()*scale;
+ MirroredGauss1D_meanFunctor mf(g, y);
+ CHECK_CLOSE(1.0, glq.integrate(mf, loc, loc+scale), eps);
+ }
+ }
+
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<Two,Three> >();
const double teps = 1.0e-12;
Two two;
Three three;
LocationScaleTransform1<Two,Three> 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<double> cumulants(4);
for (unsigned i=0; i<cumulants.size(); ++i)
cumulants[i] = i+1;
EdgeworthSeries1D s1(cumulants, EDGEWORTH_CLASSICAL, order, false);
io_test(s1);
}
TEST(EdgeworthSeries1D_0)
{
const unsigned order = 0;
std::vector<double> dummy;
double empCum[10];
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-7);
invexceedance_test(es1, 1.e-8);
get_edgeworth_cumulants(es1, empCum);
CHECK_CLOSE(0.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-7);
invexceedance_test(es2, 1.e-8);
get_edgeworth_cumulants(es2, empCum);
CHECK_CLOSE(0.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-7);
invexceedance_test(es3, 1.e-8);
get_edgeworth_cumulants(es3, empCum);
CHECK_CLOSE(0.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-7);
invexceedance_test(es4, 1.e-8);
get_edgeworth_cumulants(es4, empCum);
CHECK_CLOSE(0.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[10];
std::vector<double> c_slr(1);
c_slr[0] = k1;
std::vector<double> 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<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double fact = 1 + k1*he1(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*(k1);
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*(k1);
CHECK_CLOSE(k1, es1.cdfFactor(x), 1.0e-14);
CHECK_CLOSE(fact, es1.densityFactor(x), 1.0e-14);
CHECK_CLOSE(dens, es1.density(x), 1.0e-14);
CHECK_CLOSE(cdf, es1.cdf(x), 1.0e-14);
CHECK_CLOSE(exc, es1.exceedance(x), 1.0e-14);
}
get_edgeworth_cumulants(es1, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(1.0 - k1*k1, empCum[2], 1.e-14);
CHECK_CLOSE(2.0*k1*k1*k1, empCum[3], 1.e-14);
CHECK_CLOSE(-6*k1*k1*k1*k1, empCum[4], 1.e-14);
EdgeworthSeries1D es2(c_gen, EDGEWORTH_SEVERINI, order, false);
standard_test(es2, 6, 1.e-8);
invexceedance_test(es2, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double fact = 1 + k1*he1(x) + k3/6.0*he3(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdffact = k1*he0(x) + k3/6.0*he2(x);
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es2.cdfFactor(x), 1.0e-14);
CHECK_CLOSE(fact, es2.densityFactor(x), 1.0e-14);
CHECK_CLOSE(dens, es2.density(x), 1.0e-14);
CHECK_CLOSE(cdf, es2.cdf(x), 1.0e-14);
CHECK_CLOSE(exc, es2.exceedance(x), 1.0e-14);
}
get_edgeworth_cumulants(es2, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(1.0 - k1*k1, empCum[2], 1.e-14);
CHECK_CLOSE(2*k1*k1*k1 + k3, empCum[3], 1.e-14);
CHECK_CLOSE(-6*k1*k1*k1*k1 - 4*k1*k3, empCum[4], 1.e-14);
EdgeworthSeries1D es3(c_slr, EDGEWORTH_CLASSICAL, order, true);
standard_test(es3, 6, 1.e-8);
invexceedance_test(es3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = x0 - k1;
const double fact = 1;
const double expfact = exp(-x*x/2)/SQRTWOPIL;
const double dens = expfact*fact;
const double cdffact = 0.0;
const double cdf = ldgcdf(x) - expfact*cdffact;
const double exc = ldgexceedance(x) + expfact*cdffact;
CHECK_CLOSE(cdffact, es3.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es3.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es3.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es3.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es3.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es3, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, 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);
CHECK_CLOSE(0.0, empCum[5], 1.e-14);
CHECK_CLOSE(0.0, empCum[6], 1.e-14);
CHECK_CLOSE(0.0, empCum[7], 1.e-14);
EdgeworthSeries1D es4(c_gen, EDGEWORTH_CLASSICAL, order, false);
standard_test(es4, 6, 1.e-8);
invexceedance_test(es4, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = x0 - k1;
const double fact = 1 + k3/6.0*he3(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdffact = k3/6.0*he2(x);
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es4.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es4.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es4.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es4.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es4.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es4, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(1.0, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(0.0, empCum[4], 1.e-14);
CHECK_CLOSE(0.0, empCum[5], 1.e-14);
CHECK_CLOSE(-10*k3*k3, empCum[6], 1.e-14);
CHECK_CLOSE(0.0, empCum[7], 1.e-13);
}
TEST(EdgeworthSeries1D_2)
{
const double range = 5.0;
const unsigned ntries = 10;
const unsigned order = 2;
const double k1 = 0.3, k2 = 1.01, k3 = 0.2, k4 = 0.1;
const double k2m1 = k2 - 1.0;
const double sigma = sqrt(k2);
double empCum[10];
std::vector<double> c_slr(2);
c_slr[0] = k1;
c_slr[1] = k2;
std::vector<double> 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<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double cdffact = k1 + (k1*k1 + k2m1)/2.0*he1(x);
const double fact = 1 + k1*he1(x) + (k1*k1 + k2m1)/2.0*he2(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es1.cdfFactor(x), 1.0e-14);
CHECK_CLOSE(fact, es1.densityFactor(x), 1.0e-14);
CHECK_CLOSE(dens, es1.density(x), 1.0e-14);
CHECK_CLOSE(cdf, es1.cdf(x), 1.0e-14);
CHECK_CLOSE(exc, es1.exceedance(x), 1.0e-14);
}
get_edgeworth_cumulants(es1, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(-(k1*(k1*k1 + 3*k2m1)), empCum[3], 1.e-14);
CHECK_CLOSE(3*k1*k1*k1*k1 + 6*k1*k1*k2m1 - 3*k2m1*k2m1, empCum[4], 1.e-14);
EdgeworthSeries1D es2(c_gen, EDGEWORTH_SEVERINI, order, false);
standard_test(es2, 6, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double cdffact = (72*k1 + 36*(k1*k1 + k2m1)*he1(x) +
12*k3*he2(x) + 12*k1*k3*he3(x) +
3*k4*he3(x) + k3*k3*he5(x))/72;
const double fact = (72 + 72*k1*he1(x) + 36*(k1*k1 + k2m1)*he2(x) +
12*k3*he3(x) + 12*k1*k3*he4(x) +
3*k4*he4(x) + k3*k3*he6(x))/72;
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es2.cdfFactor(x), 1.0e-12);
CHECK_CLOSE(fact, es2.densityFactor(x), 1.0e-12);
CHECK_CLOSE(dens, es2.density(x), 1.0e-12);
CHECK_CLOSE(cdf, es2.cdf(x), 1.0e-12);
CHECK_CLOSE(exc, es2.exceedance(x), 1.0e-12);
}
get_edgeworth_cumulants(es2, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(-k1*k1*k1 - 3*k1*k2m1 + k3, empCum[3], 1.e-14);
CHECK_CLOSE(3*k1*k1*k1*k1 + 6*k1*k1*k2m1 - 3*k2m1*k2m1 + k4, empCum[4], 1.e-14);
EdgeworthSeries1D es3(c_slr, EDGEWORTH_CLASSICAL, order, true);
standard_test(es3, 6, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double fact = 1;
const double expfact = exp(-x*x/2)/SQRTWOPIL/sigma;
const double dens = expfact*fact;
const double cdffact = 0.0;
const double cdf = ldgcdf(x) - expfact*cdffact;
const double exc = ldgexceedance(x) + expfact*cdffact;
CHECK_CLOSE(cdffact, es3.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es3.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es3.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es3.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es3.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es3, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-12);
CHECK_CLOSE(k1, empCum[1], 1.e-12);
CHECK_CLOSE(k2, empCum[2], 1.e-12);
CHECK_CLOSE(0.0, empCum[3], 1.e-12);
CHECK_CLOSE(0.0, empCum[4], 1.e-12);
CHECK_CLOSE(0.0, empCum[5], 1.e-12);
CHECK_CLOSE(0.0, empCum[6], 1.e-12);
CHECK_CLOSE(0.0, empCum[7], 1.e-12);
EdgeworthSeries1D es4(c_gen, EDGEWORTH_CLASSICAL, order, false);
standard_test(es4, 6, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double cdffact = k3*he2(x)/6/sigma/sigma +
k4*he3(x)/24/pow(sigma,3) + k3*k3*he5(x)/72/pow(sigma,5);
const double fact = 1 + (k3*he3(x)/6/sigma/sigma +
k4*he4(x)/24/pow(sigma,3) + k3*k3*he6(x)/72/pow(sigma,5))/sigma;
const double dens = exp(-x*x/2)/SQRTWOPIL/sigma*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
CHECK_CLOSE(cdffact, es4.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es4.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es4.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es4.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es4.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es4, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(0.0, empCum[5], 1.e-14);
CHECK_CLOSE(0.0, empCum[6], 1.e-13);
CHECK_CLOSE(-35*k3*k4, empCum[7], 1.e-13);
}
TEST(EdgeworthSeries1D_3)
{
const double range = 5.0;
const unsigned ntries = 10;
const unsigned order = 3;
double k1 = 0.3, k2 = 1.01, k3 = 0.2, k4 = 0.1, k5 = 0.015;
double k2m1 = k2 - 1.0;
const double sigma = sqrt(k2);
double empCum[10];
std::vector<double> c_slr(3);
c_slr[0] = k1;
c_slr[1] = k2;
c_slr[2] = k3;
std::vector<double> 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, 3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double cdffact = k1 + (k1*k1*he1(x))/2 + (k2m1*he1(x))/2 + (k1*k1*k1*he2(x))/6 +
(k1*k2m1*he2(x))/2 + (k3*he2(x))/6;
const double fact = 1 + k1*he1(x) + (k1*k1*he2(x))/2 + (k2m1*he2(x))/2 + (k1*k1*k1*he3(x))/6 +
(k1*k2m1*he3(x))/2 + (k3*he3(x))/6;
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es1.cdfFactor(x), 1.0e-14);
CHECK_CLOSE(fact, es1.densityFactor(x), 1.0e-14);
CHECK_CLOSE(dens, es1.density(x), 1.0e-14);
CHECK_CLOSE(cdf, es1.cdf(x), 1.0e-14);
CHECK_CLOSE(exc, es1.exceedance(x), 1.0e-14);
}
get_edgeworth_cumulants(es1, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(-k1*k1*k1*k1 - 6*k1*k1*k2m1 - 3*k2m1*k2m1 - 4*k1*k3, empCum[4], 1.e-14);
EdgeworthSeries1D es2(c_gen, EDGEWORTH_SEVERINI, order, false);
standard_test(es2, 3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double cdffact = k1 + (k1*k1*he1(x))/2 + (k2m1*he1(x))/2 + (k1*k1*k1*he2(x))/6 +
(k1*k2m1*he2(x))/2 + (k3*he2(x))/6 + (k1*k3*he3(x))/6 +
(k4*he3(x))/24 + (k1*k1*k3*he4(x))/12 + (k2m1*k3*he4(x))/12 +
(k1*k4*he4(x))/24 + (k5*he4(x))/120 + (k3*k3*he5(x))/72 +
(k1*k3*k3*he6(x))/72 + (k3*k4*he6(x))/144 + (k3*k3*k3*he8(x))/1296;
const double fact = 1 + k1*he1(x) + (k1*k1*he2(x))/2 + (k2m1*he2(x))/2 + (k1*k1*k1*he3(x))/6 +
(k1*k2m1*he3(x))/2 + (k3*he3(x))/6 + (k1*k3*he4(x))/6 +
(k4*he4(x))/24 + (k1*k1*k3*he5(x))/12 + (k2m1*k3*he5(x))/12 +
(k1*k4*he5(x))/24 + (k5*he5(x))/120 + (k3*k3*he6(x))/72 +
(k1*k3*k3*he7(x))/72 + (k3*k4*he7(x))/144 + (k3*k3*k3*he9(x))/1296;
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es2.cdfFactor(x), 1.0e-12);
CHECK_CLOSE(fact, es2.densityFactor(x), 1.0e-12);
CHECK_CLOSE(dens, es2.density(x), 1.0e-12);
CHECK_CLOSE(cdf, es2.cdf(x), 1.0e-12);
CHECK_CLOSE(exc, es2.exceedance(x), 1.0e-12);
}
get_edgeworth_cumulants(es2, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(-k1*k1*k1*k1 - 6*k1*k1*k2m1 - 3*k2m1*k2m1 + k4, empCum[4], 1.e-14);
EdgeworthSeries1D es3(c_slr, EDGEWORTH_CLASSICAL, order, true);
standard_test(es3, 3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double cdffact = (k3/6)*he2(x)/sigma/sigma;
const double fact = 1 + (k3/6)*he3(x)/sigma/sigma/sigma;
const double expfact = exp(-x*x/2)/SQRTWOPIL/sigma;
const double dens = expfact*fact;
const double cdf = ldgcdf(x) - expfact*cdffact;
const double exc = ldgexceedance(x) + expfact*cdffact;
CHECK_CLOSE(cdffact, es3.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es3.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es3.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es3.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es3.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es3, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(0.0, empCum[4], 1.e-14);
CHECK_CLOSE(0.0, empCum[5], 1.e-14);
CHECK_CLOSE(-10*k3*k3, empCum[6], 1.e-14);
CHECK_CLOSE(0.0, empCum[7], 1.e-14);
EdgeworthSeries1D es4(c_gen, EDGEWORTH_CLASSICAL, order, false);
standard_test(es4, 3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double k1orig = k1;
const double k2m1orig = k2m1;
k1 = 0.0;
k2m1 = 0.0;
const double cdffact = (k1*k1*k1*he2(x)/powl(sigma,2))/6 +
(k1*k2m1*he2(x)/powl(sigma,2))/2 + (k3*he2(x)/powl(sigma,2))/6 + (k1*k3*he3(x)/powl(sigma,3))/6 +
(k4*he3(x)/powl(sigma,3))/24 + (k1*k1*k3*he4(x)/powl(sigma,4))/12 + (k2m1*k3*he4(x)/powl(sigma,4))/12 +
(k1*k4*he4(x)/powl(sigma,4))/24 + (k5*he4(x)/powl(sigma,4))/120 + (k3*k3*he5(x)/powl(sigma,5))/72 +
(k1*k3*k3*he6(x)/powl(sigma,6))/72 + (k3*k4*he6(x)/powl(sigma,6))/144 + (k3*k3*k3*he8(x)/powl(sigma,8))/1296;
const double fact = 1 + (k1*k1*k1*he3(x)/powl(sigma,3))/6 +
(k1*k2m1*he3(x)/powl(sigma,3))/2 + (k3*he3(x)/powl(sigma,3))/6 + (k1*k3*he4(x)/powl(sigma,4))/6 +
(k4*he4(x)/powl(sigma,4))/24 + (k1*k1*k3*he5(x)/powl(sigma,5))/12 + (k2m1*k3*he5(x)/powl(sigma,5))/12 +
(k1*k4*he5(x)/powl(sigma,5))/24 + (k5*he5(x)/powl(sigma,5))/120 + (k3*k3*he6(x)/powl(sigma,6))/72 +
(k1*k3*k3*he7(x)/powl(sigma,7))/72 + (k3*k4*he7(x)/powl(sigma,7))/144 + (k3*k3*k3*he9(x)/powl(sigma,9))/1296;
k1 = k1orig;
k2m1 = k2m1orig;
const double dens = exp(-x*x/2)/SQRTWOPIL/sigma*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
CHECK_CLOSE(cdffact, es4.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es4.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es4.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es4.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es4.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es4, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(k5, empCum[5], 1.e-14);
CHECK_CLOSE(0.0, empCum[6], 1.e-13);
CHECK_CLOSE(0.0, empCum[7], 1.e-13);
}
TEST(EdgeworthSeries1D_4)
{
const double range = 5;
const unsigned ntries = 10;
const unsigned order = 4;
double k1 = 0.3, k2 = 1.01, k3 = 0.2, k4 = 0.1, k5 = 0.015, k6 = 0.07;
double k2m1 = k2 - 1.0;
const double sigma = sqrt(k2);
double empCum[10];
std::vector<double> c_slr(4);
c_slr[0] = k1;
c_slr[1] = k2;
c_slr[2] = k3;
c_slr[3] = k4;
std::vector<double> 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, 3, 1.e-8);
invexceedance_test(es1, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x = range*(2.0*test_rng()-1.0);
const double cdffact = k1 + (k1*k1*he1(x))/2 + (k2m1*he1(x))/2 + (k1*k1*k1*he2(x))/6 +
(k1*k2m1*he2(x))/2 + (k3*he2(x))/6 +
(k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24*he3(x);
const double fact = 1 + k1*he1(x) + (k1*k1*he2(x))/2 + (k2m1*he2(x))/2 + (k1*k1*k1*he3(x))/6 +
(k1*k2m1*he3(x))/2 + (k3*he3(x))/6 +
(k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24*he4(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, es1.cdfFactor(x), 1.0e-14);
CHECK_CLOSE(fact, es1.densityFactor(x), 1.0e-14);
CHECK_CLOSE(dens, es1.density(x), 1.0e-14);
CHECK_CLOSE(cdf, es1.cdf(x), 1.0e-14);
CHECK_CLOSE(exc, es1.exceedance(x), 1.0e-14);
}
get_edgeworth_cumulants(es1, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(-pow(k1,5) - 10*k1*k1*k1*k2m1 - 15*k1*k2m1*k2m1 - 10*k1*k1*k3 - 10*k2m1*k3 - 5*k1*k4, empCum[5], 1.e-14);
EdgeworthSeries1D es2(c_gen, EDGEWORTH_SEVERINI, order, false);
standard_test(es2, 3, 1.e-8);
invexceedance_test(es2, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
double coeffs[13];
const double x = range*(2.0*test_rng()-1.0);
coeffs[0] = 1.0;
coeffs[1] = k1;
coeffs[2] = (k1*k1 + k2m1)/2;
coeffs[3] = (k1*k1*k1 + 3*k1*k2m1 + k3)/6;
coeffs[4] = (k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24;
coeffs[5] = (10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120;
coeffs[6] = (20*k1*k1*k1*k3 + 10*k3*k3 + 15*k1*k1*k4 + 15*k2m1*k4 + 6*k1*(10*k2m1*k3 + k5) + k6)/720;
coeffs[7] = (k3*(2*k1*k3 + k4))/144;
coeffs[8] = (40*k1*k1*k3*k3 + 40*k2m1*k3*k3 + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760;
coeffs[9] = k3*k3*k3/1296;
coeffs[10] = (k3*k3*(4*k1*k3 + 3*k4))/5184;
coeffs[11] = 0.0;
coeffs[12] = k3*k3*k3*k3/31104;
const double fact = hermiteSeriesSumProb(coeffs, 12, x);
const double cdffact = hermiteSeriesSumProb(coeffs+1, 11, x);
const double cdffactAlt = k1*he0(x) +
(k1*k1 + k2m1)/2*he1(x) +
(k1*k1*k1 + 3*k1*k2m1 + k3)/6*he2(x) +
(k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24*he3(x) +
(10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120*he4(x) +
(20*k1*k1*k1*k3 + 10*k3*k3 + 15*k1*k1*k4 + 15*k2m1*k4 + 6*k1*(10*k2m1*k3 + k5) + k6)/720*he5(x) +
(k3*(2*k1*k3 + k4))/144*he6(x) +
(40*k1*k1*k3*k3 + 40*k2m1*k3*k3 + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760*he7(x) +
k3*k3*k3/1296*he8(x) +
(k3*k3*(4*k1*k3 + 3*k4))/5184*he9(x) +
k3*k3*k3*k3/31104*he11(x);
const double factAlt = 1 + k1*he1(x) +
(k1*k1 + k2m1)/2*he2(x) +
(k1*k1*k1 + 3*k1*k2m1 + k3)/6*he3(x) +
(k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24*he4(x) +
(10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120*he5(x) +
(20*k1*k1*k1*k3 + 10*k3*k3 + 15*k1*k1*k4 + 15*k2m1*k4 + 6*k1*(10*k2m1*k3 + k5) + k6)/720*he6(x) +
(k3*(2*k1*k3 + k4))/144*he7(x) +
(40*k1*k1*k3*k3 + 40*k2m1*k3*k3 + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760*he8(x) +
k3*k3*k3/1296*he9(x) +
(k3*k3*(4*k1*k3 + 3*k4))/5184*he10(x) +
k3*k3*k3*k3/31104*he12(x);
const double dens = exp(-x*x/2)/SQRTWOPIL*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL*cdffact;
CHECK_CLOSE(cdffact, cdffactAlt, 1.0e-12);
CHECK_CLOSE(fact, factAlt, 1.0e-12);
CHECK_CLOSE(cdffact, es2.cdfFactor(x), 1.0e-12);
CHECK_CLOSE(fact, es2.densityFactor(x), 1.0e-12);
CHECK_CLOSE(dens, es2.density(x), 1.0e-12);
CHECK_CLOSE(cdf, es2.cdf(x), 1.0e-12);
CHECK_CLOSE(exc, es2.exceedance(x), 1.0e-12);
}
get_edgeworth_cumulants(es2, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(-pow(k1,5) - 10*k1*k1*k1*k2m1 - 15*k1*k2m1*k2m1 + k5, empCum[5], 1.e-14);
CHECK_CLOSE(5*pow(k1,6) + 45*pow(k1,4)*k2m1 + 45*k1*k1*k2m1*k2m1 -
15*k2m1*k2m1*k2m1 + k6, empCum[6], 1.e-14);
EdgeworthSeries1D es3(c_slr, EDGEWORTH_CLASSICAL, order, true);
standard_test(es3, 3, 1.e-8);
invexceedance_test(es3, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double cdffact = (k3*he2(x)/sigma/sigma)/6 +
(k4)/24*he3(x)/sigma/sigma/sigma;
const double fact = 1 + (k3*he3(x)/sigma/sigma/sigma)/6 +
(k4)/24*he4(x)/sigma/sigma/sigma/sigma;
const double expfact = exp(-x*x/2)/SQRTWOPIL/sigma;
const double dens = expfact*fact;
const double cdf = ldgcdf(x) - expfact*cdffact;
const double exc = ldgexceedance(x) + expfact*cdffact;
CHECK_CLOSE(cdffact, es3.cdfFactor(x0), 1.0e-14);
CHECK_CLOSE(fact, es3.densityFactor(x0), 1.0e-14);
CHECK_CLOSE(dens, es3.density(x0), 1.0e-14);
CHECK_CLOSE(cdf, es3.cdf(x0), 1.0e-14);
CHECK_CLOSE(exc, es3.exceedance(x0), 1.0e-14);
}
get_edgeworth_cumulants(es3, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(0.0, empCum[5], 1.e-14);
CHECK_CLOSE(-10*k3*k3, empCum[6], 1.e-14);
CHECK_CLOSE(-35*k3*k4, empCum[7], 1.e-14);
EdgeworthSeries1D es4(c_gen, EDGEWORTH_CLASSICAL, order, false);
standard_test(es4, 3, 1.e-8);
invexceedance_test(es4, 1.e-8);
for (unsigned i=0; i<ntries; ++i)
{
const double x0 = range*(2.0*test_rng()-1.0);
const double x = (x0 - k1)/sigma;
const double k1orig = k1;
const double k2m1orig = k2m1;
k1 = 0.0;
k2m1 = 0.0;
double coeffs[13];
coeffs[0] = 0.0;
coeffs[1] = k1;
coeffs[2] = (k1*k1 + k2m1)/2;
coeffs[3] = (k1*k1*k1 + 3*k1*k2m1 + k3)/6;
coeffs[4] = (k1*k1*k1*k1 + 6*k1*k1*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24;
coeffs[5] = (10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120;
coeffs[6] = (20*k1*k1*k1*k3 + 10*k3*k3 + 15*k1*k1*k4 + 15*k2m1*k4 + 6*k1*(10*k2m1*k3 + k5) + k6)/720;
coeffs[7] = (k3*(2*k1*k3 + k4))/144;
coeffs[8] = (40*k1*k1*k3*k3 + 40*k2m1*k3*k3 + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760;
coeffs[9] = k3*k3*k3/1296;
coeffs[10] = (k3*k3*(4*k1*k3 + 3*k4))/5184;
coeffs[11] = 0.0;
coeffs[12] = k3*k3*k3*k3/31104;
k1 = k1orig;
k2m1 = k2m1orig;
for (unsigned i=3; i<13; ++i)
coeffs[i] /= powl(sigma, i-1);
const double cdffact = hermiteSeriesSumProb(coeffs+1, 11, x);
const double fact = 1.0 + hermiteSeriesSumProb(coeffs, 12, x)/sigma;
const double dens = exp(-x*x/2)/SQRTWOPIL/sigma*fact;
const double cdf = ldgcdf(x) - exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
const double exc = ldgexceedance(x) + exp(-x*x/2)/SQRTWOPIL/sigma*cdffact;
CHECK_CLOSE(cdffact, es4.cdfFactor(x0), 1.0e-12);
CHECK_CLOSE(fact, es4.densityFactor(x0), 1.0e-12);
CHECK_CLOSE(dens, es4.density(x0), 1.0e-12);
CHECK_CLOSE(cdf, es4.cdf(x0), 1.0e-12);
CHECK_CLOSE(exc, es4.exceedance(x0), 1.0e-12);
}
get_edgeworth_cumulants(es4, empCum);
CHECK_CLOSE(0.0, empCum[0], 1.e-14);
CHECK_CLOSE(k1, empCum[1], 1.e-14);
CHECK_CLOSE(k2, empCum[2], 1.e-14);
CHECK_CLOSE(k3, empCum[3], 1.e-14);
CHECK_CLOSE(k4, empCum[4], 1.e-14);
CHECK_CLOSE(k5, empCum[5], 1.e-14);
CHECK_CLOSE(k6, empCum[6], 1.e-13);
CHECK_CLOSE(0.0, empCum[7], 1.e-13);
CHECK_CLOSE(0.0, empCum[8], 1.e-12);
CHECK_CLOSE(-126*k4*k5 - 84*k3*k6, empCum[9], 1.e-10);
}
TEST(EdgeworthSeries1D_5)
{
const double eps = 1.0e-12;
const unsigned ntries = 100;
const unsigned order = 4;
double k1 = 0.3, k2 = 1.01, k3 = 0.2, k4 = 0.1, k5 = 0.015, k6 = 0.07;
std::vector<double> 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_gen, EDGEWORTH_CLASSICAL, order, false);
EdgeworthSeries1DOld es2(c_gen, EDGEWORTH_CLASSICAL, order, false);
for (unsigned i=0; i<ntries; ++i)
{
const double x = -3.0 + 6.0*test_rng();
CHECK_CLOSE(es2.density(x), es1.density(x), eps);
}
}
TEST(Tabulated1D_1)
{
const double data[] = {0, 1, 2, 3, 3, 0};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<4; ++deg)
{
Tabulated1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-8);
}
}
TEST(Tabulated1D_2)
{
const double data[] = {1};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<4; ++deg)
{
Tabulated1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-8);
}
}
TEST(QuantileTable1D)
{
const double data[] = {0.2, 0.25, 0.4, 0.5, 0.55};
const double location = 1.0;
const double width = 2.5;
QuantileTable1D qt(location, width, data, sizeof(data)/sizeof(data[0]));
invcdf_test(qt, 1.e-12);
standard_test(qt, 3.5, 5.0e-3);
}
TEST(BinnedDensity1D_1)
{
const double data[] = {0, 1, 3, 3, 3, 0, 8};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<2; ++deg)
{
BinnedDensity1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-3);
}
}
TEST(BinnedDensity1D_2)
{
const double data[] = {0, 0, 1, 11, 3, 7, 0, 8, 0, 0};
const double location = 0.0;
const double width = 1.0;
const unsigned nquant = 1000;
EquidistantInLinearSpace qs(0.0, 1.0, nquant);
std::vector<double> 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<nquant; ++i)
CHECK_CLOSE(quantiles[i], t1.quantile(qs[i]), 1.0e-14);
}
TEST(BinnedDensity1D_3)
{
const double data[] = {1};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<2; ++deg)
{
BinnedDensity1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-3);
}
}
TEST(Quadratic1D_1)
{
Quadratic1D q1(-1.0, 2.0, 0.5/2.0, 0.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-12);
}
TEST(Quadratic1D_2)
{
Quadratic1D q1(-1.0, 2.0, 0.5/2.0, 0.2/2.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-12);
}
TEST(Quadratic1D_3)
{
Quadratic1D q1(10.0, 60.0, 0.16, 0.08);
CHECK_CLOSE(10.0, q1.quantile(0.0), 1.0e-14);
CHECK_CLOSE(70.0, q1.quantile(1.0), 1.0e-14);
CHECK_CLOSE(0.01391555555555556, q1.density(11.0), 1.0e-14);
CHECK_CLOSE(0.01633333333333333, q1.density(45.0), 1.0e-14);
}
TEST(LogQuadratic1D_1)
{
LogQuadratic1D q1(-1.0, 2.0, 0.5, 1.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-14);
}
TEST(LogQuadratic1D_2)
{
LogQuadratic1D q1(-1.0, 2.0, 0.5, -1.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-14);
}
TEST(LogQuadratic1D_3)
{
LogQuadratic1D q1(-1.0, 2.0, 0.5, 0.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-14);
}
TEST(LogQuadratic1D_4)
{
LogQuadratic1D q1(-1.0, 2.0, -0.5, 0.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-14);
}
TEST(LogQuadratic1D_5)
{
LogQuadratic1D q1(-1.0, 2.0, 0.0, -10.0/2.0);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-14);
}
TEST(LogQuadratic1D_6)
{
LogQuadratic1D q1(-1.0, 2.0, 0.0, 10.0/2.0);
standard_test(q1, 1.0, 1.e-9);
invcdf_test(q1, 1.e-14);
}
TEST(LogLogQuadratic1D_1)
{
StaticDistribution1DReader::registerClass<LogLogQuadratic1D>();
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<GaussianMixtureEntry> 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<GaussianMixtureEntry> 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<GaussianMixtureEntry> 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<double>(gridPoints, gridPoints+nGridPoints));
std::vector<CPP11_shared_ptr<const AbsDistribution1D> > distros;
distros.push_back(CPP11_shared_ptr<const AbsDistribution1D>(new Gauss1D(1.0, 2.0)));
distros.push_back(CPP11_shared_ptr<const AbsDistribution1D>(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);
}
}
TEST(DeltaMixture1D_1)
{
const double eps = 1.0e-12;
const unsigned nPoints = 3;
double coords[nPoints] = {0., 1., 2.};
double weights[nPoints] = {1., 3., 2.};
DeltaMixture1D dmix(5.0, 1.0, coords, weights, nPoints);
io_test(dmix);
CHECK_CLOSE(0.6666666666666666, dmix.cdf(6.0), eps);
CHECK_CLOSE(0.3333333333333333, dmix.exceedance(6.0), eps);
MersenneTwister rng;
int counts[nPoints] = {0};
for (unsigned i=0; i<6000; ++i)
{
double rnd;
dmix.random(rng, &rnd);
CHECK(rnd == 5.0 || rnd == 6.0 || rnd == 7.0);
if (rnd == 5.0)
++counts[0];
if (rnd == 6.0)
++counts[1];
if (rnd == 7.0)
++counts[2];
}
CHECK(std::abs(counts[0] - 1000) < 200);
CHECK(std::abs(counts[2] - 2000) < 300);
CHECK_CLOSE(0.5, dmix.integral(6.0, true, 6.0, true), eps);
CHECK_CLOSE(0.0, dmix.integral(6.0, false, 6.0, false), eps);
CHECK_CLOSE(0.0, dmix.integral(6.0, false, 6.0, true), eps);
CHECK_CLOSE(0.0, dmix.integral(6.0, true, 6.0, false), eps);
CHECK_CLOSE(1.0, dmix.integral(5.0, true, 7.0, true), eps);
CHECK_CLOSE(0.5, dmix.integral(5.0, false, 7.0, false), eps);
CHECK_CLOSE(0.8333333333333334, dmix.integral(5.0, false, 7.0, true), eps);
CHECK_CLOSE(0.6666666666666667, dmix.integral(5.0, true, 7.0, false), eps);
const double expectedMoments[4] = {1.0, 37/6.0, 77/2.0, 1459/6.0};
long double moments[4];
dmix.moments(0, 3, moments);
for (unsigned i=0; i<=3; ++i)
CHECK_CLOSE(expectedMoments[i], moments[i], eps);
const double expectedCMoments[4] = {1.0, 37/6.0, 17/36.0, -2/27.0};
dmix.centralMoments(3, moments);
for (unsigned i=0; i<=3; ++i)
CHECK_CLOSE(expectedCMoments[i], moments[i], eps);
}
TEST(DeltaMixture1D_2)
{
const double eps = 1.0e-12;
const unsigned nPoints = 3;
double coords[nPoints] = {0., 1., 2.};
double weights[nPoints] = {1., 3., 2.};
DeltaMixture1D dmix(5.0, 2.0, coords, weights, nPoints);
io_test(dmix);
CHECK_CLOSE(0.6666666666666666, dmix.cdf(7.0), eps);
CHECK_CLOSE(0.3333333333333333, dmix.exceedance(7.0), eps);
MersenneTwister rng;
int counts[nPoints] = {0};
for (unsigned i=0; i<6000; ++i)
{
double rnd;
dmix.random(rng, &rnd);
CHECK(rnd == 5.0 || rnd == 7.0 || rnd == 9.0);
if (rnd == 5.0)
++counts[0];
if (rnd == 7.0)
++counts[1];
if (rnd == 9.0)
++counts[2];
}
CHECK(std::abs(counts[0] - 1000) < 200);
CHECK(std::abs(counts[2] - 2000) < 300);
CHECK_CLOSE(0.5, dmix.integral(7.0, true, 7.0, true), eps);
CHECK_CLOSE(0.0, dmix.integral(7.0, false, 7.0, false), eps);
CHECK_CLOSE(0.0, dmix.integral(7.0, false, 7.0, true), eps);
CHECK_CLOSE(0.0, dmix.integral(7.0, true, 7.0, false), eps);
CHECK_CLOSE(1.0, dmix.integral(5.0, true, 9.0, true), eps);
CHECK_CLOSE(0.5, dmix.integral(5.0, false, 9.0, false), eps);
CHECK_CLOSE(0.8333333333333334, dmix.integral(5.0, false, 9.0, true), eps);
CHECK_CLOSE(0.6666666666666667, dmix.integral(5.0, true, 9.0, false), eps);
const double expectedMoments[4] = {1.0, 22/3.0, 167/3.0, 1306/3.0};
long double moments[4];
dmix.moments(0, 3, moments);
for (unsigned i=0; i<=3; ++i)
CHECK_CLOSE(expectedMoments[i], moments[i], eps);
const double expectedCMoments[4] = {1.0, 22/3.0, 17/9.0, -16/27.0};
dmix.centralMoments(3, moments);
for (unsigned i=0; i<=3; ++i)
CHECK_CLOSE(expectedCMoments[i], moments[i], eps);
}
TEST(ComparisonDistribution1D)
{
Gauss1D d1(0, 1);
Quadratic1D d2(0.0, 1.0, 0.5, 0.2);
CompositeDistribution1D d3(d2, d1);
ComparisonDistribution1D d4(d3, d1);
for (unsigned i=0; i<100; ++i)
{
const double x = test_rng();
CHECK_CLOSE(d2.density(x), d4.density(x), 1.0e-10);
}
standard_test_01(d2, 1.e-12);
standard_test_01(d4, 1.e-10);
}
TEST(ChiSquare)
{
const double eps = 1.0e-10;
const int n = 7;
const double norm = pow(2.0,-n/2.0)/Gamma(n/2.0);
Gamma1D chisq(0.0, 2.0, n/2.0);
for (unsigned i=0; i<100; ++i)
{
const double x = test_rng()*20;
const double ref = norm*pow(x,n/2.0-1)*exp(-x/2);
CHECK_CLOSE(ref, chisq.density(x), eps);
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:47 AM (22 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983139
Default Alt Text
(75 KB)

Event Timeline