distribution.hpp

Go to the documentation of this file.

#pragma once
 
#include <cstdint>
#include <iomanip>
#include <optional>
#include <random>
#include <span>
#include <sstream>
#include <string>
 
#include <randlib/math/constants.hpp>
#include <randlib/core/fwd.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/variates.hpp>
#include <randlib/engine/seeding.hpp>
#include <randlib/engine/uniform.hpp>
 
namespace randlib
{
 
template <typename T>
class ProbabilityDistribution
{
    mutable std::optional<DefaultEngine> ownEngine{};
 
protected:
    static constexpr size_t MAX_ITER_REJECTION = detail::MAX_ITER_REJECTION;
 
    std::string toStringWithPrecision(const double a_value, const int n = 6) const;
 
    ProbabilityDistribution() = default;
 
    DefaultEngine &engine() const
    {
        if (!ownEngine)
            ownEngine.emplace(random_seed());
        return *ownEngine;
    }
 
    virtual T variate(detail::EngineRef g) const = 0;
 
    virtual void sample(detail::EngineRef g, std::span<T> outputData) const
    {
        for (T &var : outputData)
            var = variate(g);
    }
 
public:
    using value_type = T;
 
    virtual ~ProbabilityDistribution() = default;
 
    [[nodiscard]] virtual std::string Name() const = 0;
 
    [[nodiscard]] virtual T MinValue() const = 0;
 
    [[nodiscard]] virtual T MaxValue() const = 0;
 
    [[nodiscard]] virtual double F(const T &x) const = 0;
 
    void CumulativeDistributionFunction(std::span<const T> x, std::span<double> y) const;
 
    [[nodiscard]] virtual double S(const T &x) const;
 
    void SurvivalFunction(std::span<const T> x, std::span<double> y) const;
 
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] T operator()(G &g) const
    {
        return variate(g);
    }
 
    [[nodiscard]] T Variate() const { return variate(engine()); }
 
    void Sample(std::span<T> outputData) const { sample(engine(), outputData); }
 
    template <std::uniform_random_bit_generator G>
    void Sample(G &g, std::span<T> outputData) const
    {
        sample(g, outputData);
    }
 
    void Reseed(std::uint64_t seed) { ownEngine.emplace(seed); }
 
protected:
    static constexpr char POSITIVITY_VIOLATION[] = "All elements should be positive";
    static constexpr char NON_NEGATIVITY_VIOLATION[] = "All elements should be non-negative";
    static constexpr char UPPER_LIMIT_VIOLATION[] = "No element should be greater than ";
    static constexpr char LOWER_LIMIT_VIOLATION[] = "No element should be smaller than ";
 
    fit_error fitError(FitErrorCode code, const std::string &explanation) const;
};
 
} // namespace randlib
 
#include <randlib/math/quadrature.hpp>
#include <randlib/math/roots.hpp>
#include <randlib/math/optimize.hpp>
 
namespace randlib
{
 
template <typename T>
std::string ProbabilityDistribution<T>::toStringWithPrecision(const double a_value, const int n) const
{
    std::ostringstream out;
    out << std::setprecision(n) << a_value;
    return out.str();
}
 
template <typename T>
void ProbabilityDistribution<T>::CumulativeDistributionFunction(std::span<const T> 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->F(x[i]);
}
 
template <typename T>
double ProbabilityDistribution<T>::S(const T &x) const
{
    return 1.0 - this->F(x);
}
 
template <typename T>
void ProbabilityDistribution<T>::SurvivalFunction(std::span<const T> 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->S(x[i]);
}
 
template <typename T>
constexpr char ProbabilityDistribution<T>::POSITIVITY_VIOLATION[];
template <typename T>
constexpr char ProbabilityDistribution<T>::NON_NEGATIVITY_VIOLATION[];
template <typename T>
constexpr char ProbabilityDistribution<T>::UPPER_LIMIT_VIOLATION[];
template <typename T>
constexpr char ProbabilityDistribution<T>::LOWER_LIMIT_VIOLATION[];
 
template <typename T>
fit_error ProbabilityDistribution<T>::fitError(FitErrorCode code, const std::string &explanation) const
{
    std::string reason;
    switch (code) {
    case FitErrorCode::WrongSample:
        reason = "Sample couldn't be generated by this distribution. ";
        break;
    case FitErrorCode::NotApplicable:
        reason = "Method cannot be applied here. ";
        break;
    case FitErrorCode::WrongReturn:
        reason = "Method returns invalid parameters. ";
        break;
    case FitErrorCode::TooFewElements:
        reason = "Sample is too small. ";
        break;
    case FitErrorCode::WrongLevel:
        reason = "Significance level should be positive and smaller than 1. ";
        break;
    case FitErrorCode::Undefined:
        reason = "Unknown type of error. ";
        break;
    }
    return fit_error(code, this->Name() + ": " + reason + explanation);
}
 
 
} // namespace randlib