discrete.hpp

Go to the documentation of this file.

#pragma once
 
#include <span>
 
#include <limits>
 
#include <randlib/math/constants.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/univariate.hpp>
 
namespace randlib
{
 
template <typename IntType>
class DiscreteDistribution : virtual public UnivariateDistribution<IntType>
{
    static_assert(std::is_integral_v<IntType> && std::is_signed_v<IntType>,
                  "Discrete distribution supports only signed integral types");
 
protected:
    DiscreteDistribution() {}
 
public:
    virtual ~DiscreteDistribution() {}
 
protected:
public:
    [[nodiscard]] virtual double P(const IntType &k) const { return std::exp(this->logP(k)); }
 
    [[nodiscard]] virtual double logP(const IntType &x) const = 0;
 
    void ProbabilityMassFunction(std::span<const IntType> x, std::span<double> y) const;
 
    void LogProbabilityMassFunction(std::span<const IntType> x, std::span<double> y) const;
 
    [[nodiscard]] IntType Mode() const override;
 
protected:
    IntType quantileImpl(double p, IntType initValue) const override;
    IntType quantileImpl(double p) const override;
    IntType quantileImpl1m(double p, IntType initValue) const override;
    IntType quantileImpl1m(double p) const override;
    long double ExpectedValue(const std::function<double(IntType)> &funPtr, IntType minPoint,
                              IntType maxPoint) const override;
    std::complex<double> CFImpl(double t) const override;
 
private:
    IntType searchQuantile(double p, bool upperTail) const;
 
public:
    [[nodiscard]] double Hazard(const IntType &x) const override;
 
    [[nodiscard]] double LikelihoodFunction(std::span<const IntType> sample) const override;
 
    [[nodiscard]] double LogLikelihoodFunction(std::span<const IntType> sample) const override;
 
    bool PearsonChiSquaredTest(std::span<const IntType> orderStatistic, double alpha, int lowerBoundary,
                               int upperBoundary, size_t numberOfEstimatedParameters = 0) const;
 
    bool PearsonChiSquaredTest(std::span<const IntType> orderStatistic, double alpha,
                               size_t numberOfEstimatedParameters = 0) const;
};
 
} // namespace randlib
 
#include <randlib/core/exponential_family.hpp>
#include <randlib/continuous/gamma.hpp>
 
namespace randlib
{
 
template <typename IntType>
void DiscreteDistribution<IntType>::ProbabilityMassFunction(std::span<const IntType> x, std::span<double> y) const
{
    const size_t size = std::min(x.size(), y.size());
    for (size_t i = 0; i != size; ++i)
        y[i] = this->P(x[i]);
}
 
template <typename IntType>
void DiscreteDistribution<IntType>::LogProbabilityMassFunction(std::span<const IntType> x, std::span<double> y) const
{
    const size_t size = std::min(x.size(), y.size());
    for (size_t i = 0; i != size; ++i)
        y[i] = this->logP(x[i]);
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::Mode() const
{
    IntType x = this->Median();
    double logProb = this->logP(x), newLogProb = this->logP(x + 1);
    if (logProb < newLogProb) {
        do {
            ++x;
            logProb = newLogProb;
            newLogProb = this->logP(x + 1);
        } while (logProb < newLogProb);
    } else {
        newLogProb = this->logP(x - 1);
        while (logProb < newLogProb) {
            --x;
            logProb = newLogProb;
            newLogProb = this->logP(x - 1);
        }
    }
    return x;
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::quantileImpl(double p, IntType initValue) const
{
    IntType down = initValue, up = down + 1;
    double fu = this->F(up), fd = this->F(down);
    while (fu < p) {
        fd = fu;
        fu = this->F(++up);
    }
    down = up - 1;
    while (fd > p) {
        fd = this->F(--down);
    }
    up = down + 1;
    return (fd < p) ? up : down;
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::quantileImpl(double p) const
{
    return searchQuantile(p, false);
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::quantileImpl1m(double p, IntType initValue) const
{
    IntType down = initValue, up = down + 1;
    double su = this->S(up), sd = this->S(down);
    while (su > p) {
        sd = su;
        su = this->S(++up);
    }
    down = up - 1;
    while (sd < p) {
        sd = this->S(--down);
    }
    up = down + 1;
 
    return (sd > p) ? up : down;
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::quantileImpl1m(double p) const
{
    return searchQuantile(p, true);
}
 
template <typename IntType>
IntType DiscreteDistribution<IntType>::searchQuantile(double p, bool upperTail) const
{
    static constexpr double slack = 8 * std::numeric_limits<double>::epsilon();
    auto reached = [this, p, upperTail](long long k) {
        return upperTail ? this->S(static_cast<IntType>(k)) <= p * (1 + slack)
                         : this->F(static_cast<IntType>(k)) >= p * (1 - slack);
    };
    const long long lo = this->MinValue(), hi = this->MaxValue();
    const bool loBounded = lo > static_cast<long long>(std::numeric_limits<IntType>::lowest());
 
    long long a, b;
    if (loBounded || !reached(0)) {
        a = loBounded ? lo - 1 : 0;
        long long step = 1;
        b = a + step;
        while (b < hi && !reached(b)) {
            a = b;
            step *= 2;
            b = (step > hi - a) ? hi : a + step;
        }
        b = std::min(b, hi);
    } else {
        b = 0;
        long long step = 1;
        a = b - step;
        while (a > lo && reached(a)) {
            b = a;
            step *= 2;
            a = (step > a - lo) ? lo : b - step;
        }
    }
    while (b - a > 1) {
        long long middle = a + (b - a) / 2;
        if (reached(middle))
            b = middle;
        else
            a = middle;
    }
    return static_cast<IntType>(b);
}
 
template <typename IntType>
long double DiscreteDistribution<IntType>::ExpectedValue(const std::function<double(IntType)> &funPtr, IntType minPoint,
                                                         IntType maxPoint) const
{
    Support suppType = this->SupportType();
    IntType k = minPoint, upperBoundary = maxPoint;
    if (suppType == Support::Finite || suppType == Support::RightSemiFinite) {
        k = std::max(k, this->MinValue());
    }
    if (suppType == Support::Finite || suppType == Support::LeftSemiFinite) {
        upperBoundary = std::min(upperBoundary, this->MaxValue());
    }
 
    double sum = 0;
    do {
        double addon = funPtr(k);
        if (addon != 0.0) {
            double prob = this->P(k);
            if (prob < math::detail::minimumTolerance)
                return sum;
            addon *= this->P(k);
            sum += addon;
        }
        ++k;
    } while (k <= upperBoundary);
    return sum;
}
 
template <typename IntType>
std::complex<double> DiscreteDistribution<IntType>::CFImpl(double t) const
{
    static constexpr long long MAX_TERMS = 1000000;
    static constexpr double EPSILON = 1e-17;
    const long long lo = this->MinValue(), hi = this->MaxValue(), start = this->Median();
    std::complex<long double> sum = 0;
    auto walk = [&](long long from, long long step, long long end, bool upper) {
        const std::complex<double> ratio = std::polar(1.0, step * t);
        const std::complex<double> w = ratio / (1.0 - ratio);
        const double damping = std::abs(1.0 - ratio);
        std::complex<double> zk = std::polar(1.0, t * static_cast<double>(from));
        double last[4] = {0, 0, 0, 0};
        long long k = from;
        for (long long terms = 0; (step > 0) ? k <= end : k >= end; k += step, ++terms) {
            double p = this->P(static_cast<IntType>(k));
            sum += std::complex<long double>(p * zk);
            zk *= ratio;
            if (terms % 64 == 63)
                zk = std::polar(1.0, t * static_cast<double>(k + step));
            last[3] = last[2];
            last[2] = last[1];
            last[1] = last[0];
            last[0] = p;
            if (p == 0.0 && terms > 16)
                return;
            if (terms < 16 || damping == 0.0)
                continue;
            double thirdDifference = std::fabs(last[0] - 3 * last[1] + 3 * last[2] - last[3]);
            double nextTerm = std::pow(std::abs(w), 3) * thirdDifference / damping;
            bool unbounded = (step > 0) ? end - k > MAX_TERMS : k - end > MAX_TERMS;
            bool decreasing = last[0] < last[1] && last[1] < last[2] && last[2] < last[3];
            bool converged = unbounded && decreasing && nextTerm < EPSILON * std::abs(sum);
            bool powerOfTwo = ((terms + 1) & terms) == 0;
            bool massGone = powerOfTwo &&
                            (upper ? this->S(static_cast<IntType>(k)) : this->F(static_cast<IntType>(k - 1))) < EPSILON;
            if (converged || massGone || terms >= MAX_TERMS)
                break;
        }
        if ((step > 0) ? k > end : k < end)
            return;
        k += step;
        double next[3];
        for (int j = 0; j < 3; ++j)
            next[j] = ((step > 0) ? k + j * step <= end : k + j * step >= end)
                          ? this->P(static_cast<IntType>(k + j * step))
                          : 0.0;
        std::complex<double> series = next[0] + w * (next[1] - next[0]) + w * w * (next[2] - 2 * next[1] + next[0]);
        sum += std::complex<long double>(zk * series / (1.0 - ratio));
    };
    walk(start, 1, hi, true);
    if (start > lo)
        walk(start - 1, -1, lo, false);
    return std::complex<double>(sum);
}
 
template <typename IntType>
double DiscreteDistribution<IntType>::Hazard(const IntType &x) const
{
    if (x < this->MinValue())
        return 0.0; 
    if (x > this->MaxValue())
        return NAN; 
    return this->P(x) / this->S(x);
}
 
template <typename IntType>
double DiscreteDistribution<IntType>::LikelihoodFunction(std::span<const IntType> sample) const
{
    long double res = 1.0;
    for (const IntType &var : sample)
        res *= this->P(var);
    return res;
}
 
template <typename IntType>
double DiscreteDistribution<IntType>::LogLikelihoodFunction(std::span<const IntType> sample) const
{
    long double res = 0.0;
    for (const IntType &var : sample)
        res += this->logP(var);
    return res;
}
 
template <typename IntType>
bool DiscreteDistribution<IntType>::PearsonChiSquaredTest(std::span<const IntType> orderStatistic, double alpha,
                                                          int lowerBoundary, int upperBoundary,
                                                          size_t numberOfEstimatedParameters) const
{
    const size_t n = orderStatistic.size();
 
    if (n == 0)
        throw std::invalid_argument("Order statistic should not be empty");
    if (!(alpha > 0 && alpha < 1))
        throw std::invalid_argument("Significance level should be positive and smaller than 1");
    if (lowerBoundary >= upperBoundary)
        throw std::invalid_argument("Lower boundary should be smaller than upper one");
    for (size_t j = 1; j != n; ++j) {
        if (orderStatistic[j] < orderStatistic[j - 1])
            throw std::invalid_argument("Order statistic should be sorted in ascending order");
    }
    if (orderStatistic[0] < this->MinValue())
        throw std::invalid_argument(
            "Some elements in the sample are too small to belong to this distribution, they should be greater than " +
            this->toStringWithPrecision(this->MinValue()));
    if (orderStatistic[n - 1] > this->MaxValue())
        throw std::invalid_argument(
            "Some elements in the sample are too large to belong to this distribution, they should be smaller than " +
            this->toStringWithPrecision(this->MaxValue()));
 
    const long long minValue = this->MinValue(), maxValue = this->MaxValue();
    const long long lower = (lowerBoundary > minValue) ? lowerBoundary : orderStatistic[0] - 1LL;
    const long long upper = (upperBoundary < maxValue) ? upperBoundary : orderStatistic[n - 1];
    const double nInv = 1.0 / n;
    double sum = 0.0;
    size_t k = 0;
    auto addCell = [&sum, &k, nInv](size_t count, double expectedProb) {
        if (!(expectedProb > 0)) {
            if (count > 0)
                sum = INFINITY;
            return;
        }
        double addon = nInv * count - expectedProb;
        sum += addon * addon / expectedProb;
        ++k;
    };
 
    size_t i = 0;
    if (lower >= minValue) {
        i = std::upper_bound(orderStatistic.begin(), orderStatistic.end(), lower) - orderStatistic.begin();
        addCell(i, this->F(static_cast<IntType>(lower)));
    }
    double observedMass = 0;
    long long distinct = 0;
    while (i < n && orderStatistic[i] <= upper) {
        const IntType x = orderStatistic[i];
        size_t count = 0;
        for (; i < n && orderStatistic[i] == x; ++i)
            ++count;
        const double probability = this->P(x);
        observedMass += probability;
        ++distinct;
        addCell(count, probability);
    }
    const long long unobserved = upper - lower - distinct;
    if (unobserved > 0) {
        const double cellMass =
            this->F(static_cast<IntType>(upper)) - ((lower >= minValue) ? this->F(static_cast<IntType>(lower)) : 0.0);
        sum += std::max(cellMass - observedMass, 0.0);
        k += static_cast<size_t>(unobserved);
    }
    if (upper < maxValue)
        addCell(n - i, this->S(static_cast<IntType>(upper)));
 
    if (k <= numberOfEstimatedParameters + 1) {
        throw std::invalid_argument("Sample is too small, number of groups (" + this->toStringWithPrecision(k) +
                                    ") should be bigger than number of estimated parameters plus 1 (" +
                                    this->toStringWithPrecision(numberOfEstimatedParameters + 1) + ")");
    }
    double statistic = n * sum;
    ChiSquaredRand X(k - 1 - numberOfEstimatedParameters);
    double q = X.Quantile1m(alpha);
    return statistic <= q;
}
 
template <typename IntType>
bool DiscreteDistribution<IntType>::PearsonChiSquaredTest(std::span<const IntType> orderStatistic, double alpha,
                                                          size_t numberOfEstimatedParameters) const
{
    return PearsonChiSquaredTest(orderStatistic, alpha, this->MinValue(), this->MaxValue(),
                                 numberOfEstimatedParameters);
}
 
} // namespace randlib