Index: trunk/tests/test_Distributions1D.cc
===================================================================
--- trunk/tests/test_Distributions1D.cc	(revision 596)
+++ trunk/tests/test_Distributions1D.cc	(revision 597)
@@ -1,1860 +1,1870 @@
 #include <complex>
 #include <iostream>
 
 #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/DeltaMixture1D.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 "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 get_cumulants_severini(const EdgeworthSeries1D& s, double cumulants[8])
+    void cumulants_from_moments(const double m[10], double cumulants[10])
+    {
+        for (unsigned i=0; 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_cumulants_severini(const EdgeworthSeries1D& s, double cumulants[10])
     {
         assert(s.method() == EDGEWORTH_SEVERINI);
-        double m[8];
+        double m[10];
         GaussHermiteQuadrature quad(128);
-        for (unsigned i=0; i<8; ++i)
+        for (unsigned i=0; i<10; ++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];
+        cumulants_from_moments(m, cumulants);
     }
 
-    void get_cumulants_classical(const EdgeworthSeries1D& s, double cumulants[8])
+    void get_cumulants_classical(const EdgeworthSeries1D& s, double cumulants[10])
     {
         assert(s.method() == EDGEWORTH_CLASSICAL);
-        double m[8];
+        double m[10];
         GaussHermiteQuadrature quad(128);
-        for (unsigned i=0; i<8; ++i)
+        for (unsigned i=0; i<10; ++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];
+        cumulants_from_moments(m, 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 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<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[8];
+        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-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];
+        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_cumulants_severini(es1, empCum);
         CHECK_CLOSE(1.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_cumulants_severini(es2, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es3, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es4, empCum);
         CHECK_CLOSE(1.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[8];
+        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_cumulants_severini(es1, empCum);
         CHECK_CLOSE(1.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_cumulants_severini(es2, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es3, empCum);
         CHECK_CLOSE(1.0, empCum[0], 1.e-14);
         CHECK_CLOSE(k1, empCum[1], 1.e-14);
         CHECK_CLOSE(k2, 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);
         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_cumulants_classical(es4, empCum);
         CHECK_CLOSE(1.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-14);
         CHECK_CLOSE(-35*k3*k4, empCum[7], 1.e-14);
     }
 
     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[8];
+        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_cumulants_severini(es1, empCum);
         CHECK_CLOSE(1.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_cumulants_severini(es2, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es3, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es4, empCum);
         CHECK_CLOSE(1.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-14);
         CHECK_CLOSE(0.0, empCum[7], 1.e-14);
     }
 
     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[8];
+        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_cumulants_severini(es1, empCum);
         CHECK_CLOSE(1.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_cumulants_severini(es2, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es3, empCum);
         CHECK_CLOSE(1.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_cumulants_classical(es4, empCum);
         CHECK_CLOSE(1.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-14);
         CHECK_CLOSE(0.0, empCum[7], 1.e-14);
+        CHECK_CLOSE(0.0, empCum[8], 1.e-12);
+        CHECK_CLOSE(-126*k4*k5 - 84*k3*k6, empCum[9], 1.e-10);
     }
 
     TEST(Tabulated1D)
     {
         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(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(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);
     }
 }