libcruft-util/test/rand/generator/normal.cpp

#include "maths.hpp"
#include "rand/distribution/normal.hpp"
#include "tap.hpp"

#include <iostream>
#include <random>


///////////////////////////////////////////////////////////////////////////////
/// Probability density function for a normal distribution with specified
/// mean and stddev at point `x`.
static
float pdf (float x, float mean, float stddev)
{
    float const power = cruft::pow2 ((x - mean) / stddev) / -2;
    float const scale = 1.f / (stddev * std::sqrt (2.f * cruft::pi<float>));
    return scale * std::exp (power);
}


///----------------------------------------------------------------------------
/// Calculate the maximum difference between a histogram and a PDF for a
/// normal distribution with a number of buckets.
static
float max_histogram_error (int buckets)
{
    // These constants weren't rigorously selected. Eyeballing the generated
    // values suggested they had some level of precision and didn't explode
    // the test's runtime.
    int const BUCKETS = buckets;
    int const ITERATIONS = BUCKETS * 10'000;
    float const MEAN   = BUCKETS / 2.f;
    float const STDDEV = BUCKETS * .15f;

    // Use _our_ normal distribution, not the stdlib one.
    cruft::rand::distribution::normal<float> g (MEAN, STDDEV);
    // We use a stdlib generator with reasonable quality just so we're not
    // testing both our generator and our distributions simultaneously.
    std::mt19937_64 rand;

    std::vector<int> counts (BUCKETS, 0);
    for (int i = 0; i < ITERATIONS; ++i) {
        auto const val = g (rand);
        if (val >= BUCKETS || val < 0)
            continue;
        counts[int (val)]++;
    }

    float max_diff = 0.f;
    for (int i = 0; i < BUCKETS; ++i) {
        float expected = pdf (i, MEAN, STDDEV);
        float actual   = counts[i] / float (ITERATIONS);
        float diff     = std::abs (expected - actual);

        max_diff = std::max (max_diff, diff);
    }

    return max_diff;
}


///////////////////////////////////////////////////////////////////////////////
int main (void)
{
    cruft::TAP::logger tap;
    tap.expect_lt (
        max_histogram_error (500),
        1.e-4f,
        "normal distribution histogram maximum relative error"
    );
    return tap.status ();
}