Page MenuHomeHEPForge
Files F8310399

No OneTemporary
Actions

View Options

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_
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);
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:22 PM (8 h, 34 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023781
Default Alt Text
(201 KB)

Event Timeline