Page MenuHomeHEPForge
Files F19252116

test_Distributions1D.cc
No OneTemporary
Actions

Size
16 KB
Referenced Files
None
Subscribers
None

test_Distributions1D.cc
View Options

#include <complex>
#include <iostream>
#include "UnitTest++.h"
#include "test_utils.hh"
#include "npstat/stat/Distributions1D.hh"
#include "npstat/stat/GaussianMixture1D.hh"
#include "npstat/stat/JohnsonCurves.hh"
#include "npstat/stat/TruncatedDistribution1D.hh"
#include "npstat/stat/LeftCensoredDistribution.hh"
#include "npstat/stat/RightCensoredDistribution.hh"
#include "npstat/stat/QuantileTable1D.hh"
#include "npstat/stat/arrayStats.hh"
#include "npstat/nm/EquidistantSequence.hh"
using namespace npstat;
using namespace std;
namespace {
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 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);
}
}
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);
}
}
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);
}
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_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_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_JohnsonSu)
{
JohnsonSu jsu1(0.0, 1.0, 1.0, 10.0);
standard_test(jsu1, 6, 1.e-8);
JohnsonSu jsu2(0.0, 1.0, -1.0, 10.0);
standard_test(jsu2, 6, 1.e-8);
}
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(Tabulated1D)
{
const double data[] = {0, 1, 2, 3, 3, 0};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<4; ++deg)
{
Tabulated1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-8);
}
}
TEST(QuantileTable1D)
{
const double data[] = {0.2, 0.25, 0.4, 0.5, 0.55};
const double location = 1.0;
const double width = 2.5;
QuantileTable1D qt(location, width, data, sizeof(data)/sizeof(data[0]));
invcdf_test(qt, 1.e-12);
standard_test(qt, 3.5, 5.0e-3);
}
TEST(BinnedDensity1D_1)
{
const double data[] = {0, 1, 3, 3, 3, 0, 8};
const double location = 1.7;
const double width = 5.9;
for (unsigned deg = 0; deg<2; ++deg)
{
BinnedDensity1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), deg);
invcdf_test(t1, 1.e-8);
double in = simpson_integral(DensityFunctor1D(t1), location,
location+width);
CHECK_CLOSE(1.0, in, 1.e-3);
}
}
TEST(BinnedDensity1D_2)
{
const double data[] = {0, 0, 1, 11, 3, 7, 0, 8, 0, 0};
const double location = 0.0;
const double width = 1.0;
const unsigned nquant = 1000;
EquidistantInLinearSpace qs(0.0, 1.0, nquant);
std::vector<double> quantiles(nquant);
arrayQuantiles1D(data, sizeof(data)/sizeof(data[0]), 0.0, 1.0,
&qs[0], &quantiles[0], nquant);
BinnedDensity1D t1(location, width,
data, sizeof(data)/sizeof(data[0]), 0);
for (unsigned i=0; i<nquant; ++i)
CHECK_CLOSE(quantiles[i], t1.quantile(qs[i]), 1.0e-14);
}
TEST(Quadratic1D_1)
{
Quadratic1D q1(-1.0, 2.0, 0.5, 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, 0.2);
standard_test(q1, 1.0, 1.e-10);
invcdf_test(q1, 1.e-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);
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);
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);
}
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(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(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-10;
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);
}
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:14 AM (4 h, 16 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6545833
Default Alt Text
test_Distributions1D.cc (16 KB)

Event Timeline

HEPForge.org · .