kolmogorov_smirnov.hpp

Go to the documentation of this file.

#pragma once
 
#include <randlib/math/constants.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/continuous.hpp>
 
namespace randlib
{
 
template <typename RealType>
class KolmogorovSmirnovRand final : public ContinuousDistribution<RealType>
{
    RealType variate(detail::EngineRef g) const override { return (*this)(g); }
 
public:
    KolmogorovSmirnovRand();
    [[nodiscard]] std::string Name() const override;
    [[nodiscard]] Support SupportType() const override { return Support::RightSemiFinite; }
    [[nodiscard]] RealType MinValue() const override { return 0; }
    [[nodiscard]] RealType MaxValue() const override { return INFINITY; }
 
private:
    static double L(RealType x);
    static double K(RealType x);
 
public:
    [[nodiscard]] double f(const RealType &x) const override;
    [[nodiscard]] double logf(const RealType &x) const override;
    [[nodiscard]] double F(const RealType &x) const override;
    [[nodiscard]] double S(const RealType &x) const override;
    double logF(const RealType &x) const;
    double logS(const RealType &x) const;
 
private:
    template <std::uniform_random_bit_generator G>
    RealType truncatedGammaVariate(G &g) const;
    template <std::uniform_random_bit_generator G>
    RealType variateForTheLeftMostInterval(G &g) const;
    template <std::uniform_random_bit_generator G>
    RealType variateForTheRightMostInterval(G &g) const;
 
public:
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] RealType operator()(G &g) const;
 
    [[nodiscard]] long double Mean() const override;
    [[nodiscard]] long double Variance() const override;
    [[nodiscard]] RealType Mode() const override;
    [[nodiscard]] RealType Median() const override;
 
private:
    RealType quantileImpl(double p) const override;
    RealType quantileImpl1m(double p) const override;
};
 
} // namespace randlib
 
#include <randlib/continuous/exponential.hpp>
#include <randlib/continuous/uniform.hpp>
 
namespace randlib
{
 
template <typename RealType>
KolmogorovSmirnovRand<RealType>::KolmogorovSmirnovRand()
{
}
 
template <typename RealType>
std::string KolmogorovSmirnovRand<RealType>::Name() const
{
    return "Kolmogorov-Smirnov";
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::L(RealType x)
{
    if (x <= 0.0)
        return 0.0;
    double sum = 0.0, addon = 0.0;
    int k = 1;
    double aux = randlib::consts::pi_sq * 0.125 / (x * x);
    do {
        int temp = (2 * k - 1);
        temp *= temp;
        addon = std::exp(-temp * aux);
        sum += addon;
        ++k;
    } while (addon > math::detail::minimumTolerance * sum);
    return randlib::consts::sqrt2pi * sum / x;
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::K(RealType x)
{
    if (x <= 0.0)
        return 1.0;
    double sum = 0.0, addon = 0.0;
    int k = 1;
    double xSq = x * x;
    do {
        int temp = 2 * k * k;
        addon = std::exp(-temp * xSq);
        sum += (k & 1) ? addon : -addon;
        ++k;
    } while (addon > math::detail::minimumTolerance * sum);
    return 2 * sum;
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::f(const RealType &x) const
{
    if (x <= 0.0)
        return 0.0;
    double sum = 0.0, addon = 0.0;
    int k = 1;
    double xSq = x * x;
    if (x < 1.0) {
        double aux = 0.125 / xSq;
        do {
            double temp = randlib::consts::pi * (2 * k - 1);
            temp *= temp;
            addon = temp - 4 * xSq;
            addon *= std::exp(-temp * aux);
            sum += addon;
            ++k;
        } while (addon > math::detail::minimumTolerance * sum);
        return randlib::consts::sqrt2pi * sum * 0.25 / std::pow(x, 4);
    }
    do {
        int temp = k * k;
        addon = std::exp(-2 * temp * xSq);
        addon *= temp;
        sum += (k & 1) ? addon : -addon;
        ++k;
    } while (addon > math::detail::minimumTolerance * sum);
    return 8 * sum * x;
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::logf(const RealType &x) const
{
    return std::log(f(x));
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::F(const RealType &x) const
{
    return (x > 1.0) ? 1.0 - K(x) : L(x);
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::S(const RealType &x) const
{
    return (x > 1.0) ? K(x) : 1.0 - L(x);
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::logF(const RealType &x) const
{
    return (x > 1.0) ? std::log1pl(-K(x)) : std::log(L(x));
}
 
template <typename RealType>
double KolmogorovSmirnovRand<RealType>::logS(const RealType &x) const
{
    return (x > 1.0) ? std::log(K(x)) : std::log1pl(-L(x));
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType KolmogorovSmirnovRand<RealType>::truncatedGammaVariate(G &g) const
{
    static constexpr long double tp = 2.193245422464302l; 
    static constexpr long double rate = 1.2952909208355123l;
    size_t iter = 0;
    do {
        RealType E0 = rate * detail::standardExponential<RealType>(g);
        RealType E1 = 2 * detail::standardExponential<RealType>(g);
        RealType gammaValue = tp + E0;
        if (E0 * E0 <= tp * E1 * (gammaValue + tp))
            return gammaValue;
        RealType Wp = E0 / tp;
        if (Wp - std::log1pl(Wp) <= E1)
            return gammaValue;
    } while (++iter <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Kolmogorov-Smirnov distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType KolmogorovSmirnovRand<RealType>::variateForTheLeftMostInterval(G &g) const
{
    size_t iter1 = 0;
    do {
        RealType gammaValue = truncatedGammaVariate(g);
        RealType X = randlib::consts::pi / std::sqrt(8 * gammaValue);
        RealType W = 0.0;
        RealType Z = 0.5 / gammaValue;
        size_t n = 1, iter2 = 0;
        RealType Q = 1.0;
        RealType U = detail::uniform01<RealType>(g);
        while (U >= W && ++iter2 <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION) {
            W += Z * Q;
            if (U >= W)
                return X;
            n += 2;
            int nSq = n * n;
            Q = std::exp(gammaValue - gammaValue * nSq);
            W -= nSq * Q;
        }
    } while (++iter1 <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Kolmogorov-Smirnov distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType KolmogorovSmirnovRand<RealType>::variateForTheRightMostInterval(G &g) const
{
    static constexpr double tSq = 0.5625; 
    size_t iter1 = 0;
    do {
        RealType E = detail::standardExponential<RealType>(g);
        RealType U = detail::uniform01<RealType>(g);
        RealType X = std::sqrt(tSq + 0.5 * E);
        RealType W = 0.0;
        size_t n = 1, iter2 = 0;
        RealType Z = -2 * X * X;
        while (U > W && ++iter2 <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION) {
            ++n;
            int nSq = n * n;
            W += nSq * std::exp(Z * (nSq - 1));
            if (U >= W)
                return X;
            ++n;
            nSq = n * n;
            W -= nSq * std::exp(Z * (nSq - 1));
        }
    } while (++iter1 <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Kolmogorov-Smirnov distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType KolmogorovSmirnovRand<RealType>::operator()(G &g) const
{
    bool isLeft = detail::uniform01<RealType>(g) < 0.3728329582237386; 
    return isLeft ? variateForTheLeftMostInterval(g) : variateForTheRightMostInterval(g);
}
 
template <typename RealType>
long double KolmogorovSmirnovRand<RealType>::Mean() const
{
    return randlib::consts::sqrtpi * randlib::consts::inv_sqrt2 * randlib::consts::ln2;
}
 
template <typename RealType>
long double KolmogorovSmirnovRand<RealType>::Variance() const
{
    long double mean = Mean();
    return randlib::consts::pi_sq / 12 - mean * mean;
}
 
template <typename RealType>
RealType KolmogorovSmirnovRand<RealType>::Mode() const
{
    return 0.735467812776958l;
}
 
template <typename RealType>
RealType KolmogorovSmirnovRand<RealType>::Median() const
{
    return 0.82757355518990761l;
}
 
template <typename RealType>
RealType KolmogorovSmirnovRand<RealType>::quantileImpl(double p) const
{
    RealType guess = std::sqrt(-0.5 * (std::log1pl(-p) - randlib::consts::ln2));
    return ContinuousDistribution<RealType>::quantileImpl(p, guess);
}
 
template <typename RealType>
RealType KolmogorovSmirnovRand<RealType>::quantileImpl1m(double p) const
{
    RealType guess = std::sqrt(-0.5 * std::log(0.5 * p));
    return ContinuousDistribution<RealType>::quantileImpl1m(p, guess);
}
 
} // namespace randlib