continuous/gamma.hpp

Go to the documentation of this file.

#pragma once
 
#include <span>
 
#include <randlib/math/constants.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/continuous.hpp>
#include <randlib/core/exponential_family.hpp>
 
namespace randlib
{
 
template <typename RealType>
class GammaDistribution : virtual public ContinuousDistribution<RealType>
{
    RealType variate(detail::EngineRef g) const override { return (*this)(g); }
    void sample(detail::EngineRef g, std::span<RealType> outputData) const override;
 
protected:
    double alpha = 1; 
    double beta = 1; 
    double theta = 1; 
    double lgammaAlpha = 0; 
    double pdfCoef = 0; 
    double logAlpha = 0; 
    double logBeta = 0; 
    double digammaAlpha = -randlib::consts::egamma; 
 
private:
    struct genCoef_t
    {
        double t, b;
    } genCoef = {0, 0};
 
protected:
    GammaDistribution(double shape, double rate);
 
public:
    virtual ~GammaDistribution() {}
 
    [[nodiscard]] Support SupportType() const override { return Support::RightSemiFinite; }
    [[nodiscard]] RealType MinValue() const override { return 0; }
    [[nodiscard]] RealType MaxValue() const override { return INFINITY; }
 
protected:
    void SetParameters(double shape, double rate);
    void SetShape(double shape);
 
public:
    [[nodiscard]] inline double GetShape() const noexcept { return alpha; }
    [[nodiscard]] inline double GetScale() const noexcept { return theta; }
    [[nodiscard]] inline double GetRate() const noexcept { return beta; }
    [[nodiscard]] inline double GetLogGammaShape() const noexcept { return lgammaAlpha; }
    [[nodiscard]] inline double GetLogShape() const noexcept { return logAlpha; }
    [[nodiscard]] inline double GetDigammaShape() const noexcept { return digammaAlpha; }
    [[nodiscard]] inline double GetLogRate() const noexcept { return logBeta; }
 
    [[nodiscard]] double f(const RealType &x) const override;
    [[nodiscard]] double logf(const RealType &x) const override;
    [[nodiscard]] double F(const RealType &x) const override;
    double logF(const RealType &x) const;
    [[nodiscard]] double S(const RealType &x) const override;
    double logS(const RealType &x) const;
 
private:
    enum GENERATOR_ID
    {
        INTEGER_SHAPE, 
        ONE_AND_A_HALF_SHAPE, 
        SMALL_SHAPE, 
        FISHMAN, 
        MARSAGLIA_TSANG 
    };
 
    static GENERATOR_ID getIdOfUsedGenerator(double shape)
    {
        if (shape < 0.34)
            return SMALL_SHAPE;
        if (shape <= 3.0 && math::detail::areClose(shape, std::round(shape)))
            return INTEGER_SHAPE;
        if (math::detail::areClose(shape, 1.5))
            return ONE_AND_A_HALF_SHAPE;
        if (shape > 1.0 && shape < 1.2)
            return FISHMAN;
        return MARSAGLIA_TSANG;
    }
 
    template <std::uniform_random_bit_generator G>
    static RealType variateThroughExponentialSum(int shape, G &g);
    template <std::uniform_random_bit_generator G>
    static RealType variateForShapeOneAndAHalf(G &g);
    template <std::uniform_random_bit_generator G>
    RealType variateBest(G &g) const;
    template <std::uniform_random_bit_generator G>
    static RealType variateAhrensDieter(double shape, G &g);
    template <std::uniform_random_bit_generator G>
    static RealType variateFishman(double shape, G &g);
    template <std::uniform_random_bit_generator G>
    static RealType variateMarsagliaTsang(double shape, G &g);
 
public:
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] static RealType StandardVariate(double shape, G &g);
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] static RealType Variate(double shape, double rate, G &g);
 
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] RealType operator()(G &g) const;
 
    [[nodiscard]] long double Mean() const override;
    long double GeometricMean() const;
    [[nodiscard]] long double Variance() const override;
    long double GeometricVariance() const;
    [[nodiscard]] RealType Mode() const override;
    [[nodiscard]] RealType Median() const override;
    [[nodiscard]] long double Skewness() const override;
    [[nodiscard]] long double ExcessKurtosis() const override;
 
private:
    RealType initRootForSmallP(double r) const;
    RealType initRootForLargeP(double logQ) const;
    RealType initRootForLargeShape(double p) const;
    RealType initRootForLargeShape1m(double p) const;
    RealType quantileInitialGuess(double p) const;
    RealType quantileInitialGuess1m(double p) const;
 
    double df(RealType x) const;
    double dfDivf(RealType x) const;
    RealType quantileImpl(double p, RealType initValue) const override;
    RealType quantileImpl(double p) const override;
    RealType quantileImpl1m(double p, RealType initValue) const override;
    RealType quantileImpl1m(double p) const override;
 
    std::complex<double> CFImpl(double t) const override;
};
 
template <typename RealType>
class GammaRand;
 
template <typename RealType>
class FreeRateGammaDistribution : public GammaDistribution<RealType>
{
protected:
    FreeRateGammaDistribution(double shape, double rate) : GammaDistribution<RealType>(shape, rate) {}
 
public:
    void SetRate(double rate);
 
    void SetScale(double scale);
 
    RANDLIB_DEPRECATED_FIT_MEMBER void FitRate(std::span<const RealType> sample, bool unbiased = false);
 
    RANDLIB_DEPRECATED_FIT_MEMBER GammaRand<RealType> FitRateBayes(std::span<const RealType> sample,
                                                                   const GammaDistribution<RealType> &priorDistribution,
                                                                   bool MAP = false);
};
 
template <typename RealType>
class GammaRand final : public FreeRateGammaDistribution<RealType>, public ExponentialFamily<RealType, DoublePair>
{
public:
    GammaRand(double shape = 1, double rate = 1) : FreeRateGammaDistribution<RealType>(shape, rate) {}
 
    [[nodiscard]] static GammaRand FromShapeRate(double shape, double rate) { return GammaRand(shape, rate); }
 
    [[nodiscard]] static GammaRand FromShapeScale(double shape, double scale)
    {
        if (!(scale > 0.0))
            throw std::invalid_argument("Gamma distribution: scale should be positive");
        return GammaRand(shape, 1.0 / scale);
    }
    [[nodiscard]] std::string Name() const override;
 
    using GammaDistribution<RealType>::SetParameters;
    using GammaDistribution<RealType>::SetShape;
 
    RANDLIB_DEPRECATED_FIT_MEMBER void FitShape(std::span<const RealType> sample);
 
    RANDLIB_DEPRECATED_FIT_MEMBER void Fit(std::span<const RealType> sample);
 
    DoublePair SufficientStatistic(RealType x) const override;
    DoublePair SourceParameters() const override;
    DoublePair SourceToNatural(DoublePair sourceParameters) const override;
    double LogNormalizer(DoublePair parameters) const override;
    DoublePair LogNormalizerGradient(DoublePair parameters) const override;
    double CarrierMeasure(RealType) const override;
    double CrossEntropyAdjusted(DoublePair parameters) const override;
    double EntropyAdjusted() const override;
};
 
template <typename RealType>
class ChiSquaredRand final : public GammaDistribution<RealType>
{
public:
    explicit ChiSquaredRand(size_t degree = 1) : GammaDistribution<RealType>(0.5 * degree, 0.5) {}
    [[nodiscard]] std::string Name() const override;
    void SetDegree(size_t degree);
    [[nodiscard]] inline size_t GetDegree() const { return static_cast<int>(2 * this->alpha); }
};
 
template <typename RealType>
class ErlangRand final : public FreeRateGammaDistribution<RealType>
{
public:
    ErlangRand(int shape = 1, double rate = 1) : FreeRateGammaDistribution<RealType>(shape, rate) {}
 
    [[nodiscard]] static ErlangRand FromShapeRate(int shape, double rate) { return ErlangRand(shape, rate); }
 
    [[nodiscard]] static ErlangRand FromShapeScale(int shape, double scale)
    {
        if (!(scale > 0.0))
            throw std::invalid_argument("Erlang distribution: scale should be positive");
        return ErlangRand(shape, 1.0 / scale);
    }
    [[nodiscard]] std::string Name() const override;
    void SetParameters(size_t shape, double rate);
    void SetShape(size_t shape);
};
 
} // namespace randlib
 
#include <randlib/continuous/uniform.hpp>
#include <randlib/continuous/exponential.hpp>
#include <randlib/continuous/normal.hpp>
 
namespace randlib
{
 
template <typename RealType>
GammaDistribution<RealType>::GammaDistribution(double shape, double rate)
{
    SetParameters(shape, rate);
}
 
template <typename RealType>
void GammaDistribution<RealType>::SetParameters(double shape, double rate)
{
    if (!(shape > 0.0))
        throw std::invalid_argument("Gamma distribution: shape should be positive, but it's equal to " +
                                    std::to_string(shape));
    if (!(rate > 0.0))
        throw std::invalid_argument("Gamma distribution: rate should be positive, but it's equal to " +
                                    std::to_string(rate));
 
    this->alpha = shape > 0 ? shape : 1.0;
 
    this->beta = (rate > 0.0) ? rate : 1.0;
    theta = 1.0 / this->beta;
 
    lgammaAlpha = std::lgammal(this->alpha);
    logAlpha = std::log(this->alpha);
    this->logBeta = std::log(this->beta);
    pdfCoef = -lgammaAlpha + this->alpha * this->logBeta;
 
    if (getIdOfUsedGenerator(this->alpha) == SMALL_SHAPE) {
        genCoef.t = 0.5 * std::log1pl(-alpha);
        genCoef.t = 0.07 + 0.75 * std::exp(genCoef.t);
        genCoef.b = 1.0 + std::exp(-genCoef.t) * this->alpha / genCoef.t;
    }
}
 
template <typename RealType>
void GammaDistribution<RealType>::SetShape(double shape)
{
    SetParameters(shape, this->beta);
}
 
template <typename RealType>
double GammaDistribution<RealType>::f(const RealType &x) const
{
    if (x < 0.0)
        return 0.0;
    if (x == 0.0) {
        if (this->alpha > 1.0)
            return 0.0;
        return (this->alpha == 1.0) ? this->beta : INFINITY;
    }
    return std::exp(logf(x));
}
 
template <typename RealType>
double GammaDistribution<RealType>::logf(const RealType &x) const
{
    if (x < 0.0)
        return -INFINITY;
    if (x == 0.0) {
        if (this->alpha > 1.0)
            return -INFINITY;
        return (this->alpha == 1.0) ? this->logBeta : INFINITY;
    }
    double y = (this->alpha - 1.0) * std::log(x);
    y -= x * this->beta;
    y += pdfCoef;
    return y;
}
 
template <typename RealType>
double GammaDistribution<RealType>::F(const RealType &x) const
{
    return (x > 0.0) ? math::pgamma(this->alpha, x * this->beta, logAlpha, lgammaAlpha) : 0.0;
}
 
template <typename RealType>
double GammaDistribution<RealType>::logF(const RealType &x) const
{
    return (x > 0.0) ? math::lpgamma(this->alpha, x * this->beta, logAlpha, lgammaAlpha) : -INFINITY;
}
 
template <typename RealType>
double GammaDistribution<RealType>::S(const RealType &x) const
{
    return (x > 0.0) ? math::qgamma(this->alpha, x * this->beta, logAlpha, lgammaAlpha) : 1.0;
}
 
template <typename RealType>
double GammaDistribution<RealType>::logS(const RealType &x) const
{
    return (x > 0.0) ? math::lqgamma(this->alpha, x * this->beta, logAlpha, lgammaAlpha) : 0.0;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateThroughExponentialSum(int shape, G &g)
{
    double X = 0.0;
    for (int i = 0; i < shape; ++i)
        X += detail::standardExponential<RealType>(g);
    return X;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateForShapeOneAndAHalf(G &g)
{
    RealType W = detail::standardExponential<RealType>(g);
    RealType N = detail::standardNormal<RealType>(g);
    return W + 0.5 * N * N;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateBest(G &g) const
{
    double X = 0;
    size_t iter = 0;
    do {
        double V = genCoef.b * detail::uniform01<RealType>(g);
        double W = detail::uniform01<RealType>(g);
        if (V <= 1) {
            X = genCoef.t * std::pow(V, 1.0 / this->alpha);
            if (W <= (2.0 - X) / (2.0 + X) || W <= std::exp(-X))
                return X;
        } else {
            X = -std::log(genCoef.t * (genCoef.b - V) / this->alpha);
            double Y = X / genCoef.t;
            if (W * (this->alpha + Y - this->alpha * Y) <= 1 || W <= std::pow(Y, this->alpha - 1))
                return X;
        }
    } while (++iter <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Gamma distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateAhrensDieter(double shape, G &g)
{
    double X = 0;
    size_t iter = 0;
    double shapeInv = 1.0 / shape;
    double t = shapeInv + randlib::consts::inv_e;
    do {
        double U = detail::uniform01<RealType>(g);
        double p = shape * t * U;
        double W = detail::standardExponential<RealType>(g);
        if (p <= 1) {
            X = std::pow(p, shapeInv);
            if (X <= W)
                return X;
        } else {
            X = -std::log(t * (1 - U));
            if ((1 - shape) * std::log(X) <= W)
                return X;
        }
    } while (++iter <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Gamma distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateFishman(double shape, G &g)
{
    double W1, W2;
    double shapem1 = shape - 1;
    do {
        W1 = detail::standardExponential<RealType>(g);
        W2 = detail::standardExponential<RealType>(g);
    } while (W2 < shapem1 * (W1 - std::log(W1) - 1));
    return shape * W1;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::variateMarsagliaTsang(double shape, G &g)
{
    RealType d = shape - 1.0 / 3;
    RealType c = 3 * std::sqrt(d);
    size_t iter = 0;
    do {
        RealType N;
        do {
            N = detail::standardNormal<RealType>(g);
        } while (N <= -c);
        RealType v = 1 + N / c;
        v = v * v * v;
        N *= N;
        RealType U = detail::uniform01<RealType>(g);
        if (U < 1.0 - 0.331 * N * N || std::log(U) < 0.5 * N + d * (1.0 - v + std::log(v))) {
            return d * v;
        }
    } while (++iter <= ProbabilityDistribution<RealType>::MAX_ITER_REJECTION);
    throw std::runtime_error("Gamma distribution: sampling failed");
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::StandardVariate(double shape, G &g)
{
    if (!(shape > 0))
        throw std::invalid_argument("Gamma distribution: shape should be positive");
 
    GENERATOR_ID genId = getIdOfUsedGenerator(shape);
 
    switch (genId) {
    case INTEGER_SHAPE:
        return variateThroughExponentialSum(std::round(shape), g);
    case ONE_AND_A_HALF_SHAPE:
        return variateForShapeOneAndAHalf(g);
    case SMALL_SHAPE:
        return variateAhrensDieter(shape, g);
    case FISHMAN:
        return variateFishman(shape, g);
    case MARSAGLIA_TSANG:
        return variateMarsagliaTsang(shape, g);
    default:
        throw std::runtime_error("Gamma distribution: invalid generator id");
    }
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::Variate(double shape, double rate, G &g)
{
    if (!(shape > 0.0))
        throw std::invalid_argument("Gamma distribution: shape should be positive");
    if (!(rate > 0.0))
        throw std::invalid_argument("Gamma distribution: rate should be positive");
    return StandardVariate(shape, g) / rate;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType GammaDistribution<RealType>::operator()(G &g) const
{
    GENERATOR_ID genId = getIdOfUsedGenerator(this->alpha);
 
    switch (genId) {
    case INTEGER_SHAPE:
        return theta * variateThroughExponentialSum(this->alpha, g);
    case ONE_AND_A_HALF_SHAPE:
        return theta * variateForShapeOneAndAHalf(g);
    case SMALL_SHAPE:
        return theta * variateBest(g);
    case FISHMAN:
        return theta * variateFishman(this->alpha, g);
    case MARSAGLIA_TSANG:
        return theta * variateMarsagliaTsang(this->alpha, g);
    default:
        throw std::runtime_error("Gamma distribution: invalid generator id");
    }
}
 
template <typename RealType>
void GammaDistribution<RealType>::sample(detail::EngineRef g, std::span<RealType> outputData) const
{
    GENERATOR_ID genId = getIdOfUsedGenerator(this->alpha);
 
    switch (genId) {
    case INTEGER_SHAPE:
        for (RealType &var : outputData)
            var = theta * variateThroughExponentialSum(this->alpha, g);
        break;
    case ONE_AND_A_HALF_SHAPE:
        for (RealType &var : outputData)
            var = theta * variateForShapeOneAndAHalf(g);
        break;
    case SMALL_SHAPE:
        for (RealType &var : outputData)
            var = theta * variateBest(g);
        break;
    case FISHMAN:
        for (RealType &var : outputData)
            var = theta * variateFishman(this->alpha, g);
        break;
    case MARSAGLIA_TSANG:
        for (RealType &var : outputData)
            var = theta * variateMarsagliaTsang(this->alpha, g);
        break;
    default:
        return;
    }
}
 
template <typename RealType>
long double GammaDistribution<RealType>::Mean() const
{
    return this->alpha * theta;
}
 
template <typename RealType>
long double GammaDistribution<RealType>::GeometricMean() const
{
    return this->digammaAlpha - this->logBeta;
}
 
template <typename RealType>
long double GammaDistribution<RealType>::Variance() const
{
    return this->alpha * theta * theta;
}
 
template <typename RealType>
long double GammaDistribution<RealType>::GeometricVariance() const
{
    return math::trigamma(this->alpha);
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::Mode() const
{
    return (this->alpha <= 1) ? 0 : (this->alpha - 1) * theta;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::Median() const
{
    return (this->alpha == 1.0) ? theta * randlib::consts::ln2 : quantileImpl(0.5);
}
 
template <typename RealType>
long double GammaDistribution<RealType>::Skewness() const
{
    return 2.0l / std::sqrt(this->alpha);
}
 
template <typename RealType>
long double GammaDistribution<RealType>::ExcessKurtosis() const
{
    return 6.0l / this->alpha;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::initRootForSmallP(double r) const
{
    double root = 0;
    double c[5];
    c[4] = 1;
    double denominator = this->alpha + 1;
    c[3] = 1.0 / denominator;
    denominator *= denominator;
    denominator *= this->alpha + 2;
    c[2] = 0.5 * (3 * this->alpha + 5) / denominator;
    denominator *= 3 * (this->alpha + 1) * (this->alpha + 3);
    c[1] = 8 * this->alpha + 33;
    c[1] *= this->alpha;
    c[1] += 31;
    c[1] /= denominator;
    denominator *= 8 * (this->alpha + 1) * (this->alpha + 2) * (this->alpha + 4);
    c[0] = 125 * this->alpha + 1179;
    c[0] *= this->alpha;
    c[0] += 3971;
    c[0] *= this->alpha;
    c[0] += 5661;
    c[0] *= this->alpha;
    c[0] += 2888;
    c[0] /= denominator;
    for (int i = 0; i != 5; ++i) {
        root += c[i];
        root *= r;
    }
    return root;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::initRootForLargeP(double logQ) const
{
    double x = (logQ + lgammaAlpha) / this->alpha;
    x = -std::exp(x) / this->alpha;
    return -alpha * math::Wm1Lambert(x);
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::initRootForLargeShape(double p) const
{
    if (p == 0.5)
        return this->alpha;
    double x = math::erfcinv(2 * p);
    double lambda = x * x / this->alpha + 1;
    lambda = -std::exp(-lambda);
    lambda = (x < 0) ? -math::Wm1Lambert(lambda) : -math::W0Lambert(lambda);
    return lambda * this->alpha;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::initRootForLargeShape1m(double p) const
{
    if (p == 0.5)
        return this->alpha;
    double x = -math::erfcinv(2 * p);
    double lambda = x * x / this->alpha + 1;
    lambda = -std::exp(-lambda);
    lambda = (x < 0) ? -math::Wm1Lambert(lambda) : -math::W0Lambert(lambda);
    return lambda * this->alpha;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileInitialGuess(double p) const
{
    double guess = 0;
    if (this->alpha < 10) {
        double r = std::log(p * this->alpha) + lgammaAlpha;
        r = std::exp(r / this->alpha);
        if (r < 0.2 * (this->alpha + 1)) {
            guess = initRootForSmallP(r);
        } else {
            double logQ = std::log1pl(-p);
            double maxBoundary1 = -0.5 * this->alpha - logAlpha - lgammaAlpha;
            double maxBoundary2 = this->alpha * (logAlpha - 1) - lgammaAlpha;
            if (logQ < std::min(maxBoundary1, maxBoundary2))
                guess = initRootForLargeP(logQ);
            else if (this->alpha < 1)
                guess = r;
            else
                guess = initRootForLargeShape(p);
        }
    } else
        guess = initRootForLargeShape(p);
    return guess / this->beta;
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileInitialGuess1m(double p) const
{
    if (this->alpha < 10) {
        double logQ = std::log(p);
        double maxBoundary1 = -0.5 * this->alpha - logAlpha - lgammaAlpha;
        double maxBoundary2 = this->alpha * (logAlpha - 1) - lgammaAlpha;
        if (logQ < std::min(maxBoundary1, maxBoundary2))
            return initRootForLargeP(logQ) / this->beta;
    } else {
        return initRootForLargeShape1m(p) / this->beta;
    }
    return quantileInitialGuess(1.0 - p);
}
 
template <typename RealType>
double GammaDistribution<RealType>::df(RealType x) const
{
    double z = (this->alpha - 1) - this->beta * x;
    double y = (this->alpha - 2) * std::log(x);
    y -= this->beta * x;
    y += pdfCoef;
    return z * std::exp(y);
}
 
template <typename RealType>
double GammaDistribution<RealType>::dfDivf(RealType x) const
{
    return x / (this->alpha - 1 - this->beta * x);
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileImpl(double p, RealType initValue) const
{
    return ContinuousDistribution<RealType>::quantileImpl(p, initValue);
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileImpl(double p) const
{
    return (this->alpha == 1.0) ? -theta * std::log1pl(-p) : quantileImpl(p, quantileInitialGuess(p));
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileImpl1m(double p, RealType initValue) const
{
    return ContinuousDistribution<RealType>::quantileImpl1m(p, initValue);
}
 
template <typename RealType>
RealType GammaDistribution<RealType>::quantileImpl1m(double p) const
{
    return (this->alpha == 1.0) ? -theta * std::log(p) : quantileImpl1m(p, quantileInitialGuess1m(p));
}
 
template <typename RealType>
std::complex<double> GammaDistribution<RealType>::CFImpl(double t) const
{
    return std::pow(std::complex<double>(1.0, -theta * t), -alpha);
}
// FREE-SCALE-GAMMA-DISTRIBUTION
 
template <typename RealType>
void FreeRateGammaDistribution<RealType>::SetRate(double rate)
{
    this->SetParameters(this->alpha, rate);
}
 
template <typename RealType>
void FreeRateGammaDistribution<RealType>::SetScale(double scale)
{
    SetRate(1.0 / scale);
}
 
template <typename RealType>
void FreeRateGammaDistribution<RealType>::FitRate(std::span<const RealType> sample, bool unbiased)
{
    if (!this->allElementsArePositive(sample))
        throw this->fitError(FitErrorCode::WrongSample, this->POSITIVITY_VIOLATION);
    double mean = this->GetSampleMean(sample);
    double coef = this->alpha - (unbiased ? 1.0 / sample.size() : 0.0);
    this->SetParameters(this->alpha, coef / mean);
}
 
template <typename RealType>
GammaRand<RealType>
FreeRateGammaDistribution<RealType>::FitRateBayes(std::span<const RealType> sample,
                                                  const GammaDistribution<RealType> &priorDistribution, bool MAP)
{
    if (!this->allElementsArePositive(sample))
        throw this->fitError(FitErrorCode::WrongSample, this->POSITIVITY_VIOLATION);
    double kappa = priorDistribution.GetShape();
    double gamma = priorDistribution.GetRate();
    double newShape = this->alpha * sample.size() + kappa;
    double newRate = this->GetSampleSum(sample) + gamma;
    GammaRand<RealType> posteriorDistribution(newShape, newRate);
    this->SetParameters(this->alpha, MAP ? posteriorDistribution.Mode() : posteriorDistribution.Mean());
    return posteriorDistribution;
}
// GAMMARAND
 
template <typename RealType>
std::string GammaRand<RealType>::Name() const
{
    return "Gamma(" + this->toStringWithPrecision(this->GetShape()) + ", " +
           this->toStringWithPrecision(this->GetRate()) + ")";
}
 
template <typename RealType>
void GammaRand<RealType>::FitShape(std::span<const RealType> sample)
{
    if (!this->allElementsArePositive(sample))
        throw this->fitError(FitErrorCode::WrongSample, this->POSITIVITY_VIOLATION);
 
    double shape = this->GetSampleMean(sample) * this->beta;
    double s = this->GetSampleLogMean(sample) + this->logBeta;
    if (!math::findRootNewtonFirstOrder<double>(
            [s](double x) {
                double first = math::digamma(x) - s;
                double second = math::trigamma(x);
                return DoublePair(first, second);
            },
            shape))
        throw this->fitError(FitErrorCode::Undefined, "Error in root-finding procedure");
    SetParameters(shape, this->beta);
}
 
template <typename RealType>
void GammaRand<RealType>::Fit(std::span<const RealType> sample)
{
    if (!this->allElementsArePositive(sample))
        throw this->fitError(FitErrorCode::WrongSample, this->POSITIVITY_VIOLATION);
 
    double average = this->GetSampleMean(sample);
    double s = std::log(average) - this->GetSampleLogMean(sample);
    double sm3 = s - 3.0, sp12 = 12.0 * s;
    double shape = sm3 * sm3 + 2 * sp12;
    shape = std::sqrt(shape);
    shape -= sm3;
    shape /= sp12;
 
    if (!math::findRootNewtonFirstOrder<double>(
            [s](double x) {
                double first = math::digammamLog(x) + s;
                double second = math::trigamma(x) - 1.0 / x;
                return DoublePair(first, second);
            },
            shape))
        throw this->fitError(FitErrorCode::Undefined, "Error in root-finding procedure");
 
    SetParameters(shape, shape / average);
}
 
template <typename RealType>
DoublePair GammaRand<RealType>::SufficientStatistic(RealType x) const
{
    return {std::log(x), x};
}
 
template <typename RealType>
DoublePair GammaRand<RealType>::SourceParameters() const
{
    return {this->alpha, this->beta};
}
 
template <typename RealType>
DoublePair GammaRand<RealType>::SourceToNatural(DoublePair sourceParameters) const
{
    double shape = sourceParameters.first;
    double rate = sourceParameters.second;
    return {shape - 1, -rate};
}
 
template <typename RealType>
double GammaRand<RealType>::LogNormalizer(DoublePair parameters) const
{
    double shape = parameters.first + 1, rate = -parameters.second;
    double F = std::lgamma(shape);
    F -= shape * std::log(rate);
    return F;
}
 
template <typename RealType>
DoublePair GammaRand<RealType>::LogNormalizerGradient(DoublePair parameters) const
{
    double shape = parameters.first + 1, rate = -parameters.second;
    double gradF1 = math::digamma(shape) - std::log(rate);
    double gradF2 = shape / rate;
    return {gradF1, gradF2};
}
 
template <typename RealType>
double GammaRand<RealType>::CarrierMeasure(RealType) const
{
    return 0;
}
 
template <typename RealType>
double GammaRand<RealType>::CrossEntropyAdjusted(DoublePair parameters) const
{
    double shapeq = parameters.first + 1, rateq = -parameters.second;
    double H = std::lgamma(shapeq);
    H -= shapeq * std::log(rateq);
    H -= (shapeq - 1) * (this->digammaAlpha - this->logBeta);
    H += rateq * this->alpha / this->beta;
    return H;
}
 
template <typename RealType>
double GammaRand<RealType>::EntropyAdjusted() const
{
    double H = this->lgammaAlpha;
    H -= this->logBeta;
    H -= (this->alpha - 1) * this->digammaAlpha;
    H += this->alpha;
    return H;
}
// CHI-SQUAREDRAND
 
template <typename RealType>
std::string ChiSquaredRand<RealType>::Name() const
{
    return "Chi-squared(" + this->toStringWithPrecision(GetDegree()) + ")";
}
 
template <typename RealType>
void ChiSquaredRand<RealType>::SetDegree(size_t degree)
{
    GammaDistribution<RealType>::SetParameters(0.5 * degree, 0.5);
}
// ERLANGRAND
 
template <typename RealType>
std::string ErlangRand<RealType>::Name() const
{
    return "Erlang(" + this->toStringWithPrecision(this->GetShape()) + ", " +
           this->toStringWithPrecision(this->GetRate()) + ")";
}
 
template <typename RealType>
void ErlangRand<RealType>::SetParameters(size_t shape, double rate)
{
    GammaDistribution<RealType>::SetParameters(shape, rate);
}
 
template <typename RealType>
void ErlangRand<RealType>::SetShape(size_t shape)
{
    SetParameters(shape, this->beta);
}
 
} // namespace randlib