negative_hypergeometric.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/continuous/beta.hpp>
#include <randlib/core/discrete.hpp>
 
namespace randlib
{
 
template <typename IntType>
class NegativeHyperGeometricRand final : public DiscreteDistribution<IntType>
{
    IntType variate(detail::EngineRef g) const override { return (*this)(g); }
 
    IntType N = 1; 
    IntType M = 1; 
    IntType m = 1; 
    double pmfCoef = 1; 
    BetaRand<double> successFraction{}; 
 
public:
    NegativeHyperGeometricRand(IntType totalSize = 1, IntType totalSuccessesNum = 1, IntType limitSuccessesNum = 1);
    [[nodiscard]] std::string Name() const override;
    [[nodiscard]] Support SupportType() const override { return Support::Finite; }
    [[nodiscard]] IntType MinValue() const override { return 0; }
    [[nodiscard]] IntType MaxValue() const override { return N - M; }
 
    void SetParameters(IntType totalSize, IntType totalSuccessesNum, IntType limitSuccessesNum);
    [[nodiscard]] inline IntType GetTotalSize() { return N; }
    [[nodiscard]] inline IntType GetTotalSuccessesNum() { return M; }
    [[nodiscard]] inline IntType GetLimitSuccessesNum() { return m; }
 
    [[nodiscard]] double P(const IntType &k) const override;
    [[nodiscard]] double logP(const IntType &k) const override;
    [[nodiscard]] double F(const IntType &k) const override;
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] IntType operator()(G &g) const;
 
    [[nodiscard]] long double Mean() const override;
    [[nodiscard]] long double Variance() const override;
 
private:
    static constexpr double SMALL_URN = 16;
    template <std::uniform_random_bit_generator G>
    IntType variateByUrn(G &g) const;
};
 
} // namespace randlib
 
#include <randlib/discrete/bernoulli.hpp>
#include <randlib/discrete/binomial.hpp>
 
namespace randlib
{
 
template <typename IntType>
NegativeHyperGeometricRand<IntType>::NegativeHyperGeometricRand(IntType totalSize, IntType totalSuccessesNum,
                                                                IntType limitSuccessesNum)
{
    SetParameters(totalSize, totalSuccessesNum, limitSuccessesNum);
}
 
template <typename IntType>
std::string NegativeHyperGeometricRand<IntType>::Name() const
{
    return "Negative hypergeometric(" + this->toStringWithPrecision(N) + ", " + this->toStringWithPrecision(M) + ", " +
           this->toStringWithPrecision(m) + ")";
}
 
template <typename IntType>
void NegativeHyperGeometricRand<IntType>::SetParameters(IntType totalSize, IntType totalSuccessesNum,
                                                        IntType limitSuccessesNum)
{
    if (!(totalSize > 0 && totalSuccessesNum > 0 && limitSuccessesNum > 0))
        throw std::invalid_argument("Negative-HyperGeometric distribution: all parameters should be positive");
    if (!(totalSuccessesNum <= totalSize))
        throw std::invalid_argument(
            "Negative-HyperGeometric distribution: total size shouldn't be smaller than total successes number");
    if (!(limitSuccessesNum <= totalSuccessesNum))
        throw std::invalid_argument("Negative-HyperGeometric distribution: total successes number shouldn't be smaller "
                                    "than limit successes number");
 
    N = totalSize;
    M = totalSuccessesNum;
    m = limitSuccessesNum;
 
    successFraction.SetShapes(m, M - m + 1);
    pmfCoef = math::lfact(M);
    pmfCoef += math::lfact(N - M);
    pmfCoef -= math::lfact(m - 1);
    pmfCoef -= math::lfact(M - m);
    pmfCoef -= math::lfact(N);
}
 
template <typename IntType>
double NegativeHyperGeometricRand<IntType>::P(const IntType &k) const
{
    return (k < MinValue() || k > MaxValue()) ? 0.0 : std::exp(logP(k));
}
 
template <typename IntType>
double NegativeHyperGeometricRand<IntType>::logP(const IntType &k) const
{
    if (k < MinValue() || k > MaxValue())
        return -INFINITY;
    double p = math::lfact(k + m - 1);
    p += math::lfact(N - m - k);
    p -= math::lfact(k);
    p -= math::lfact(N - M - k);
    return p + pmfCoef;
}
 
template <typename IntType>
double NegativeHyperGeometricRand<IntType>::F(const IntType &k) const
{
    // relation with hypergeometric distribution can be used here instead
    if (k < MinValue())
        return 0.0;
    IntType maxVal = MaxValue();
    if (k >= maxVal)
        return 1.0;
    if (k <= 0.5 * maxVal) {
        double sum = 0;
        for (IntType i = 0; i <= k; ++i)
            sum += P(i);
        return sum;
    }
    double sum = 1.0;
    for (IntType i = k + 1; i <= maxVal; ++i)
        sum -= P(i);
    return sum;
}
 
template <typename IntType>
template <std::uniform_random_bit_generator G>
IntType NegativeHyperGeometricRand<IntType>::operator()(G &g) const
{
    const double expectedDraws = m * (static_cast<double>(N) + 1) / (static_cast<double>(M) + 1);
    if (expectedDraws < SMALL_URN)
        return variateByUrn(g);
    return BinomialDistribution<IntType>::Variate(N - M, successFraction(g), g);
}
 
template <typename IntType>
template <std::uniform_random_bit_generator G>
IntType NegativeHyperGeometricRand<IntType>::variateByUrn(G &g) const
{
    double p = static_cast<double>(M) / N;
    IntType successesNum = 0;
    IntType num = 0;
    while (successesNum < m) {
        ++num;
        if (BernoulliRand::Variate(p, g) && ++successesNum == num - N + M)
            return N - M;
        p = M - successesNum;
        p /= N - num;
    }
    return num - successesNum;
}
 
template <typename IntType>
long double NegativeHyperGeometricRand<IntType>::Mean() const
{
    long double mean = m;
    mean *= N - M;
    mean /= static_cast<long double>(M) + 1.0L;
    return mean;
}
 
template <typename IntType>
long double NegativeHyperGeometricRand<IntType>::Variance() const
{
    long double Mp1 = static_cast<long double>(M) + 1.0L;
    long double var = 1.0L - m / Mp1;
    var *= static_cast<long double>(N) + 1.0L;
    var *= N - M;
    var /= Mp1 * (Mp1 + 1);
    return m * var;
}
 
} // namespace randlib