Index: trunk/tests/test_Distributions1D.cc
===================================================================
--- trunk/tests/test_Distributions1D.cc	(revision 743)
+++ trunk/tests/test_Distributions1D.cc	(revision 744)
@@ -1,2084 +1,2136 @@
 #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/PolynomialDistro1D.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);
     }
 
+    class DummyFunct : public Functor1<double,double>
+    {
+    public:
+        inline DummyFunct(double a, double b)
+            : a_(a), b_(b) {}
+
+        inline double operator()(const double& x) const
+        {
+            const double twoxm1 = 2*x - 1;
+            const double p1 = twoxm1;
+            const double p2 = 0.5*(3.0*twoxm1*twoxm1 - 1.0);
+            return exp(a_*p1 + b_*p2);
+        }
+
+    private:
+        double a_;
+        double b_;
+    };
+    
     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_beta_symbeta_compare)
     {
         const double eps = 1.0e-14;
 
         SymmetricBeta1D sb(0.0, 2.0, 3);
         Beta1D b(-2.0, 4.0, 4.0, 4.0);
 
         for (unsigned i=0; i<100; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             CHECK_CLOSE(b.density(x), sb.density(x), eps);
         }
     }
 
     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(PolynomialDistro1D_1)
     {
         const double coeffs[] = {0.2, 0.1};
         PolynomialDistro1D pd(-1.0, 2.0, coeffs, sizeof(coeffs)/sizeof(coeffs[0]));
         standard_test(pd, 1.0, 1.e-10);
         invcdf_test(pd, 1.e-12);
     }
 
     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);
+        for (unsigned i=0; i<10; ++i)
+        {
+            const double a = (test_rng()-0.5)/10.0;
+            const double b = (test_rng()-0.5)/10.0;
+            Quadratic1D q1(0.0, 1.0, a, b);
+            for (unsigned j=0; j<10; ++j)
+            {
+                const double x = test_rng();
+                const double twoxm1 = 2*x - 1;
+                const double p1 = twoxm1;
+                const double p2 = 0.5*(3.0*twoxm1*twoxm1 - 1.0);
+                const double ref = 1.0 + a*p1 + b*p2;
+                CHECK_CLOSE(ref, q1.density(x), 1.0e-12);
+            }
+        }
     }
 
     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(LogQuadratic1D_7)
+    {
+        for (unsigned i=0; i<10; ++i)
+        {
+            const double a = (test_rng()-0.5)/10.0;
+            const double b = (test_rng()-0.5)/10.0;
+            LogQuadratic1D q1(0.0, 1.0, a, b);
+
+            GaussLegendreQuadrature glq(1024);
+            const double norm = glq.integrate(DummyFunct(a, b), 0.0, 1.0);
+
+            for (unsigned j=0; j<10; ++j)
+            {
+                const double x = test_rng();
+                const double twoxm1 = 2*x - 1;
+                const double p1 = twoxm1;
+                const double p2 = 0.5*(3.0*twoxm1*twoxm1 - 1.0);
+                const double ref = exp(a*p1 + b*p2)/norm;
+                CHECK_CLOSE(ref, q1.density(x), 1.0e-12);
+            }
+        }
+    }
+
     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);
         }
     }
 }
Index: trunk/npstat/stat/Distributions1D.cc
===================================================================
--- trunk/npstat/stat/Distributions1D.cc	(revision 743)
+++ trunk/npstat/stat/Distributions1D.cc	(revision 744)
@@ -1,2820 +1,2820 @@
 #include <cmath>
 #include <cfloat>
 #include <algorithm>
 #include <stdexcept>
 
 #include "geners/binaryIO.hh"
 
 #include "npstat/nm/MathUtils.hh"
 #include "npstat/nm/SpecialFunctions.hh"
 #include "npstat/nm/interpolate.hh"
 
 #include "npstat/stat/Distributions1D.hh"
 #include "npstat/stat/StatUtils.hh"
 #include "npstat/stat/distributionReadError.hh"
 
 #define SQR2PI 2.5066282746310005
 #define SQRT2  1.41421356237309505
 #define SQRPI  1.77245385090551603
 
 #define SQRT2L 1.414213562373095048801689L
 #define SQRPIL 1.77245385090551602729816748L
 #define TWOPIL 6.28318530717958647692528676656L
 
 static long double inverseErf(const long double fval)
 {
     long double x = npstat::inverseGaussCdf((fval + 1.0L)/2.0L)/SQRT2L;
     for (unsigned i=0; i<2; ++i)
     {
         const long double guessed = erfl(x);
         const long double deri = 2.0L/SQRPIL*expl(-x*x);
         x += (fval - guessed)/deri;
     }
     return x;
 }
 
 static unsigned improved_random(npstat::AbsRandomGenerator& g,
                                 long double* generatedRandom)
 {
     const long double extra = sqrt(DBL_EPSILON);
 
     long double u = 0.0L;
     unsigned callcount = 0;
     while (u <= 0.0L || u >= 1.0L)
     {
         u = g()*(1.0L + extra) - extra/2.0L;
         u += (g() - 0.5L)*extra;
         callcount += 2U;
     }
     *generatedRandom = u;
     return callcount;
 }
 
 static unsigned gauss_random(const double mean, const double sigma,
                              npstat::AbsRandomGenerator& g,
                              double* generatedRandom)
 {
     assert(generatedRandom);
     long double r1 = 0.0L, r2 = 0.0L;
     const unsigned calls = improved_random(g, &r1) + improved_random(g, &r2);
     *generatedRandom = mean + sigma*sqrtl(-2.0L*logl(r1))*sinl(TWOPIL*(r2-0.5L));
     return calls;
 }
 
 // static unsigned gauss_random(const double mean, const double sigma,
 //                              npstat::AbsRandomGenerator& g,
 //                              double* generatedRandom)
 // {
 //     assert(generatedRandom);
 //     long double r1 = 0.0L;
 //     const unsigned count = improved_random(g, &r1);
 //     *generatedRandom = mean + sigma*SQRT2*inverseErf(2.0L*r1 - 1.0L);
 //     return count;
 // }
 
 namespace npstat {
     bool SymmetricBeta1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, n_);
         return !os.fail();
     }
 
     SymmetricBeta1D* SymmetricBeta1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<SymmetricBeta1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             double n;
             gs::read_pod(in, &n);
             if (!in.fail())
                 return new SymmetricBeta1D(location, scale, n);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool SymmetricBeta1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const SymmetricBeta1D& r = 
             static_cast<const SymmetricBeta1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && n_ == r.n_;
     }
 
     SymmetricBeta1D::SymmetricBeta1D(const double location,
                                      const double scale,
                                      const double power)
         : AbsScalableDistribution1D(location, scale),
           n_(power), intpow_(-1)
     {
         norm_ = calculateNorm();
     }
 
     SymmetricBeta1D::SymmetricBeta1D(const double location,
                                      const double scale,
                                      const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           n_(params[0]), intpow_(-1)
     {
         norm_ = calculateNorm();
     }
 
     double SymmetricBeta1D::calculateNorm()
     {
         static const double normcoeffs[11] = {
             0.5, 0.75, 0.9375, 1.09375, 1.23046875, 1.353515625,
             1.46630859375, 1.571044921875, 1.6692352294921875,
             1.76197052001953125, 1.85006904602050781};
 
         if (n_ <= -1.0) throw std::invalid_argument(
             "In npstat::SymmetricBeta1D::calculateNorm: "
             "invalid power parameter");
 
         const int intpow = static_cast<int>(floor(n_));
         if (static_cast<double>(intpow) == n_ &&
             intpow >= 0 && intpow <= 10)
         {
             intpow_ = intpow;
             return normcoeffs[intpow];
         }
         else
             return Gamma(1.5 + n_)/sqrt(M_PI)/Gamma(1.0 + n_);
     }
 
     double SymmetricBeta1D::unscaledDensity(const double x) const
     {
         const double oneminusrsq = 1.0 - x*x;
         if (oneminusrsq <= 0.0)
             return 0.0;
         else
         {
             double fcn;
             switch (intpow_)
             {
             case 0:
                 fcn = 1.0;
                 break;
             case 1:
                 fcn = oneminusrsq;
                 break;
             case 2:
                 fcn = oneminusrsq*oneminusrsq;
                 break;
             case 3:
                 fcn = oneminusrsq*oneminusrsq*oneminusrsq;
                 break;
             case 4:
                 {
                     const double tmp2 = oneminusrsq*oneminusrsq;
                     fcn = tmp2*tmp2;
                 }
                 break;
             case 5:
                 {
                     const double tmp2 = oneminusrsq*oneminusrsq;
                     fcn = tmp2*tmp2*oneminusrsq;
                 }
                 break;
             case 6:
                 {
                     const double tmp2 = oneminusrsq*oneminusrsq;
                     fcn = tmp2*tmp2*tmp2;
                 }
                 break;
             default:
                 fcn = pow(oneminusrsq, n_);
                 break;
             }
             return norm_*fcn;
         }
     }
 
     double SymmetricBeta1D::unscaledCdf(const double x) const
     {
         if (x >= 1.0)
             return 1.0;
         else if (x <= -1.0)
             return 0.0;
         else if (n_ == 0.0)
             return (x + 1.0)/2.0;
         else
             return incompleteBeta(n_+1.0, n_+1.0, (x + 1.0)/2.0);
     }
 
     double SymmetricBeta1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::SymmetricBeta1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return -1.0;
         else if (r1 == 1.0)
             return 1.0;
         else
         {
             double r;
             if (n_ == 0.0)
                 r = r1*2.0 - 1.0;
             else
                 r = 2.0*inverseIncompleteBeta(n_+1.0, n_+1.0, r1) - 1.0;
             if (r < -1.0)
                 r = -1.0;
             else if (r > 1.0)
                 r = 1.0;
             return r;
         }
     }
 
     bool Beta1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, alpha_);
         gs::write_pod(os, beta_);
         return !os.fail();
     }
 
     Beta1D* Beta1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Beta1D>());
         current.ensureSameId(id);
 
         double location, scale, a, b;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &a);
             gs::read_pod(in, &b);
             if (!in.fail())
                 return new Beta1D(location, scale, a, b);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool Beta1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Beta1D& r = static_cast<const Beta1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                alpha_ == r.alpha_ && beta_ == r.beta_;
     }
 
     Beta1D::Beta1D(const double location, const double scale,
                    const double pa, const double pb)
         : AbsScalableDistribution1D(location, scale),
           alpha_(pa),
           beta_(pb)
     {
         norm_ = calculateNorm();
     }
 
     Beta1D::Beta1D(const double location,
                    const double scale,
                    const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           alpha_(params[0]),
           beta_(params[1])
     {
         norm_ = calculateNorm();
     }
 
     double Beta1D::calculateNorm() const
     {
         if (!(alpha_ > 0.0 && beta_ > 0.0)) throw std::invalid_argument(
             "In npstat::Beta1D::calculateNorm: invalid power parameters");
         return Gamma(alpha_ + beta_)/Gamma(alpha_)/Gamma(beta_);
     }
 
     double Beta1D::unscaledDensity(const double x) const
     {
         if (x <= 0.0 || x >= 1.0)
             return 0.0;
         else if (alpha_ == 1.0 && beta_ == 1.0)
             return 1.0;
         else
             return norm_*pow(x, alpha_-1.0)*pow(1.0-x, beta_-1.0);
     }
 
     double Beta1D::unscaledCdf(const double x) const
     {
         if (x >= 1.0)
             return 1.0;
         else if (x <= 0.0)
             return 0.0;
         else if (alpha_ == 1.0 && beta_ == 1.0)
             return x;
         else
             return incompleteBeta(alpha_, beta_, x);
     }
 
     double Beta1D::unscaledExceedance(const double x) const
     {
         return 1.0 - unscaledCdf(x);
     }
 
     double Beta1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Beta1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         else if (r1 == 1.0)
             return 1.0;
         else if (alpha_ == 1.0 && beta_ == 1.0)
             return r1;
         else
             return inverseIncompleteBeta(alpha_, beta_, r1);
     }
 
     Gamma1D::Gamma1D(const double location, const double scale,
                      const double a)
         : AbsScalableDistribution1D(location, scale),
           alpha_(a)
     {
         initialize();
     }
 
     Gamma1D::Gamma1D(double location, double scale,
                      const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           alpha_(params[0])
     {
         initialize();
     }
 
     void Gamma1D::initialize()
     {
         if (!(alpha_ > 0.0)) throw std::invalid_argument(
             "In npstat::Gamma1D::initialize: invalid power parameter");
         norm_ = 1.0/Gamma(alpha_);
     }
 
     bool Gamma1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Gamma1D& r = static_cast<const Gamma1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && alpha_ == r.alpha_;
     }
 
     double Gamma1D::unscaledDensity(const double x) const
     {
         if (x > 0.0)
             return norm_*pow(x, alpha_-1.0)*exp(-x);
         else
             return 0.0;
     }
 
     double Gamma1D::unscaledCdf(const double x) const
     {
         if (x > 0.0)
             return incompleteGamma(alpha_, x);
         else
             return 0.0;
     }
 
     double Gamma1D::unscaledExceedance(const double x) const
     {
         if (x > 0.0)
             return incompleteGammaC(alpha_, x);
         else
             return 1.0;
     }
 
     double Gamma1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Gamma1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         return inverseIncompleteGamma(alpha_, r1);
     }
 
     bool Gamma1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, alpha_);
         return !os.fail();
     }
 
     Gamma1D* Gamma1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Gamma1D>());
         current.ensureSameId(id);
 
         double location, scale, a;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &a);
             if (!in.fail())
                 return new Gamma1D(location, scale, a);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     Gauss1D::Gauss1D(const double location, const double scale)
         : AbsScalableDistribution1D(location, scale),
           xmin_(inverseGaussCdf(0.0)), xmax_(inverseGaussCdf(1.0))
     {
     }
 
     Gauss1D::Gauss1D(const double location, const double scale,
                      const std::vector<double>& /* params */)
         : AbsScalableDistribution1D(location, scale),
           xmin_(inverseGaussCdf(0.0)), xmax_(inverseGaussCdf(1.0))
     {
     }
 
     double Gauss1D::unscaledDensity(const double x) const
     {
         if (x < xmin_ || x > xmax_)
             return 0.0;
         else
             return exp(-x*x/2.0)/SQR2PI;
     }
 
     unsigned Gauss1D::random(AbsRandomGenerator& g,
                              double* generatedRandom) const
     {
         return gauss_random(location(), scale(), g, generatedRandom);
     }
 
     bool Gauss1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Gauss1D* Gauss1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Gauss1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Gauss1D(location, scale);
     }
 
     double Gauss1D::unscaledCdf(const double x) const
     {
         if (x <= xmin_)
             return 0.0;
         if (x >= xmax_)
             return 1.0;
         if (x < 0.0)
             return erfc(-x/SQRT2)/2.0;
         else
             return (1.0 + erf(x/SQRT2))/2.0;
     }
 
     double Gauss1D::unscaledExceedance(const double x) const
     {
         if (x <= xmin_)
             return 1.0;
         if (x >= xmax_)
             return 0.0;
         if (x > 0.0)
             return erfc(x/SQRT2)/2.0;
         else
             return (1.0 - erf(x/SQRT2))/2.0;
     }
 
     double Gauss1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Gauss1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         return inverseGaussCdf(r1);
     }
 
     double Uniform1D::unscaledDensity(const double x) const
     {
         if (x >= 0.0 && x <= 1.0)
             return 1.0;
         else
             return 0.0;
     }
 
     bool Uniform1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Uniform1D* Uniform1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Uniform1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Uniform1D(location, scale);
     }
 
     double Uniform1D::unscaledCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0;
         else if (x >= 1.0)
             return 1.0;
         else
             return x;
     }
 
     double Uniform1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Uniform1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         return r1;
     }
 
     double IsoscelesTriangle1D::unscaledDensity(const double x) const
     {
         if (x > -1.0 && x < 1.0)
             return 1.0 - fabs(x);
         else
             return 0.0;
     }
 
     bool IsoscelesTriangle1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
     
     IsoscelesTriangle1D* IsoscelesTriangle1D::read(
         const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<IsoscelesTriangle1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new IsoscelesTriangle1D(location, scale);
     }
 
     double IsoscelesTriangle1D::unscaledCdf(const double x) const
     {
         if (x <= -1.0)
             return 0.0;
         else if (x >= 1.0)
             return 1.0;
         else if (x <= 0.0)
         {
             const double tmp = 1.0 + x;
             return 0.5*tmp*tmp;
         }
         else
         {
             const double tmp = 1.0 - x;
             return 1.0 - 0.5*tmp*tmp;
         }
     }
 
     double IsoscelesTriangle1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::IsoscelesTriangle1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return -1.0;
         else if (r1 == 1.0)
             return 1.0;
         else if (r1 <= 0.5)
             return sqrt(2.0*r1) - 1.0;
         else
             return 1.0 - sqrt((1.0 - r1)*2.0);
     }
 
     bool Exponential1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Exponential1D* Exponential1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Exponential1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Exponential1D(location, scale);
     }
 
     double Exponential1D::unscaledDensity(const double x) const
     {
         if (x < 0.0)
             return 0.0;
         const double eval = exp(-x);
         return eval < DBL_MIN ? 0.0 : eval;
     }
 
     double Exponential1D::unscaledCdf(const double x) const
     {
         return x > 0.0 ? 1.0 - exp(-x) : 0.0;
     }
 
     double Exponential1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Exponential1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 1.0)
             return -log(DBL_MIN);
         else
             return -log(1.0 - r1);
     }
 
     double Exponential1D::unscaledExceedance(const double x) const
     {
         if (x < 0.0)
             return 1.0;
         const double eval = exp(-x);
         return eval < DBL_MIN ? 0.0 : eval;
     }
 
     bool Logistic1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Logistic1D* Logistic1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Logistic1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Logistic1D(location, scale);
     }
 
     double Logistic1D::unscaledDensity(const double x) const
     {
         const double eval = exp(-x);
         if (eval < DBL_MIN)
             return 0.0;
         else
         {
             const double tmp = 1.0 + eval;
             return eval/tmp/tmp;
         }
     }
 
     double Logistic1D::unscaledCdf(const double x) const
     {
         const double lmax = -log(DBL_MIN);
         if (x <= -lmax)
             return 0.0;
         else if (x >= lmax)
             return 1.0;
         else
             return 1.0/(1.0 + exp(-x));
     }
 
     double Logistic1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Logistic1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         const double lmax = -log(DBL_MIN);
         if (r1 == 0.0)
             return -lmax;
         else if (r1 == 1.0)
             return lmax;
         else
             return log(r1/(1.0 - r1));
     }
 
     double Logistic1D::unscaledExceedance(const double x) const
     {
         const double lmax = -log(DBL_MIN);
         if (x >= lmax)
             return 0.0;
         else if (x <= -lmax)
             return 1.0;
         else
         {
             const double eval = exp(-x);
             return eval/(1.0 + eval);
         }
     }
 
     bool Quadratic1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, a_);
         gs::write_pod(os, b_);
         return !os.fail();
     }
 
     Quadratic1D* Quadratic1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Quadratic1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             double a, b;
             gs::read_pod(in, &a);
             gs::read_pod(in, &b);
             if (!in.fail())
                 return new Quadratic1D(location, scale, a, b);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool Quadratic1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Quadratic1D& r = static_cast<const Quadratic1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) &&
                a_ == r.a_ && b_ == r.b_;
     }
 
     Quadratic1D::Quadratic1D(const double location, const double scale,
                              const double a, const double b)
         : AbsScalableDistribution1D(location, scale), a_(a), b_(b)
     {
         verifyNonNegative();
     }
 
     Quadratic1D::Quadratic1D(const double location, const double scale,
                              const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           a_(params[0]),
           b_(params[1])
     {
         verifyNonNegative();
     }
 
     void Quadratic1D::verifyNonNegative()
     {
-        const double a = 2.0*a_;
-        const double b = 2.0*b_;
+        const double a = a_;
+        const double b = b_;
         if (b == 0.0)
         {
             if (fabs(a) > 1.0)
                 throw std::invalid_argument(
                     "In npstat::Quadratic1D::verifyNonNegative:"
                     " invalid distribution parameters");
         }
         else
         {
             double x1 = 0.0, x2 = 0.0;
             const double sixb = 6*b;
             if (solveQuadratic((2*a-sixb)/sixb, (1-a+b)/sixb, &x1, &x2))
             {
                 if (!(fabs(x1 - 0.5) >= 0.5 && fabs(x2 - 0.5) >= 0.5))
                     throw std::invalid_argument(
                         "In npstat::Quadratic1D::verifyNonNegative:"
                         " invalid distribution parameters");
             }
         }
     }
 
     double Quadratic1D::unscaledDensity(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
         else
-            return 1.0 + 2.0*(b_ - a_ + x*(2.0*a_ + 6.0*b_*(x - 1.0)));
+            return 1.0 + b_ - a_ + x*(2.0*a_ + 6.0*b_*(x - 1.0));
     }
 
     double Quadratic1D::unscaledCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0;
         else if (x >= 1.0)
             return 1.0;
         else
-            return x*(1.0 + 2.0*(b_ - a_ + x*(a_ - 3.0*b_ + 2.0*b_*x)));
+            return x*(1.0 + b_ - a_ + x*(a_ - 3.0*b_ + 2.0*b_*x));
     }
 
     double Quadratic1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Quadratic1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         else if (r1 == 1.0)
             return 1.0;
         else
         {
-            const double a = 2.0*a_;
-            const double b = 2.0*b_;
+            const double a = a_;
+            const double b = b_;
             if (b == 0.0)
             {
                 if (a == 0.0)
                     return r1;
                 else
                 {
                     double x0 = 0.0, x1 = 0.0;
                     const unsigned n = solveQuadratic(
                         (1.0 - a)/a, -r1/a, &x0, &x1);
                     if (!n) throw std::runtime_error(
                         "In npstat::Quadratic1D::unscaledQuantile: "
                         "no solutions");
                     if (fabs(x0 - 0.5) < fabs(x1 - 0.5))
                         return x0;
                     else
                         return x1;
                 }
             }
             else
             {
                 const double twob = 2*b;
                 double x[3] = {0.0};
                 const unsigned n = solveCubic(
                     (a - 3*b)/twob, (1 - a + b)/twob, -r1/twob, x);
                 if (n == 1U)
                     return x[0];
                 else
                 {
                     unsigned ibest = 0;
                     double dbest = fabs(x[0] - 0.5);
                     for (unsigned i=1; i<n; ++i)
                         if (fabs(x[i] - 0.5) < dbest)
                         {
                             ibest = i;
                             dbest = fabs(x[i] - 0.5);
                         }
                     return x[ibest];
                 }
             }
         }
     }
 
     bool LogQuadratic1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, a_);
         gs::write_pod(os, b_);
         return !os.fail();
     }
 
     LogQuadratic1D* LogQuadratic1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<LogQuadratic1D>());
         current.ensureSameId(id);
 
         double location, scale, a, b;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &a);
             gs::read_pod(in, &b);
             if (!in.fail())
                 return new LogQuadratic1D(location, scale, a, b);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool LogQuadratic1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const LogQuadratic1D& r = static_cast<const LogQuadratic1D&>(
             otherBase);
         return AbsScalableDistribution1D::isEqual(r) &&
                a_ == r.a_ && b_ == r.b_;
     }
 
     LogQuadratic1D::LogQuadratic1D(const double location, const double scale,
                                    const double a, const double b)
         : AbsScalableDistribution1D(location, scale), a_(a), b_(b)
     {
         normalize();
     }
 
     LogQuadratic1D::LogQuadratic1D(const double location, const double scale,
                                    const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           a_(params[0]),
           b_(params[1])
     {
         normalize();
     }
 
     inline long double LogQuadratic1D::quadInteg(const long double x) const
     {
         return dawsonIntegral(x)*expl(x*x);
     }
 
     void LogQuadratic1D::normalize()
     {
         ref_ = 0.0L;
         range_ = 1.0L;
         k_ = 0.0;
         s_ = 0.0;
         norm_ = 1.0;
 
-        const double b = 2.0*b_;
-        const double a = 2.0*a_;
+        const double b = b_;
+        const double a = a_;
 
         if (b > DBL_EPSILON)
         {
             k_ = sqrt(6.0*b);
             s_ = 0.5 - a/(6.0*b);
             ref_ = quadInteg(k_*s_);
             range_ = ref_ - quadInteg(k_*(s_ - 1.0));
             norm_ = k_/range_;
         }
         else if (b < -DBL_EPSILON)
         {
             k_ = sqrt(-6.0*b);
             s_ = 0.5 - a/(6.0*b);
             ref_ = erfl(k_*s_);
             range_ = ref_ - erfl(k_*(s_ - 1.0));
             norm_ = 2.0*k_/SQRPI/range_;
         }
         else if (fabs(a) > DBL_EPSILON)
         {
             range_ = expl(2.0L*a) - 1.0L;
             if (fabs(a) > 1.e-10)
                 norm_ = a/sinh(a);
         }
     }
 
     double LogQuadratic1D::unscaledDensity(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
 
-        const double b = 2.0*b_;
+        const double b = b_;
         if (fabs(b) > DBL_EPSILON)
         {
             const double delta = x - s_;
             return norm_*exp(6.0*b*delta*delta);
         }
 
-        const double a = 2.0*a_;
+        const double a = a_;
         if (fabs(a) > DBL_EPSILON)
             return norm_*exp((2.0*x - 1.0)*a);
         else
             return 1.0;
     }
 
     double LogQuadratic1D::unscaledCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0;
         else if (x >= 1.0)
             return 1.0;
         else
         {
-            const double b = 2.0*b_;
-            const double a = 2.0*a_;
+            const double b = b_;
+            const double a = a_;
             if (b > DBL_EPSILON)
                 return (ref_ - quadInteg(k_*(s_ - x)))/range_;
             else if (b < -DBL_EPSILON)
                 return (ref_ - erfl(k_*(s_ - x)))/range_;
             else if (fabs(a) > DBL_EPSILON)
                 return (expl(2.0L*a*x) - 1.0L)/range_;
             else
                 return x;
         }
     }
 
     double LogQuadratic1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::LogQuadratic1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         else if (r1 == 1.0)
             return 1.0;
         else
         {
-            const double b = 2.0*b_;
-            const double a = 2.0*a_;
+            const double b = b_;
+            const double a = a_;
             double q = 0.0;
             if (b > DBL_EPSILON)
                 q = s_ - inverseExpsqIntegral(ref_ - r1*range_)/k_;
             else if (b < -DBL_EPSILON)
                 q = s_ - inverseErf(ref_ - r1*range_)/k_;
             else if (fabs(a) > DBL_EPSILON)
                 q = logl(r1*range_ + 1.0L)/2.0/a;
             else
                 q = r1;
 
             if (q < 0.0)
                 q = 0.0;
             else if (q > 1.0)
                 q = 1.0;
             return q;
         }
     }
 
     bool TruncatedGauss1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, nsigma_);
         return !os.fail();
     }
 
     TruncatedGauss1D* TruncatedGauss1D::read(
         const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<TruncatedGauss1D>());
         current.ensureSameId(id);
 
         double location, scale, nsig;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &nsig);
             if (!in.fail())
                 return new TruncatedGauss1D(location, scale, nsig);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool TruncatedGauss1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const TruncatedGauss1D& r = 
             static_cast<const TruncatedGauss1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && nsigma_ == r.nsigma_;
     }
 
     void TruncatedGauss1D::initialize()
     {
         if (nsigma_ <= 0.0) throw std::invalid_argument(
             "In npstat::TruncatedGauss1D::initialize: "
             "invalid truncation parameter");
         const double maxSig = inverseGaussCdf(1.0);
         if (nsigma_ >= maxSig)
         {
             nsigma_ = maxSig;
             cdf0_ = 0.0;
             norm_ = 1.0;
         }
         else
         {
             cdf0_ = erfc(nsigma_/SQRT2)/2.0;
             const double u = (1.0 + erf(nsigma_/SQRT2))/2.0;
             norm_ = 1.0/(u - cdf0_);
         }
     }
 
     TruncatedGauss1D::TruncatedGauss1D(const double location,
                                        const double scale,
                                        const double i_nsigma)
         : AbsScalableDistribution1D(location, scale),
           nsigma_(i_nsigma)
     {
         initialize();
     }
 
     TruncatedGauss1D::TruncatedGauss1D(const double location,
                                        const double scale,
                                        const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           nsigma_(params[0])
     {
         initialize();
     }
 
     double TruncatedGauss1D::unscaledDensity(const double x) const
     {
         if (fabs(x) > nsigma_)
             return 0.0;
         else
             return norm_*exp(-x*x/2.0)/SQR2PI;
     }
 
     unsigned TruncatedGauss1D::random(AbsRandomGenerator& g,
                                       double* generatedRandom) const
     {
         const double m = location();
         const double s = scale();
         unsigned cnt = gauss_random(m, s, g, generatedRandom);
         while (fabs(*generatedRandom - m) > nsigma_*s)
             cnt += gauss_random(m, s, g, generatedRandom);
         return cnt;
     }
 
     double TruncatedGauss1D::unscaledCdf(const double x) const
     {
         if (x <= -nsigma_)
             return 0.0;
         else if (x >= nsigma_)
             return 1.0;
         else if (x < 0.0)
             return (erfc(-x/SQRT2)/2.0 - cdf0_)*norm_;
         else
             return ((1.0 + erf(x/SQRT2))/2.0 - cdf0_)*norm_;
     }
 
     double TruncatedGauss1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::TruncatedGauss1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return -nsigma_;
         else if (r1 == 1.0)
             return nsigma_;
         else
             return inverseGaussCdf(r1/norm_ + cdf0_);
     }
 
     bool MirroredGauss1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const MirroredGauss1D& r = 
             static_cast<const MirroredGauss1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
             mu0_ == r.mu0_ && sigma0_ == r.sigma0_;
     }
 
     MirroredGauss1D::MirroredGauss1D(const double location, const double scale,
                                      const double mean, double const sigma)
         : AbsScalableDistribution1D(location, scale),
           mu0_(mean),
           sigma0_(sigma)
     {
         validate();
     }
 
     MirroredGauss1D::MirroredGauss1D(const double location,
                                      const double scale,
                                      const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           mu0_(params[0]),
           sigma0_(params[1])
     {
         validate();
     }
 
     bool MirroredGauss1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, mu0_);
         gs::write_pod(os, sigma0_);
         return !os.fail();
     }
 
     MirroredGauss1D* MirroredGauss1D::read(
         const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<MirroredGauss1D>());
         current.ensureSameId(id);
 
         double location, scale, mu, sig;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &mu);
             gs::read_pod(in, &sig);
             if (!in.fail())
                 return new MirroredGauss1D(location, scale, mu, sig);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     double MirroredGauss1D::unscaledDensity(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
         Gauss1D g(x, sigma0_);
         long double acc = g.density(mu0_)*1.0L + g.density(-mu0_);
         for (unsigned k=1; ; ++k)
         {
             const long double old = acc;
             acc += g.density(2.0*k + mu0_);
             acc += g.density(2.0*k - mu0_);
             acc += g.density(-2.0*k + mu0_);
             acc += g.density(-2.0*k - mu0_);
             if (old == acc)
                 break;
         }
         return acc;
     }
 
     long double MirroredGauss1D::ldCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0L;
         else if (x >= 1.0)
             return 1.0L;
         else
         {
             long double acc = 0.0L;
             {
                 Gauss1D g(mu0_, sigma0_);
                 acc += (g.cdf(x) - g.cdf(0.0));
             }
             {
                 Gauss1D g(-mu0_, sigma0_);
                 acc += (g.cdf(x) - g.cdf(0.0));
             }
             for (unsigned k=1; ; ++k)
             {
                 const long double old = acc;
                 {
                     Gauss1D g(2.0*k + mu0_, sigma0_);
                     acc += (g.cdf(x) - g.cdf(0.0));
                 }
                 {
                     Gauss1D g(2.0*k - mu0_, sigma0_);
                     acc += (g.cdf(x) - g.cdf(0.0));
                 }
                 {
                     Gauss1D g(-2.0*k + mu0_, sigma0_);
                     acc -= (g.exceedance(x) - g.exceedance(0.0));
                 }
                 {
                     Gauss1D g(-2.0*k - mu0_, sigma0_);
                     acc -= (g.exceedance(x) - g.exceedance(0.0));
                 }
                 if (old == acc)
                     break;
             }
             return acc;
         }
     }
 
     double MirroredGauss1D::unscaledCdf(const double x) const
     {
         return ldCdf(x);
     }
 
     double MirroredGauss1D::unscaledExceedance(const double x) const
     {
         return 1.0L - ldCdf(x);
     }
 
     double MirroredGauss1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::MirroredGauss1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         else if (r1 == 1.0)
             return 1.0;
         else
         {
             const long double ldr1 = r1;
             double xmin = 0.0;
             double xmax = 1.0;
             for (unsigned i=0; i<1000; ++i)
             {
                 if ((xmax - xmin)/xmax <= 2.0*DBL_EPSILON)
                     break;
                 const double xtry = (xmin + xmax)/2.0;
                 const long double ld = ldCdf(xtry);
                 if (ld == ldr1)
                     return xtry;
                 else if (ld > ldr1)
                     xmax = xtry;
                 else
                     xmin = xtry;
             }
             return (xmin + xmax)/2.0;
         }
     }
 
     void MirroredGauss1D::validate()
     {
         if (mu0_ < 0.0 || mu0_ > 1.0) throw std::invalid_argument(
             "In MirroredGauss1D::validate: interval mean must be within [0, 1]");
         if (sigma0_ <= 0.0) throw std::invalid_argument(
             "In MirroredGauss1D::validate: interval sigma must be positive");
     }
 
     BifurcatedGauss1D::BifurcatedGauss1D(
         const double location, const double scale,
         const double i_leftSigmaFraction,
         const double i_nSigmasLeft, const double i_nSigmasRight)
         : AbsScalableDistribution1D(location, scale),
           leftSigma_(i_leftSigmaFraction*2.0),
           rightSigma_(2.0 - leftSigma_),
           nSigmasLeft_(i_nSigmasLeft),
           nSigmasRight_(i_nSigmasRight)
     {
         initialize();
     }
 
     BifurcatedGauss1D::BifurcatedGauss1D(const double location,
                                          const double scale,
                                          const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           leftSigma_(params[0]*2.0),
           rightSigma_(2.0 - leftSigma_),
           nSigmasLeft_(params[1]),
           nSigmasRight_(params[2])
     {
         initialize();
     }
 
     void BifurcatedGauss1D::initialize()
     {
         if (leftSigma_ < 0.0 || rightSigma_ < 0.0)
             throw std::invalid_argument(
                 "In npstat::BifurcatedGauss1D::initialize: "
                 "invalid left sigma fraction");
         if (nSigmasLeft_ < 0.0) throw std::invalid_argument(
             "In npstat::BifurcatedGauss1D::initialize: "
             "invalid left truncation parameter");
         if (nSigmasRight_ < 0.0) throw std::invalid_argument(
             "In npstat::BifurcatedGauss1D::initialize: "
             "invalid right truncation parameter");
         if (nSigmasLeft_ + nSigmasRight_ == 0.0) throw std::invalid_argument(
             "In npstat::BifurcatedGauss1D::initialize: "
             "both truncation parameters can not be 0");
 
         const double maxNSigma = inverseGaussCdf(1.0);
         if (nSigmasRight_ > maxNSigma)
             nSigmasRight_ = maxNSigma;
         if (nSigmasLeft_ > maxNSigma)
             nSigmasLeft_ = maxNSigma;
 
         cdf0Left_ = erfc(nSigmasLeft_/SQRT2)/2.0;
         cdf0Right_ = (1.0 + erf(nSigmasRight_/SQRT2))/2.0;
         assert(cdf0Right_ > cdf0Left_);
 
         const double leftArea = (0.5 - cdf0Left_)*leftSigma_;
         const double rightArea = (cdf0Right_ - 0.5)*rightSigma_;
         norm_ = 1.0/(leftArea + rightArea);
         leftCdfFrac_ = leftArea/(leftArea + rightArea);
     }
 
     double BifurcatedGauss1D::unscaledDensity(const double x) const
     {
         if (x == 0.0)
             return norm_/SQR2PI;
         else if (x > 0.0)
         {
             if (x > rightSigma_*nSigmasRight_)
                 return 0.0;
             else
             {
                 const double dx = x/rightSigma_;
                 return norm_*exp(-dx*dx/2.0)/SQR2PI;
             }
         }
         else
         {
             if (x < -leftSigma_*nSigmasLeft_)
                 return 0.0;
             else
             {
                 const double dx = x/leftSigma_;
                 return norm_*exp(-dx*dx/2.0)/SQR2PI;
             }
         }
     }
 
     double BifurcatedGauss1D::unscaledCdf(const double x) const
     {
         if (x == 0.0)
             return leftCdfFrac_;
         else if (x > 0.0)
         {
             if (x > rightSigma_*nSigmasRight_)
                 return 1.0;
             else
             {
                 const double dx = x/rightSigma_;
                 const double cdfDelta = (1.0 + erf(dx/SQRT2))/2.0 - cdf0Right_;
                 return 1.0 + cdfDelta*(1.0 - leftCdfFrac_)/(cdf0Right_ - 0.5);
             }
         }
         else
         {
             if (x < -leftSigma_*nSigmasLeft_)
                 return 0.0;
             else
             {
                 const double dx = x/leftSigma_;
                 const double cdfDelta = erfc(-dx/SQRT2)/2.0 - cdf0Left_;
                 return cdfDelta*leftCdfFrac_/(0.5 - cdf0Left_);
             }
         }
     }
 
     double BifurcatedGauss1D::unscaledExceedance(const double x) const
     {
         return 1.0 - unscaledCdf(x);
     }
 
     double BifurcatedGauss1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::BifurcatedGauss1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return -nSigmasLeft_*leftSigma_;
         else if (r1 == 1.0)
             return nSigmasRight_*rightSigma_;
         else
         {
             if (r1 == leftCdfFrac_)
                 return 0.0;
             else if (r1 < leftCdfFrac_)
             {
                 // Map 0 into cdf0Left_ and leftCdfFrac_ into 0.5
                 const double arg = r1/leftCdfFrac_*(0.5-cdf0Left_) + cdf0Left_;
                 return leftSigma_*inverseGaussCdf(arg);
             }
             else
             {
                 // Map leftCdfFrac_ into 0.5 and 1.0 into cdf0Right_
                 const double d = (r1 - leftCdfFrac_)/(1.0 - leftCdfFrac_);
                 const double arg = 0.5 + d*(cdf0Right_ - 0.5);
                 return rightSigma_*inverseGaussCdf(arg);
             }
         }
     }
 
     bool BifurcatedGauss1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const BifurcatedGauss1D& r =
             static_cast<const BifurcatedGauss1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) &&
             leftSigma_ == r.leftSigma_ &&
             rightSigma_ == r.rightSigma_ &&
             nSigmasLeft_ == r.nSigmasLeft_ &&
             nSigmasRight_ == r.nSigmasRight_ &&
             norm_ == r.norm_ &&
             leftCdfFrac_ == r.leftCdfFrac_ &&
             cdf0Left_ == r.cdf0Left_ &&
             cdf0Right_ == r.cdf0Right_;
     }
 
     bool BifurcatedGauss1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, leftSigma_);
         gs::write_pod(os, rightSigma_);
         gs::write_pod(os, nSigmasLeft_);
         gs::write_pod(os, nSigmasRight_);
         gs::write_pod(os, norm_);
         gs::write_pod(os, leftCdfFrac_);
         gs::write_pod(os, cdf0Left_);
         gs::write_pod(os, cdf0Right_);
         return !os.fail();
     }
 
     BifurcatedGauss1D::BifurcatedGauss1D(const double location,
                                          const double scale)
         : AbsScalableDistribution1D(location, scale)
     {
     }
 
     BifurcatedGauss1D* BifurcatedGauss1D::read(
         const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<BifurcatedGauss1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             BifurcatedGauss1D* ptr = new BifurcatedGauss1D(location, scale);
             gs::read_pod(in, &ptr->leftSigma_);
             gs::read_pod(in, &ptr->rightSigma_);
             gs::read_pod(in, &ptr->nSigmasLeft_);
             gs::read_pod(in, &ptr->nSigmasRight_);
             gs::read_pod(in, &ptr->norm_);
             gs::read_pod(in, &ptr->leftCdfFrac_);
             gs::read_pod(in, &ptr->cdf0Left_);
             gs::read_pod(in, &ptr->cdf0Right_);
             if (!in.fail() &&
                 ptr->leftSigma_ >= 0.0 &&
                 ptr->rightSigma_ >= 0.0 &&
                 ptr->leftSigma_ + ptr->rightSigma_ > 0.0 &&
                 ptr->nSigmasLeft_ >= 0.0 &&
                 ptr->nSigmasRight_ >= 0.0 &&
                 ptr->nSigmasLeft_ + ptr->nSigmasRight_ > 0.0 &&
                 ptr->norm_ > 0.0 &&
                 ptr->leftCdfFrac_ >= 0.0 &&
                 ptr->cdf0Left_ >= 0.0 &&
                 ptr->cdf0Right_ > ptr->cdf0Left_)
                 return ptr;
             else
                 delete ptr;
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     Cauchy1D::Cauchy1D(const double location, const double scale,
                        const std::vector<double>& /* params */)
         : AbsScalableDistribution1D(location, scale),
           support_(sqrt(DBL_MAX/M_PI))
     {
     }
 
     Cauchy1D::Cauchy1D(const double location, const double scale)
         : AbsScalableDistribution1D(location, scale),
           support_(sqrt(DBL_MAX/M_PI))
     {
     }
 
     bool Cauchy1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Cauchy1D* Cauchy1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Cauchy1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Cauchy1D(location, scale);
     }
 
     double Cauchy1D::unscaledDensity(const double x) const
     {
         if (fabs(x) < support_)
             return 1.0/M_PI/(1.0 + x*x);
         else
             return 0.0;
     }
 
     double Cauchy1D::unscaledCdf(const double x) const
     {
         if (x < -support_)
             return 0.0;
         else if (x > support_)
             return 1.0;
         else
             return atan(x)/M_PI + 0.5;
     }
 
     double Cauchy1D::unscaledQuantile(const double x) const
     {
         if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error(
             "In npstat::Cauchy1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (x == 0.0)
             return -support_;
         else if (x == 1.0)
             return support_;
         else
             return tan(M_PI*(x - 0.5));
     }
 
     bool LogNormal::isEqual(const AbsDistribution1D& otherBase) const
     {
         const LogNormal& r = static_cast<const LogNormal&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                skew_ == r.skew_;
     }
 
     bool LogNormal::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, skew_);
         return !os.fail();
     }
 
     LogNormal* LogNormal::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<LogNormal>());
         current.ensureSameId(id);
 
         double location, scale, skew;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &skew);
             if (!in.fail())
                 return new LogNormal(location, scale, skew);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     void LogNormal::initialize()
     {
         logw_ = 0.0;
         s_ = 0.0;
         xi_ = 0.0;
         emgamovd_ = 0.0;
 
         if (skew_)
         {
             const double b1 = skew_*skew_;
             const double tmp = pow((2.0+b1+sqrt(b1*(4.0+b1)))/2.0, 1.0/3.0);
             const double w = tmp + 1.0/tmp - 1.0;
             logw_ = log(w);
             if (logw_ > 0.0)
             {
                 s_ = sqrt(logw_);
                 emgamovd_ = 1.0/sqrt(w*(w-1.0));
                 xi_ = -emgamovd_*sqrt(w);
             }
             else
             {
                 // This is not different from a Gaussian within
                 // the numerical precision of our calculations
                 logw_ = 0.0;
                 skew_ = 0.0;
             }
         }
     }
 
     LogNormal::LogNormal(const double mean, const double stdev,
                          const double skewness)
         : AbsScalableDistribution1D(mean, stdev),
           skew_(skewness)
     {
         initialize();
     }
 
     LogNormal::LogNormal(const double mean, const double stdev,
                          const std::vector<double>& params)
         : AbsScalableDistribution1D(mean, stdev),
           skew_(params[0])
     {
         initialize();
     }
 
     double LogNormal::unscaledDensity(const double x) const
     {
         if (skew_)
         {
             const double diff = skew_ > 0.0 ? x - xi_ : -x - xi_;
             if (diff <= 0.0)
                 return 0.0;
             else
             {
                 const double lg = log(diff/emgamovd_);
                 return exp(-lg*lg/2.0/logw_)/s_/SQR2PI/diff;
             }
         }
         else
         {
             // This is a Gaussian
             return exp(-x*x/2.0)/SQR2PI;
         }
     }
 
     double LogNormal::unscaledCdf(const double x) const
     {
         if (skew_)
         {
             const double diff = skew_ > 0.0 ? x - xi_ : -x - xi_;
             double posCdf = 0.0;
             if (diff > 0.0)
                 posCdf = (1.0 + erf(log(diff/emgamovd_)/s_/SQRT2))/2.0;
             return skew_ > 0.0 ? posCdf : 1.0 - posCdf;
         }
         else
             return (1.0 + erf(x/SQRT2))/2.0;
     }
 
     double LogNormal::unscaledExceedance(const double x) const
     {
         return 1.0 - unscaledCdf(x);
     }
 
     double LogNormal::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::LogNormal::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         const double g = inverseGaussCdf(skew_ >= 0.0 ? r1 : 1.0 - r1);
         if (skew_)
         {
             const double v = emgamovd_*exp(s_*g) + xi_;
             return skew_ > 0.0 ? v : -v;
         }
         else
             return g;
     }
 
     Moyal1D::Moyal1D(const double location, const double scale)
         : AbsScalableDistribution1D(location, scale),
           xmax_(-2.0*log(DBL_MIN*SQR2PI)),
           xmin_(-log(xmax_))
     {
     }
 
     Moyal1D::Moyal1D(const double location, const double scale,
                      const std::vector<double>& /* params */)
         : AbsScalableDistribution1D(location, scale),
           xmax_(-2.0*log(DBL_MIN*SQR2PI)),
           xmin_(-log(xmax_))
     {
     }
 
     double Moyal1D::unscaledDensity(const double x) const
     {
         if (x <= xmin_ || x >= xmax_)
             return 0.0;
         else
             return exp(-0.5*(x + exp(-x)))/SQR2PI;
     }
 
     double Moyal1D::unscaledCdf(const double x) const
     {
         if (x <= xmin_)
             return 0.0;
         else if (x >= xmax_)
             return 1.0;
         else
             return incompleteGammaC(0.5, 0.5*exp(-x));
     }
 
     double Moyal1D::unscaledExceedance(const double x) const
     {
         if (x <= xmin_)
             return 1.0;
         else if (x >= xmax_)
             return 0.0;
         else
             return incompleteGamma(0.5, 0.5*exp(-x));
     }
 
     double Moyal1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Moyal1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return xmin_;
         else if (r1 == 1.0)
             return xmax_;
         else
         {
             const double d = inverseIncompleteGammaC(0.5, r1);
             return -log(2.0*d);
         }
     }
 
     bool Moyal1D::write(std::ostream& os) const
     {
         return AbsScalableDistribution1D::write(os);
     }
 
     Moyal1D* Moyal1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Moyal1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (!AbsScalableDistribution1D::read(in, &location, &scale))
         {
             distributionReadError(in, classname());
             return 0;
         }
         return new Moyal1D(location, scale);
     }
 
     bool Pareto1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, c_);
         return !os.fail();
     }
 
     Pareto1D* Pareto1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Pareto1D>());
         current.ensureSameId(id);
 
         double location, scale, c;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &c);
             if (!in.fail())
                 return new Pareto1D(location, scale, c);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool Pareto1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Pareto1D& r = static_cast<const Pareto1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                c_ == r.c_;
     }
 
     void Pareto1D::initialize()
     {
         if (c_ <= 0.0) throw std::invalid_argument(
             "In npstat::Pareto1D::initialize: power parameter must be positive");
         if (c_ > 1.0)
             support_ = pow(1.0/DBL_MIN, 1.0/c_);
         else
             support_ = 1.0/DBL_MIN;
     }
 
     Pareto1D::Pareto1D(const double location, const double scale,
                        const double powerParam)
         : AbsScalableDistribution1D(location, scale),
           c_(powerParam)
     {
         initialize();
     }
 
     Pareto1D::Pareto1D(const double location, const double scale,
                        const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           c_(params[0])
     {
         initialize();
     }
 
     double Pareto1D::unscaledDensity(const double x) const
     {
         if (x < 1.0 || x > support_)
             return 0.0;
         else
             return c_/pow(x, c_ + 1.0);
     }
 
     double Pareto1D::unscaledCdf(const double x) const
     {
         if (x <= 1.0)
             return 0.0;
         else if (x >= support_)
             return 1.0;
         else
             return 1.0 - 1.0/pow(x, c_);
     }
 
     double Pareto1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Pareto1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 1.0;
         else if (r1 == 1.0)
             return support_;
         else
             return pow(1.0 - r1, -1.0/c_);
     }
 
     double Pareto1D::unscaledExceedance(const double x) const
     {
         if (x <= 1.0)
             return 1.0;
         else if (x >= support_)
             return 0.0;
         else
             return 1.0/pow(x, c_);
     }
 
     bool UniPareto1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, c_);
         return !os.fail();
     }
 
     UniPareto1D* UniPareto1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<UniPareto1D>());
         current.ensureSameId(id);
 
         double location, scale, c;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &c);
             if (!in.fail())
                 return new UniPareto1D(location, scale, c);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool UniPareto1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const UniPareto1D& r = static_cast<const UniPareto1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                c_ == r.c_;
     }
 
     void UniPareto1D::initialize()
     {
         if (c_ <= 0.0) throw std::invalid_argument(
             "In npstat::UniPareto1D::initialize: power parameter must be positive");
         if (c_ > 1.0)
             support_ = pow(1.0/DBL_MIN, 1.0/c_);
         else
             support_ = 1.0/DBL_MIN;
         amplitude_ = c_/(c_ + 1.0);
     }
 
     UniPareto1D::UniPareto1D(const double location, const double scale,
                              const double powerParam)
         : AbsScalableDistribution1D(location, scale),
           c_(powerParam)
     {
         initialize();
     }
 
     UniPareto1D::UniPareto1D(const double location, const double scale,
                              const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           c_(params[0])
     {
         initialize();
     }
 
     double UniPareto1D::unscaledDensity(const double x) const
     {
         if (x < 0.0 || x > support_)
             return 0.0;
         else if (x <= 1.0)
             return amplitude_;
         else
             return amplitude_/pow(x, c_ + 1.0);
     }
 
     double UniPareto1D::unscaledCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0;
         else if (x >= support_)
             return 1.0;
         else if (x <= 1.0)
             return x*amplitude_;
         else
             return amplitude_ + (1.0 - amplitude_)*(1.0 - 1.0/pow(x, c_));
     }
 
     double UniPareto1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::UniPareto1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 <= amplitude_)
             return r1/amplitude_;
         else if (r1 == 1.0)
             return support_;
         else
             return pow(1.0 - (r1 - amplitude_)/(1.0 - amplitude_), -1.0/c_);
     }
 
     double UniPareto1D::unscaledExceedance(const double x) const
     {
         if (x > 1.0)
             return (1.0 - amplitude_)/pow(x, c_);
         else
             return 1.0 - unscaledCdf(x);
     }
 
     bool Huber1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, tailWeight_);
         return !os.fail();
     }
 
     Huber1D* Huber1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Huber1D>());
         current.ensureSameId(id);
 
         double location, scale, tailWeight;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             gs::read_pod(in, &tailWeight);
             if (!in.fail())
                 return new Huber1D(location, scale, tailWeight);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool Huber1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Huber1D& r = static_cast<const Huber1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                tailWeight_ == r.tailWeight_;
     }
 
     void Huber1D::initialize()
     {
         if (!(tailWeight_ >= 0.0 && tailWeight_ < 1.0))
             throw std::invalid_argument(
                 "In npstat::Huber1D::initialize: "
                 "tail weight not inside [0, 1) interval");
 
         if (tailWeight_ == 0.0)
         {
             // Pure Gaussian
             a_ = DBL_MAX;
             normfactor_ = 1.0/sqrt(2.0*M_PI);
             support_ = inverseGaussCdf(1.0);
             cdf0_ = 0.0;
         }
         else
         {
             // Solve the equation for "a" by bisection
             const double eps = 2.0*DBL_EPSILON;
             double c = -2.0*log(tailWeight_);
             assert(c > 0.0);
             assert(weight(c) <= tailWeight_);
             double b = 0.0;
             while ((c - b)/(c + b) > eps)
             {
                 const double half = (c + b)/2.0;
                 if (weight(half) >= tailWeight_)
                     b = half;
                 else
                     c = half;
             }
             a_ = (c + b)/2.0;
             normfactor_ = 0.5/(exp(-a_*a_/2.0)/a_ + 
                                sqrt(M_PI/2.0)*erf(a_/SQRT2));
             support_ = a_/2.0 - log(DBL_MIN)/a_;
             cdf0_ = (1.0 + erf(-a_/SQRT2))/2.0;
         }
     }
 
     Huber1D::Huber1D(const double location, const double scale,
                      const double tailWeight)
         : AbsScalableDistribution1D(location, scale),
           tailWeight_(tailWeight)
     {
         initialize();
     }
 
     Huber1D::Huber1D(const double location, const double scale,
                      const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           tailWeight_(params[0])
     {
         initialize();
     }
 
     double Huber1D::unscaledDensity(const double x) const
     {
         const double absx = fabs(x);
         if (absx <= a_)
             return normfactor_*exp(-x*x/2.0);
         else
             return normfactor_*exp(a_*(a_/2.0 - absx));
     }
 
     double Huber1D::unscaledCdf(const double x) const
     {
         if (tailWeight_ == 0.0)
             return (1.0 + erf(x/SQRT2))/2.0;
         if (x < -a_)
             return normfactor_*exp((a_*(a_ + 2*x))/2.0)/a_;
         else if (x <= a_)
         {
             static const double sq1 = sqrt(M_PI/2.0);
             static const double sq2 = sqrt(2.0);
             return normfactor_*sq1*(erf(a_/sq2) + erf(x/sq2)) + tailWeight_/2;
         }
         else
             return 1.0 - normfactor_*exp((a_*(a_ - 2*x))/2.0)/a_;
     }
 
     double Huber1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Huber1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (tailWeight_ == 0.0)
             return inverseGaussCdf(r1);
         if (r1 == 0.0)
             return -support_;
         else if (r1 == 1.0)
             return support_;
         else if (r1 <= tailWeight_/2.0)
             return log(r1*a_/normfactor_)/a_ - 0.5*a_;
         else if (r1 < 1.0 - tailWeight_/2.0)
         {
             const double t = (r1 - tailWeight_/2.0)/normfactor_/SQR2PI + cdf0_;
             return inverseGaussCdf(t);
         }
         else
             return 0.5*a_ - log((1.0 - r1)*a_/normfactor_)/a_;
     }
 
     double Huber1D::weight(const double a) const
     {
         static const double sq1 = sqrt(M_PI/2.0);
         return 1.0/(1.0 + a*exp(a*a/2.0)*sq1*erf(a/SQRT2));
     }
 
     bool Tabulated1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, deg_);
         gs::write_pod_vector(os, table_);
         return !os.fail();
     }
 
     Tabulated1D* Tabulated1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<Tabulated1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             unsigned deg = 0;
             gs::read_pod(in, &deg);
             std::vector<double> table;
             gs::read_pod_vector(in, &table);
             if (!in.fail() && table.size())
                 return new Tabulated1D(location, scale,
                                        &table[0], table.size(), deg);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool Tabulated1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const Tabulated1D& r = static_cast<const Tabulated1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
                table_ == r.table_ && deg_ == r.deg_;
     }
 
     Tabulated1D::Tabulated1D(const double location, const double scale,
                              const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale)
     {
         const unsigned npara = params.size();
         if (npara)
             initialize(&params[0], npara, std::min(3U, npara-1U));
         else
             initialize(static_cast<double*>(0), 0U, 0U);
     }
 
     void Tabulated1D::normalize()
     {
         cdf_.clear();
         cdf_.reserve(len_);
         cdf_.push_back(0.0);
         long double sum = 0.0L;
         for (unsigned i=0; i<len_-1; ++i)
         {
             sum += intervalInteg(i);
             cdf_.push_back(static_cast<double>(sum));
         }
         double norm = cdf_[len_ - 1];
         assert(norm > 0.0);
         for (unsigned i=0; i<len_; ++i)
         {
             table_[i] /= norm;
             cdf_[i] /= norm;
         }
 
         exceed_.resize(len_);
         exceed_[len_ - 1U] = 0.0;
         sum = 0.0L;
         for (unsigned i=len_-1; i>0; --i)
         {
             sum += intervalInteg(i-1);
             exceed_[i-1] = static_cast<double>(sum);
         }
         norm = exceed_[0];
         for (unsigned i=0; i<len_; ++i)
             exceed_[i] /= norm;
     }
 
     double Tabulated1D::intervalInteg(const unsigned i) const
     {
         // The formula used here is exact for cubic polynomials
         static const double legendreRootOver2 = 0.5/sqrt(3.0);
         const double x0 = step_*(i + 0.5);
         const double v0 = unscaledDensity(x0 - legendreRootOver2*step_);
         const double v1 = unscaledDensity(x0 + legendreRootOver2*step_);
         return step_*(v0 + v1)/2.0;
     }
 
     double Tabulated1D::interpolate(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
         if (x == 0.0)
             return table_[0];
         if (x == 1.0)
             return table_[len_ - 1];
 
         unsigned idx = static_cast<unsigned>(x/step_);
         if (idx >= len_ - 1)
             idx = len_ - 2;
         const double dx = x/step_ - idx;
 
         switch (deg_)
         {
         case 0:
             if (dx < 0.5)
                 return table_[idx];
             else
                 return table_[idx + 1];
 
         case 1:
             return interpolate_linear(dx, table_[idx], table_[idx + 1]);
 
         case 2:
             if (idx == 0)
                 return interpolate_quadratic(dx, table_[idx], table_[idx + 1],
                                           table_[idx + 2]);
             else if (idx == len_ - 2)
                 return interpolate_quadratic(dx+1.0, table_[idx - 1],
                                           table_[idx], table_[idx + 1]);
             else
             {
                 const double v0 = interpolate_quadratic(
                     dx, table_[idx], table_[idx + 1], table_[idx + 2]);
                 const double v1 = interpolate_quadratic(
                     dx+1.0, table_[idx - 1], table_[idx], table_[idx + 1]);
                 return (v0 + v1)/2.0;
             }
 
         case 3:
             if (idx == 0)
                 return interpolate_cubic(dx, table_[idx], table_[idx + 1],
                                          table_[idx + 2], table_[idx + 3]);
             else if (idx == len_ - 2)
                 return interpolate_cubic(dx+2.0, table_[idx-2], table_[idx-1],
                                          table_[idx], table_[idx + 1]);
             else
                 return interpolate_cubic(dx+1.0, table_[idx - 1], table_[idx],
                                          table_[idx + 1], table_[idx + 2]);
 
         default:
             assert(0);
             return 0.0;
         }
     }
 
     double Tabulated1D::unscaledDensity(const double x) const
     {
         const double v = this->interpolate(x);
         if (v >= 0.0)
             return v;
         else
             return 0.0;
     }
 
     double Tabulated1D::unscaledCdf(double x) const
     {
         if (x <= 0.0)
             return 0.0;
         if (x >= 1.0)
             return 1.0;
 
         unsigned idx = static_cast<unsigned>(x/step_);
         if (idx >= len_ - 1)
             idx = len_ - 2;
 
         double v;
         switch (deg_)
         {
         case 0:
         {
             const double dx = x/step_ - idx;
             if (dx < 0.5)
                 v = table_[idx]*dx*step_;
             else
                 v = (table_[idx]*0.5 + table_[idx + 1]*(dx - 0.5))*step_;
         }
         break;
 
         default:
             v = interpIntegral(step_*idx, x);
         }
 
         return cdf_[idx] + v;
     }
 
     double Tabulated1D::unscaledExceedance(double x) const
     {
         if (x <= 0.0)
             return 1.0;
         if (x >= 1.0)
             return 0.0;
 
         unsigned idx = static_cast<unsigned>(x/step_);
         if (idx >= len_ - 1)
             idx = len_ - 2;
 
         double v;
         switch (deg_)
         {
         case 0:
         {
             const double dx = x/step_ - idx;
             if (dx < 0.5)
                 v = table_[idx]*dx*step_;
             else
                 v = (table_[idx]*0.5 + table_[idx + 1]*(dx - 0.5))*step_;
         }
         break;
 
         default:
             v = interpIntegral(step_*idx, x);
         }
 
         return exceed_[idx] - v;
     }
 
     double Tabulated1D::interpIntegral(const double a, const double b) const
     {
         static const double legendreRootOver2 = 0.5/sqrt(3.0);
         const double x0 = (b + a)/2.0;
         const double step = b - a;
         const double v0 = unscaledDensity(x0 - legendreRootOver2*step);
         const double v1 = unscaledDensity(x0 + legendreRootOver2*step);
         return step*(v0 + v1)/2.0;
     }
 
     double Tabulated1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::Tabulated1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         if (r1 == 1.0)
             return 1.0;
 
         unsigned idx = std::lower_bound(cdf_.begin(), cdf_.end(), r1) -
                                         cdf_.begin() - 1U;
         double xlo = step_*idx;
         const double dcdf = r1 - cdf_[idx];
         assert(dcdf > 0.0);
 
         switch (deg_)
         {
         case 0:
         {
             const double c1 = table_[idx]*0.5*step_;
             if (dcdf <= c1)
             {
                 assert(table_[idx] > 0.0);
                 return xlo + dcdf/table_[idx];
             }
             else
             {
                 assert(table_[idx+1] > 0.0);
                 return xlo + 0.5*step_ + (dcdf - c1)/table_[idx+1];
             }
         }
 
         case 1:
         {
             const double a = (table_[idx+1] - table_[idx])/step_/2.0;
             if (a == 0.0)
             {
                 assert(table_[idx] > 0.0);
                 return xlo + dcdf/table_[idx];
             }
             else
             {
                 double x1, x2;
                 const unsigned nroots = solveQuadratic(
                     table_[idx]/a, -dcdf/a, &x1, &x2);
                 if (nroots != 2U) throw std::runtime_error(
                     "In npstat::Tabulated1D::unscaledQuantile: "
                     "unexpected number of solutions");
                 if (fabs(x1 - 0.5*step_) < fabs(x2 - 0.5*step_))
                     return xlo + x1;
                 else
                     return xlo + x2;
             }
         }
 
         default:
         {
             double xhi = xlo + step_;
             const double eps = 2.0*DBL_EPSILON;
             while ((xhi - xlo)/(xhi + xlo) > eps)
             {
                 const double med = (xhi + xlo)/2.0;
                 if (unscaledCdf(med) >= r1)
                     xhi = med;
                 else
                     xlo = med;
             }
             return (xhi + xlo)/2.0;
         }
         }
     }
 
     bool BinnedDensity1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const double eps = 1.0e-12;
 
         const BinnedDensity1D& r = 
             static_cast<const BinnedDensity1D&>(otherBase);
         if (!AbsScalableDistribution1D::isEqual(r))
             return false;
         if (!(deg_ == r.deg_))
             return false;
         const unsigned long n = table_.size();
         if (!(n == r.table_.size()))
             return false;
         for (unsigned long i=0; i<n; ++i)
             if (fabs(table_[i] - r.table_[i])/
                 ((fabs(table_[i]) + fabs(r.table_[i]))/2.0 + 1.0) > eps)
                 return false;
         return true;
     }
 
     BinnedDensity1D::BinnedDensity1D(const double location, const double scale,
                                      const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale)
     {
         const unsigned npara = params.size();
         if (npara)
             initialize(&params[0], npara, std::min(1U, npara-1U));
         else
             initialize(static_cast<double*>(0), 0U, 0U);
     }
 
     void BinnedDensity1D::normalize()
     {
         cdf_.clear();
         cdf_.reserve(len_);
         long double sum = 0.0L;
 
         switch (deg_)
         {
         case 0U:
             for (unsigned i=0; i<len_; ++i)
             {
                 sum += table_[i];
                 cdf_.push_back(static_cast<double>(sum));
             }
             break;
 
         case 1U:
         {
             double oldval = 0.0;
             double* data = &table_[0];
             for (unsigned i=0; i<len_; ++i)
             {
                 sum += (data[i] + oldval)*0.5;
                 oldval = data[i];
                 cdf_.push_back(static_cast<double>(sum));
             }
             sum += oldval*0.5;
         }
         break;
 
         default:
             assert(0);
         }
         
         const double norm = static_cast<double>(sum);
         assert(norm > 0.0);
         const double integ = norm*step_;
         for (unsigned i=0; i<len_; ++i)
         {
             table_[i] /= integ;
             cdf_[i] /= norm;
         }
     }
 
     double BinnedDensity1D::unscaledDensity(const double x) const
     {
         double v = interpolate(x);
         if (v < 0.0)
             v = 0.0;
         return v;
     }
 
     double BinnedDensity1D::interpolate(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
 
         switch (deg_)
         {
         case 0:
         {
             unsigned idx = static_cast<unsigned>(x/step_);
             if (idx > len_ - 1)
                 idx = len_ - 1;
             return table_[idx];
         }
 
         case 1:
         {
             const double xs = x - step_/2.0;
             if (xs <= 0.0)
                 return table_[0];
             const unsigned idx = static_cast<unsigned>(xs/step_);
             if (idx > len_ - 2)
                 return table_[len_ - 1];
             const double dx = xs/step_ - idx;
             return interpolate_linear(dx, table_[idx], table_[idx + 1]);
         }
 
         default:
             assert(0);
             return 0.0;
         }
     }
 
     double BinnedDensity1D::unscaledCdf(double x) const
     {
         if (x <= 0.0)
             return 0.0;
         if (x >= 1.0)
             return 1.0;
 
         double v = 0.0;
         switch (deg_)
         {
         case 0:
         {            
             unsigned idx = static_cast<unsigned>(x/step_);
             if (idx > len_ - 1)
                 idx = len_ - 1;
             v = (idx ? cdf_[idx - 1] : 0.0) + table_[idx]*(x - idx*step_);
         }
         break;
 
         case 1:
         {
             const double xs = x - step_/2.0;
             if (xs <= 0.0)
                 v = table_[0]*x;
             else 
             {
                 const unsigned idx = static_cast<unsigned>(xs/step_);
                 if (idx > len_ - 2)
                     v = 1.0 - table_[len_ - 1]*(1.0 - x);
                 else
                 {
                     const double dx = xs - idx*step_;
                     const double slope = (table_[idx+1] - table_[idx])/step_;
                     v = cdf_[idx] + table_[idx]*dx + slope*dx*dx/2.0;
                 }
             }
         }
         break;
 
         default:
             assert(0);
             break;
         }
         if (v < 0.0)
             v = 0.0;
         else if (v > 1.0)
             v = 1.0;
         return v;
     }
 
     double BinnedDensity1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::BinnedDensity1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0 && deg_)
             return 0.0;
         if (r1 == 1.0 && deg_)
             return 1.0;
 
         double v = 0.0;
         switch (deg_)
         {
         case 0:
         {
             if (r1 == 0)
                 v = firstNonZeroBin_;
             else if (r1 == 1.0)
                 v = lastNonZeroBin_ + 1U;
             else if (r1 <= cdf_[0])
                 v = r1/cdf_[0];
             else
             {
                 double rem;
                 const unsigned bin = quantileBinFromCdf(&cdf_[0],len_,r1,&rem) + 1U;
                 assert(bin < len_);
                 v = bin + rem;
             }
         }
         break;
 
         case 1:
         {
             if (r1 <= cdf_[0])
                 v = 0.5*r1/cdf_[0];
             else if (r1 >= cdf_[len_ - 1])
                 v = (len_-0.5+0.5*(r1-cdf_[len_-1])/(1.0-cdf_[len_-1]));
             else
             {
                 const unsigned idx = std::lower_bound(cdf_.begin(),cdf_.end(),r1) -
                                      cdf_.begin() - 1U;
                 assert(idx < len_ - 1);
                 const double k = (table_[idx+1] - table_[idx])/step_;
                 const double y = r1 - cdf_[idx];
                 double x;
                 if (fabs(k) < 1.e-10*table_[idx])
                     x = y/table_[idx];
                 else
                 {
                     const double b = 2.0*table_[idx]/k;
                     const double c = -2.0*y/k;
                     double x1, x2;
                     if (solveQuadratic(b, c, &x1, &x2))
                     {
                         if (fabs(x1 - step_*0.5) < fabs(x2 - step_*0.5))
                             x = x1;
                         else
                             x = x2;
                     }
                     else
                     {
                         // This can happen due to various round-off problems.
                         // Assume that the quadratic equation determinant
                         // should have been 0 instead of negative.
                         x = -b/2.0;
                     }
                 }
                 if (x < 0.0)
                     x = 0.0;
                 else if (x > step_)
                     x = step_;
                 v = x/step_ + idx + 0.5;
             }
         }
         break;
 
         default:
             assert(0);
             return 0.0;
         }
         v *= step_;
         if (v < 0.0)
             v = 0.0;
         else if (v > 1.0)
             v = 1.0;
         return v;
     }
 
     bool BinnedDensity1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, deg_);
         gs::write_pod_vector(os, table_);
         return !os.fail();
     }
 
     BinnedDensity1D* BinnedDensity1D::read(const gs::ClassId& id,
                                            std::istream& in)
     {
         static const gs::ClassId current(
             gs::ClassId::makeId<BinnedDensity1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             unsigned deg = 0;
             gs::read_pod(in, &deg);
             std::vector<double> table;
             gs::read_pod_vector(in, &table);
             if (!in.fail() && table.size())
                 return new BinnedDensity1D(location, scale,
                                            &table[0], table.size(), deg);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool StudentsT1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod(os, nDoF_);
         return !os.fail();
     }
 
     StudentsT1D* StudentsT1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<StudentsT1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             double nDoF;
             gs::read_pod(in, &nDoF);
             if (!in.fail())
                 return new StudentsT1D(location, scale, nDoF);
         }
         distributionReadError(in, classname());
         return 0;
     }
 
     bool StudentsT1D::isEqual(const AbsDistribution1D& otherBase) const
     {
         const StudentsT1D& r = static_cast<const StudentsT1D&>(otherBase);
         return AbsScalableDistribution1D::isEqual(r) && 
             nDoF_ == r.nDoF_;
     }
 
     StudentsT1D::StudentsT1D(const double location, const double scale,
                              const double degreesOfFreedom)
         : AbsScalableDistribution1D(location, scale),
           nDoF_(degreesOfFreedom)
     {
         initialize();
     }
     
     StudentsT1D::StudentsT1D(const double location, const double scale,
                              const std::vector<double>& params)
         : AbsScalableDistribution1D(location, scale),
           nDoF_(params[0])
     {
         initialize();
     }
 
     void StudentsT1D::initialize()
     {
         if (nDoF_ <= 0.0) throw std::invalid_argument(
             "In npstat::StudentsT1D::initialize: invalid number "
             "of degrees of freedom");
         power_ = (nDoF_ + 1.0)/2.0;
         bignum_ = 0.0;
         normfactor_ = Gamma(power_)/Gamma(nDoF_/2.0)/sqrt(nDoF_)/SQRPI;
     }
 
     double StudentsT1D::unscaledDensity(const double x) const
     {
         return normfactor_*pow(1.0 + x*x/nDoF_, -power_);
     }
 
     double StudentsT1D::unscaledCdf(const double t) const
     {
         const double s = sqrt(t*t + nDoF_);
         return incompleteBeta(nDoF_/2.0, nDoF_/2.0, (t + s)/(2.0*s));
     }
 
     double StudentsT1D::unscaledExceedance(const double x) const
     {
         return 1.0 - unscaledCdf(x);
     }
 
     double StudentsT1D::unscaledQuantile(const double r1) const
     {
         if (r1 == 0.5)
             return 0.0;
         const double c = inverseIncompleteBeta(nDoF_/2.0, nDoF_/2.0, r1);
         const double tmp = 2.0*c - 1.0;
         const double a = tmp*tmp;
         const double denom = 1.0 - a;
         if (denom > 0.0)
         {
             const double sqroot = sqrt(nDoF_*a/denom);
             return r1 > 0.5 ? sqroot : -sqroot;
         }
         else
         {
             if (bignum_ == 0.0)
                 (const_cast<StudentsT1D*>(this))->bignum_ = effectiveSupport();
             return r1 > 0.5 ? bignum_ : -bignum_;
         }
     }
 
     double StudentsT1D::effectiveSupport() const
     {
         // Figure out at which (positive) values of the argument
         // the density becomes effectively indistinguishable from 0
         const double biglog = (log(normfactor_) - log(DBL_MIN) + 
                                power_*log(nDoF_))/(2.0*power_);
         if (biglog >= log(DBL_MAX))
             return DBL_MAX;
         else
             return exp(biglog);
     }
 }
Index: trunk/npstat/stat/Distributions1D.hh
===================================================================
--- trunk/npstat/stat/Distributions1D.hh	(revision 743)
+++ trunk/npstat/stat/Distributions1D.hh	(revision 744)
@@ -1,1067 +1,1067 @@
 #ifndef NPSTAT_DISTRIBUTIONS1D_HH_
 #define NPSTAT_DISTRIBUTIONS1D_HH_
 
 /*!
 // \file Distributions1D.hh
 //
 // \brief A number of useful 1-d continuous statistical distributions
 //
 // Author: I. Volobouev
 //
 // November 2009
 */
 
 #include "npstat/stat/Distribution1DFactory.hh"
 
 namespace npstat {
     /**
     // The uniform distribution is defined here by a constant density
     // equal to 1 between 0 and 1 and equal to 0 everywhere else
     */
     class Uniform1D : public AbsScalableDistribution1D
     {
     public:
         inline Uniform1D(const double location, const double scale)
             : AbsScalableDistribution1D(location, scale) {}
         inline virtual Uniform1D* clone() const {return new Uniform1D(*this);}
 
         inline virtual ~Uniform1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Uniform1D";}
         static inline unsigned version() {return 1;}
         static Uniform1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Uniform1D>;
 
         inline Uniform1D(const double location, const double scale,
                          const std::vector<double>& /* params */)
             : AbsScalableDistribution1D(location, scale) {}
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return 1.0 - unscaledCdf(x);}
     };
 
     /** Isosceles triangle distribution: 1 - |x| supported on [-1, 1] */
     class IsoscelesTriangle1D : public AbsScalableDistribution1D
     {
     public:
         inline IsoscelesTriangle1D(const double location, const double scale)
             : AbsScalableDistribution1D(location, scale) {}
         inline virtual IsoscelesTriangle1D* clone() const
             {return new IsoscelesTriangle1D(*this);}
 
         inline virtual ~IsoscelesTriangle1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname()
             {return "npstat::IsoscelesTriangle1D";}
         static inline unsigned version() {return 1;}
         static IsoscelesTriangle1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<IsoscelesTriangle1D>;
 
         inline IsoscelesTriangle1D(const double location, const double scale,
                                    const std::vector<double>& /* params */)
             : AbsScalableDistribution1D(location, scale) {}
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return 1.0 - unscaledCdf(x);}
     };
 
     /** Exponential distribution. "scale" is the decay time. */
     class Exponential1D : public AbsScalableDistribution1D
     {
     public:
         inline Exponential1D(const double location, const double scale)
             : AbsScalableDistribution1D(location, scale) {}
         inline virtual Exponential1D* clone() const
             {return new Exponential1D(*this);}
 
         inline virtual ~Exponential1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Exponential1D";}
         static inline unsigned version() {return 1;}
         static Exponential1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Exponential1D>;
         
         inline Exponential1D(const double location, const double scale,
                              const std::vector<double>& /* params */)
             : AbsScalableDistribution1D(location, scale) {}
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
     };
 
     /** Logistic distribution */
     class Logistic1D : public AbsScalableDistribution1D
     {
     public:
         inline Logistic1D(const double location, const double scale)
             : AbsScalableDistribution1D(location, scale) {}
         inline virtual Logistic1D* clone() const
             {return new Logistic1D(*this);}
 
         inline virtual ~Logistic1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Logistic1D";}
         static inline unsigned version() {return 1;}
         static Logistic1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Logistic1D>;
         
         inline Logistic1D(const double location, const double scale,
                           const std::vector<double>& /* params */)
             : AbsScalableDistribution1D(location, scale) {}
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
     };
 
     /**
     // A distribution whose density has a simple quadratic shape.
     // The support is from 0 to 1, and the coefficients "a" and "b"
     // are the coefficients for the Legendre polynomials of 1st
     // and 2nd degree translated to the support region. Note that
     // only those values of "a" and "b" that guarantee non-negativity
     // of the density are allowed, otherwise the code will generate
     // a run-time error.
     */
     class Quadratic1D : public AbsScalableDistribution1D
     {
     public:
         Quadratic1D(double location, double scale, double a, double b);
         inline virtual Quadratic1D* clone() const
             {return new Quadratic1D(*this);}
 
         inline virtual ~Quadratic1D() {}
 
         inline double a() const {return a_;}
         inline double b() const {return b_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Quadratic1D";}
-        static inline unsigned version() {return 2;}
+        static inline unsigned version() {return 3;}
         static Quadratic1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Quadratic1D>;
         
         Quadratic1D(double location, double scale,
                     const std::vector<double>& params);
         inline static int nParameters() {return 2;}
 
         void verifyNonNegative();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return 1.0 - unscaledCdf(x);}
 
         double a_;
         double b_;
     };
 
     /**
     // A distribution whose density logarithm has a simple quadratic
     // shape. The support is from 0 to 1, and the coefficients "a" and "b"
     // are the coefficients for the Legendre polynomials of 1st and 2nd
     // degree translated to the support region.
     */
     class LogQuadratic1D : public AbsScalableDistribution1D
     {
     public:
         LogQuadratic1D(double location, double scale, double a, double b);
         inline virtual LogQuadratic1D* clone() const
             {return new LogQuadratic1D(*this);}
 
         inline virtual ~LogQuadratic1D() {}
 
         inline double a() const {return a_;}
         inline double b() const {return b_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::LogQuadratic1D";}
-        static inline unsigned version() {return 2;}
+        static inline unsigned version() {return 3;}
         static LogQuadratic1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<LogQuadratic1D>;
 
         LogQuadratic1D(double location, double scale,
                        const std::vector<double>& params);
         inline static int nParameters() {return 2;}
 
         void normalize();
         long double quadInteg(long double x) const;
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return 1.0 - unscaledCdf(x);}
 
         long double ref_;
         long double range_;
         double a_;
         double b_;
         double k_;
         double s_;
         double norm_;
     };
 
     /** The Gaussian (or Normal) distribution */
     class Gauss1D : public AbsScalableDistribution1D
     {
     public:
         Gauss1D(double location, double scale);
         inline virtual Gauss1D* clone() const {return new Gauss1D(*this);}
 
         inline virtual ~Gauss1D() {}
 
         // Higher quality generator than the one provided by
         // the quantile function
         virtual unsigned random(AbsRandomGenerator& g,
                                 double* generatedRandom) const;
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Gauss1D";}
         static inline unsigned version() {return 1;}
         static Gauss1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Gauss1D>;
 
         Gauss1D(const double location, const double scale,
                 const std::vector<double>& params);
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         double xmin_;
         double xmax_;
     };
 
     /** Gaussian distribution truncated at some number of sigmas */
     class TruncatedGauss1D : public AbsScalableDistribution1D
     {
     public:
         TruncatedGauss1D(double location, double scale, double nsigma);
         inline virtual TruncatedGauss1D* clone() const
             {return new TruncatedGauss1D(*this);}
 
         inline virtual ~TruncatedGauss1D() {}
 
         inline double nsigma() const {return nsigma_;}
 
         // Higher quality generator than the one provided by
         // the quantile function
         virtual unsigned random(AbsRandomGenerator& g,
                                 double* generatedRandom) const;
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::TruncatedGauss1D";}
         static inline unsigned version() {return 1;}
         static TruncatedGauss1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<TruncatedGauss1D>;
 
         TruncatedGauss1D(double location, double scale,
                          const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return unscaledCdf(-x);}
 
         double nsigma_;
         double norm_;
         double cdf0_;
     };
 
     /**
     // Gaussian distribution on the [0, 1] interval, mirrored at the boundaries.
     // This is the Green's function of the diffusion equation on [0, 1]. The
     // interval can be shifted and scaled as for the uniform distribution.
     */
     class MirroredGauss1D : public AbsScalableDistribution1D
     {
     public:
         MirroredGauss1D(double location, double scale,
                         double meanOn0_1, double sigmaOn0_1);
 
         inline virtual MirroredGauss1D* clone() const
             {return new MirroredGauss1D(*this);}
 
         inline virtual ~MirroredGauss1D() {}
 
         inline double meanOn0_1() const {return mu0_;}
         inline double sigmaOn0_1() const {return sigma0_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::MirroredGauss1D";}
         static inline unsigned version() {return 1;}
         static MirroredGauss1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<MirroredGauss1D>;
 
         MirroredGauss1D(double location, double scale,
                         const std::vector<double>& params);
         inline static int nParameters() {return 2;}
 
         void validate();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
         long double ldCdf(double x) const;
 
         double mu0_;
         double sigma0_;
     };
 
     /**
     // Bifurcated Gaussian distribution. Different sigmas
     // can be used on the left and on the right, with constructor
     // parameter "leftSigmaFraction" specifying the ratio of
     // the left sigma to the sum of sigmas (this ratio must be
     // between 0 and 1). Different truncations in terms of the
     // number of sigmas can be used as well.
     */
     class BifurcatedGauss1D : public AbsScalableDistribution1D
     {
     public:
         BifurcatedGauss1D(double location, double scale,
                           double leftSigmaFraction,
                           double nSigmasLeft, double nSigmasRight);
         inline virtual BifurcatedGauss1D* clone() const
             {return new BifurcatedGauss1D(*this);}
 
         inline virtual ~BifurcatedGauss1D() {}
 
         inline double leftSigmaFraction() const
             {return leftSigma_/(leftSigma_ + rightSigma_);}
         inline double nSigmasLeft() const
             {return nSigmasLeft_;}
         inline double nSigmasRight() const
             {return nSigmasRight_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::BifurcatedGauss1D";}
         static inline unsigned version() {return 1;}
         static BifurcatedGauss1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         BifurcatedGauss1D(double location, double scale);
 
         friend class ScalableDistribution1DFactory<BifurcatedGauss1D>;
 
         BifurcatedGauss1D(double location, double scale,
                           const std::vector<double>& params);
         inline static int nParameters() {return 3;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         double leftSigma_;
         double rightSigma_;
         double nSigmasLeft_;
         double nSigmasRight_;
         double norm_;
         double leftCdfFrac_;
         double cdf0Left_;
         double cdf0Right_;
     };
 
     /** Symmetric beta distribution */
     class SymmetricBeta1D : public AbsScalableDistribution1D
     {
     public:
         SymmetricBeta1D(double location, double scale, double power);
 
         inline virtual SymmetricBeta1D* clone() const
             {return new SymmetricBeta1D(*this);}
 
         inline virtual ~SymmetricBeta1D() {}
 
         inline double power() const {return n_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::SymmetricBeta1D";}
         static inline unsigned version() {return 1;}
         static SymmetricBeta1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<SymmetricBeta1D>;
 
         SymmetricBeta1D(double location, double scale,
                         const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return unscaledCdf(-x);}
 
         double calculateNorm();
 
         double n_;
         double norm_;
         int intpow_;
     };
 
     /** Beta1D density is proportional to x^(apha-1) * (1-x)^(beta-1) */
     class Beta1D : public AbsScalableDistribution1D
     {
     public:
         Beta1D(double location, double scale, double alpha, double beta);
         inline virtual Beta1D* clone() const
             {return new Beta1D(*this);}
 
         inline virtual ~Beta1D() {}
 
         inline double alpha() const {return alpha_;}
         inline double beta() const {return beta_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Beta1D";}
         static inline unsigned version() {return 1;}
         static Beta1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Beta1D>;
 
         Beta1D(double location, double scale,
                const std::vector<double>& params);
         inline static int nParameters() {return 2;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledExceedance(double x) const;
         double unscaledQuantile(double x) const;
 
         double calculateNorm() const;
 
         double alpha_;
         double beta_;
         double norm_;
     };
 
     /** Shiftable gamma distribution */
     class Gamma1D : public AbsScalableDistribution1D
     {
     public:
         Gamma1D(double location, double scale, double alpha);
         inline virtual Gamma1D* clone() const
             {return new Gamma1D(*this);}
 
         inline virtual ~Gamma1D() {}
 
         inline double alpha() const {return alpha_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Gamma1D";}
         static inline unsigned version() {return 1;}
         static Gamma1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Gamma1D>;
 
         Gamma1D(double location, double scale,
                 const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledExceedance(double x) const;
         double unscaledQuantile(double x) const;
 
         void initialize();
 
         double alpha_;
         double norm_;
     };
 
     /**
     // Pareto distribution. Location parameter is location of 0, scale
     // parameter is the distance between 0 and the start of the density
     // (like the normal Pareto distribution location parameter).
     */
     class Pareto1D : public AbsScalableDistribution1D
     {
     public:
         Pareto1D(double location, double scale, double powerParameter);
         inline virtual Pareto1D* clone() const {return new Pareto1D(*this);}
 
         inline virtual ~Pareto1D() {}
 
         inline double powerParameter() const {return c_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Pareto1D";}
         static inline unsigned version() {return 1;}
         static Pareto1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Pareto1D>;
 
         Pareto1D(double location, double scale,
                  const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         double c_;
         double support_;
     };
 
     /**
     // Uniform distribution with Pareto tail attached to the right, where
     // the support of the uniform would normally end. Location parameter
     // is location of 0, scale parameter is the width of the uniform part
     // (like the normal Pareto distribution location parameter).
     */
     class UniPareto1D : public AbsScalableDistribution1D
     {
     public:
         UniPareto1D(double location, double scale, double powerParameter);
         inline virtual UniPareto1D* clone() const {return new UniPareto1D(*this);}
 
         inline virtual ~UniPareto1D() {}
 
         inline double powerParameter() const {return c_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::UniPareto1D";}
         static inline unsigned version() {return 1;}
         static UniPareto1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<UniPareto1D>;
 
         UniPareto1D(double location, double scale,
                     const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         double c_;
         double support_;
         double amplitude_;
     };
 
     /** "Huber" distribution */
     class Huber1D : public AbsScalableDistribution1D
     {
     public:
         Huber1D(double location, double scale, double tailWeight);
         inline virtual Huber1D* clone() const {return new Huber1D(*this);}
 
         inline virtual ~Huber1D() {}
 
         inline double tailWeight() const {return tailWeight_;}
         inline double tailStart() const {return a_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Huber1D";}
         static inline unsigned version() {return 1;}
         static Huber1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Huber1D>;
 
         Huber1D(double location, double scale,
                 const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return unscaledCdf(-x);}
 
         double weight(double a) const;
 
         double tailWeight_;
         double a_;
         double normfactor_;
         double support_;
         double cdf0_;
     };
 
     /** Cauchy (or Breit-Wigner) distribution */
     class Cauchy1D : public AbsScalableDistribution1D
     {
     public:
         Cauchy1D(double location, double scale);
         inline virtual Cauchy1D* clone() const {return new Cauchy1D(*this);}
 
         inline virtual ~Cauchy1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Cauchy1D";}
         static inline unsigned version() {return 1;}
         static Cauchy1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Cauchy1D>;
 
         Cauchy1D(const double location, const double scale,
                  const std::vector<double>& params);
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return unscaledCdf(-x);}
 
         double support_;
     };
 
     /**
     // Log-normal distribution represented by its mean, standard
     // deviation, and skewness. This representation is more useful
     // than other representations encountered in statistical literature.
     */
     class LogNormal : public AbsScalableDistribution1D
     {
     public:
         LogNormal(double mean, double stdev, double skewness);
         inline virtual LogNormal* clone() const {return new LogNormal(*this);}
 
         inline virtual ~LogNormal() {}
 
         inline double skewness() const {return skew_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::LogNormal";}
         static inline unsigned version() {return 1;}
         static LogNormal* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<LogNormal>;
 
         LogNormal(double location, double scale,
                   const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledExceedance(double x) const;
         double unscaledQuantile(double x) const;
 
         double skew_;
 
         double logw_;
         double s_;
         double xi_;
         double emgamovd_;
     };
 
     /**
     // Moyal Distribution (originally derived by Moyal as an approximation
     // to the Landau distribution)
     */
     class Moyal1D : public AbsScalableDistribution1D
     {
     public:
         Moyal1D(double location, double scale);
         inline virtual Moyal1D* clone() const {return new Moyal1D(*this);}
 
         inline virtual ~Moyal1D() {}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Moyal1D";}
         static inline unsigned version() {return 1;}
         static Moyal1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D& r) const
             {return AbsScalableDistribution1D::isEqual(r);}
 
     private:
         friend class ScalableDistribution1DFactory<Moyal1D>;
 
         Moyal1D(double location, double scale,
                 const std::vector<double>& params);
 
         inline static int nParameters() {return 0;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         double xmax_;
         double xmin_;
     };
 
     /** Student's t-distribution */
     class StudentsT1D : public AbsScalableDistribution1D
     {
     public:
         StudentsT1D(double location, double scale, double nDegreesOfFreedom);
         inline virtual StudentsT1D* clone() const
             {return new StudentsT1D(*this);}
 
         inline virtual ~StudentsT1D() {}
 
         inline double nDegreesOfFreedom() const {return nDoF_;}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::StudentsT1D";}
         static inline unsigned version() {return 1;}
         static StudentsT1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<StudentsT1D>;
 
         StudentsT1D(double location, double scale,
                     const std::vector<double>& params);
         inline static int nParameters() {return 1;}
 
         void initialize();
         double effectiveSupport() const;
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledExceedance(double x) const;
         double unscaledQuantile(double x) const;
 
         double nDoF_;
         double normfactor_;
         double power_;
         double bignum_;
     };
 
     /** Distribution defined by an interpolation table */
     class Tabulated1D : public AbsScalableDistribution1D
     {
     public:
         // The "data" array gives (unnormalized) density values at
         // equidistant intervals. data[0] is density at 0.0, and
         // data[dataLen-1] is density at 1.0. If "dataLen" is less
         // than 2, uniform distribution will be created. Internally,
         // the data is kept in double precision.
         //
         // "interpolationDegree" must be less than 4 and less than "dataLen".
         //
         template <typename Real>
         Tabulated1D(double location, double scale,
                     const Real* data, unsigned dataLen,
                     unsigned interpolationDegree);
 
         inline Tabulated1D(const double location, const double scale,
                            const std::vector<double>& table,
                            const unsigned interpolationDegree)
             : AbsScalableDistribution1D(location, scale)
         {
             const unsigned long sz = table.size();
             initialize(sz ? &table[0] : (double*)0, sz, interpolationDegree);
         }
 
         inline virtual Tabulated1D* clone() const
             {return new Tabulated1D(*this);}
 
         inline virtual ~Tabulated1D() {}
 
         inline unsigned interpolationDegree() const {return deg_;}
         inline unsigned tableLength() const {return len_;}
         inline const double* tableData() const {return &table_[0];}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::Tabulated1D";}
         static inline unsigned version() {return 1;}
         static Tabulated1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<Tabulated1D>;
 
         // The following constructor creates interpolator
         // of maximum degree possible
         Tabulated1D(double location, double scale,
                     const std::vector<double>& params);
         inline static int nParameters() {return -1;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
         double unscaledExceedance(double x) const;
 
         template <typename Real> void initialize(
             const Real* data, unsigned dataLen, unsigned interpolationDegree);
 
         void normalize();
         double interpolate(double x) const;
         double intervalInteg(unsigned intervalNumber) const;
         double interpIntegral(double x0, double x1) const;
 
         std::vector<double> table_;
         std::vector<double> cdf_;
         std::vector<double> exceed_;
         double step_;
         unsigned len_;
         unsigned deg_;
     };
 
     /**
     // Another interpolated distribution. For this one, we will assume
     // that the coordinates correspond to 1-d histogram bin centers.
     */
     class BinnedDensity1D : public AbsScalableDistribution1D
     {
     public:
         // The "data" array gives density values at equidistant intervals.
         // data[0] is density at 0.5/dataLen, and data[dataLen-1] is density
         // at 1.0 - 0.5/dataLen.
         //
         // "interpolationDegree" must be less than 2.
         //
         template <typename Real>
         BinnedDensity1D(double location, double scale,
                         const Real* data, unsigned dataLen,
                         unsigned interpolationDegree);
 
         inline BinnedDensity1D(const double location, const double scale,
                                const std::vector<double>& table,
                                const unsigned interpolationDegree)
             : AbsScalableDistribution1D(location, scale)
         {
             const unsigned long sz = table.size();
             initialize(sz ? &table[0] : (double*)0, sz, interpolationDegree);
         }
 
         inline virtual BinnedDensity1D* clone() const
             {return new BinnedDensity1D(*this);}
 
         inline virtual ~BinnedDensity1D() {}
 
         inline unsigned interpolationDegree() const {return deg_;}
         inline unsigned tableLength() const {return len_;}
         inline const double* tableData() const {return &table_[0];}
 
         // Methods needed for I/O
         virtual gs::ClassId classId() const {return gs::ClassId(*this);}
         virtual bool write(std::ostream& os) const;
 
         static inline const char* classname() {return "npstat::BinnedDensity1D";}
         static inline unsigned version() {return 1;}
         static BinnedDensity1D* read(const gs::ClassId& id, std::istream& in);
 
     protected:
         virtual bool isEqual(const AbsDistribution1D&) const;
 
     private:
         friend class ScalableDistribution1DFactory<BinnedDensity1D>;
 
         // The following constructor creates interpolator
         // of maximum degree possible
         BinnedDensity1D(double location, double scale,
                         const std::vector<double>& params);
         inline static int nParameters() {return -1;}
 
         double unscaledDensity(double x) const;
         double unscaledCdf(double x) const;
         double unscaledQuantile(double x) const;
 
         inline double unscaledExceedance(const double x) const
             {return 1.0 - unscaledCdf(x);}
 
         template <typename Real> void initialize(
             const Real* data, unsigned dataLen, unsigned interpolationDegree);
 
         void normalize();
         double interpolate(double x) const;
 
         std::vector<double> table_;
         std::vector<double> cdf_;
         double step_;
         unsigned len_;
         unsigned deg_;
         unsigned firstNonZeroBin_;
         unsigned lastNonZeroBin_;
     };
 }
 
 #include "npstat/stat/Distributions1D.icc"
 
 #endif // NPSTAT_DISTRIBUTIONS1D_HH_