boost/random/generalized_inverse_gaussian_distribution.hpp

0001 /* boost random/generalized_inverse_gaussian_distribution.hpp header file
0002  *
0003  * Copyright Young Geun Kim 2025
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * $Id$
0011  */
0012
0013 #ifndef BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP
0014 #define BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP
0015
0016 #include <boost/config/no_tr1/cmath.hpp>
0017 #include <istream>
0018 #include <iosfwd>
0019 #include <limits>
0020 #include <boost/assert.hpp>
0021 #include <boost/limits.hpp>
0022 #include <boost/random/detail/config.hpp>
0023 #include <boost/random/detail/operators.hpp>
0024 #include <boost/random/uniform_01.hpp>
0025
0026 namespace boost {
0027 namespace random {
0028
0029 /**
0030  * The generalized inverse gaussian distribution is a real-valued distribution with
0031  * three parameters p, a, and b. It produced values > 0.
0032  *
0033  * It has
0034  * \f$\displaystyle p(x) = \frac{(a / b)^{p / 2}}{2 K_{p}(\sqrt{a b})} x^{p - 1} e^{-(a x + b / 2) / 2}\f$.
0035  * where \f$\displaystyle K_{p}\f$ is a modified Bessel function of the second kind.
0036  *
0037  * The algorithm used is from
0038  *
0039  * @blockquote
0040  * "Random variate generation for the generalized inverse Gaussian distribution",
0041  * Luc Devroye,
0042  * Statistics and Computing,
0043  * Volume 24, 2014, Pages 236 - 246
0044  * @endblockquote
0045  */
0046 template<class RealType = double>
0047 class generalized_inverse_gaussian_distribution
0048 {
0049 public:
0050     typedef RealType result_type;
0051     typedef RealType input_type;
0052
0053     class param_type {
0054     public:
0055         typedef generalized_inverse_gaussian_distribution distribution_type;
0056
0057         /**
0058          * Constructs a @c param_type object from the "p", "a", and "b"
0059          * parameters.
0060          *
0061          * Requires:
0062              * a > 0 && b >= 0 if p > 0,
0063              * a > 0 && b > 0 if p == 0,
0064              * a >= 0 && b > 0 if p < 0
0065          */
0066         explicit param_type(RealType p_arg = RealType(1.0),
0067                             RealType a_arg = RealType(1.0),
0068                             RealType b_arg = RealType(1.0))
0069         : _p(p_arg), _a(a_arg), _b(b_arg)
0070         {
0071             BOOST_ASSERT(
0072                 (p_arg > RealType(0) && a_arg > RealType(0) && b_arg >= RealType(0)) ||
0073                 (p_arg == RealType(0) && a_arg > RealType(0) && b_arg > RealType(0)) ||
0074                 (p_arg < RealType(0) && a_arg >= RealType(0) && b_arg > RealType(0))
0075             );
0076         }
0077
0078         /** Returns the "p" parameter of the distribution. */
0079         RealType p() const { return _p; }
0080         /** Returns the "a" parameter of the distribution. */
0081         RealType a() const { return _a; }
0082         /** Returns the "b" parameter of the distribution. */
0083         RealType b() const { return _b; }
0084
0085         /** Writes a @c param_type to a @c std::ostream. */
0086         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0087         {
0088             os << parm._p << ' ' << parm._a << ' ' << parm._b;
0089             return os;
0090         }
0091
0092         /** Reads a @c param_type from a @c std::istream. */
0093         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0094         {
0095             is >> parm._p >> std::ws >> parm._a >> std::ws >> parm._b;
0096             return is;
0097         }
0098
0099         /** Returns true if the two sets of parameters are the same. */
0100         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0101         { return lhs._p == rhs._p && lhs._a == rhs._a && lhs._b == rhs._b; }
0102
0103         /** Returns true if the two sets of parameters are different. */
0104         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0105
0106     private:
0107         RealType _p;
0108         RealType _a;
0109         RealType _b;
0110     };
0111
0112 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
0113     BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
0114 #endif
0115
0116     /**
0117      * Constructs an @c generalized_inverse_gaussian_distribution from its "p", "a", and "b" parameters.
0118      *
0119      * Requires:
0120      * a > 0 && b >= 0 if p > 0,
0121      * a > 0 && b > 0 if p == 0,
0122      * a >= 0 && b > 0 if p < 0
0123      */
0124     explicit generalized_inverse_gaussian_distribution(RealType p_arg = RealType(1.0),
0125                                                        RealType a_arg = RealType(1.0),
0126                                                        RealType b_arg = RealType(1.0))
0127     : _p(p_arg), _a(a_arg), _b(b_arg)
0128     {
0129         BOOST_ASSERT(
0130             (p_arg > RealType(0) && a_arg > RealType(0) && b_arg >= RealType(0)) ||
0131             (p_arg == RealType(0) && a_arg > RealType(0) && b_arg > RealType(0)) ||
0132             (p_arg < RealType(0) && a_arg >= RealType(0) && b_arg > RealType(0))
0133         );
0134         init();
0135     }
0136     /** Constructs an @c generalized_inverse_gaussian_distribution from its parameters. */
0137     explicit generalized_inverse_gaussian_distribution(const param_type& parm)
0138     : _p(parm.p()), _a(parm.a()), _b(parm.b())
0139     {
0140         init();
0141     }
0142
0143     /**
0144      * Returns a random variate distributed according to the
0145      * generalized inverse gaussian distribution.
0146      */
0147     template<class URNG>
0148     RealType operator()(URNG& urng) const
0149     {
0150 #ifndef BOOST_NO_STDC_NAMESPACE
0151         using std::sqrt;
0152         using std::log;
0153         using std::min;
0154         using std::exp;
0155         using std::cosh;
0156 #endif
0157         RealType t = result_type(1);
0158         RealType s = result_type(1);
0159         RealType log_concave = -psi(result_type(1));
0160         if (log_concave >= result_type(.5) && log_concave <= result_type(2)) {
0161             t = result_type(1);
0162         } else if (log_concave > result_type(2)) {
0163             t = sqrt(result_type(2) / (_alpha + _abs_p));
0164         } else if (log_concave < result_type(.5)) {
0165             t = log(result_type(4) / (_alpha + result_type(2) * _abs_p));
0166         }
0167         log_concave = -psi(result_type(-1));
0168         if (log_concave >= result_type(.5) && log_concave <= result_type(2)) {
0169             s = result_type(1);
0170         } else if (log_concave > result_type(2)) {
0171             s = sqrt(result_type(4) / (_alpha * cosh(1) + _abs_p));
0172         } else if (log_concave < result_type(.5)) {
0173             s = min(result_type(1) / _abs_p, log(result_type(1) + result_type(1) / _alpha + sqrt(result_type(1) / (_alpha * _alpha) + result_type(2) / _alpha)));
0174         }
0175         RealType eta = -psi(t);
0176         RealType zeta = -psi_deriv(t);
0177         RealType theta = -psi(-s);
0178         RealType xi = psi_deriv(-s);
0179         RealType p = result_type(1) / xi;
0180         RealType r = result_type(1) / zeta;
0181         RealType t_deriv = t - r * eta;
0182         RealType s_deriv = s - p * theta;
0183         RealType q = t_deriv + s_deriv;
0184         RealType u = result_type(0);
0185         RealType v = result_type(0);
0186         RealType w = result_type(0);
0187         RealType cand = result_type(0);
0188         do
0189         {
0190             u = uniform_01<RealType>()(urng);
0191             v = uniform_01<RealType>()(urng);
0192             w = uniform_01<RealType>()(urng);
0193             if (u < q / (p + q + r)) {
0194                 cand = -s_deriv + q * v;
0195             } else if (u < (q + r) / (p + q + r)) {
0196                 cand = t_deriv - r * log(v);
0197             } else {
0198                 cand = -s_deriv + p * log(v);
0199             }
0200         } while (w * chi(cand, s, t, s_deriv, t_deriv, eta, zeta, theta, xi) > exp(psi(cand)));
0201         cand = (_abs_p / _omega + sqrt(result_type(1) + _abs_p * _abs_p / (_omega * _omega))) * exp(cand);
0202         return _p > 0 ? cand * sqrt(_b / _a) : sqrt(_b / _a) / cand;
0203     }
0204
0205     /**
0206      * Returns a random variate distributed accordint to the beta
0207      * distribution with parameters specified by @c param.
0208      */
0209     template<class URNG>
0210     result_type operator()(URNG& urng, const param_type& parm) const
0211     {
0212         return generalized_inverse_gaussian_distribution(parm)(urng);
0213     }
0214
0215     /** Returns the "p" parameter of the distribution. */
0216     RealType p() const { return _p; }
0217     /** Returns the "a" parameter of the distribution. */
0218     RealType a() const { return _a; }
0219     /** Returns the "b" parameter of the distribution. */
0220     RealType b() const { return _b; }
0221
0222     /** Returns the smallest value that the distribution can produce. */
0223     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0224     { return RealType(0.0); }
0225     /** Returns the largest value that the distribution can produce. */
0226     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0227     { return (std::numeric_limits<RealType>::infinity)(); }
0228
0229     /** Returns the parameters of the distribution. */
0230     param_type param() const { return param_type(_p, _a, _b); }
0231     /** Sets the parameters of the distribution. */
0232     void param(const param_type& parm)
0233     {
0234         _p = parm.p();
0235         _a = parm.a();
0236         _b = parm.b();
0237         init();
0238     }
0239
0240     /**
0241      * Effects: Subsequent uses of the distribution do not depend
0242      * on values produced by any engine prior to invoking reset.
0243      */
0244     void reset() { }
0245
0246     /** Writes an @c generalized_inverse_gaussian_distribution to a @c std::ostream. */
0247     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, generalized_inverse_gaussian_distribution, wd)
0248     {
0249         os << wd.param();
0250         return os;
0251     }
0252
0253     /** Reads an @c generalized_inverse_gaussian_distribution from a @c std::istream. */
0254     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, generalized_inverse_gaussian_distribution, wd)
0255     {
0256         param_type parm;
0257         if(is >> parm) {
0258             wd.param(parm);
0259         }
0260         return is;
0261     }
0262
0263     /**
0264      * Returns true if the two instances of @c generalized_inverse_gaussian_distribution will
0265      * return identical sequences of values given equal generators.
0266      */
0267     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(generalized_inverse_gaussian_distribution, lhs, rhs)
0268     { return lhs._p == rhs._p && lhs._a == rhs._a && lhs._b == rhs._b; }
0269
0270     /**
0271      * Returns true if the two instances of @c generalized_inverse_gaussian_distribution will
0272      * return different sequences of values given equal generators.
0273      */
0274     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(generalized_inverse_gaussian_distribution)
0275
0276 private:
0277     RealType _p;
0278     RealType _a;
0279     RealType _b;
0280     // some data precomputed from the parameters
0281     RealType _abs_p;
0282     RealType _omega;
0283     RealType _alpha;
0284
0285     /// \cond hide_private_members
0286     void init()
0287     {
0288 #ifndef BOOST_NO_STDC_NAMESPACE
0289         using std::abs;
0290         using std::sqrt;
0291 #endif
0292         _abs_p = abs(_p);
0293         _omega = sqrt(_a * _b); // two-parameter representation (p, omega)
0294         _alpha = sqrt(_omega * _omega + _abs_p * _abs_p) - _abs_p;
0295     }
0296
0297     result_type psi(const RealType& x) const
0298     {
0299 #ifndef BOOST_NO_STDC_NAMESPACE
0300         using std::cosh;
0301         using std::exp;
0302 #endif
0303         return -_alpha * (cosh(x) - result_type(1)) - _abs_p * (exp(x) - x - result_type(1));
0304     }
0305
0306     result_type psi_deriv(const RealType& x) const
0307     {
0308 #ifndef BOOST_NO_STDC_NAMESPACE
0309         using std::sinh;
0310         using std::exp;
0311 #endif
0312         return -_alpha * sinh(x) - _abs_p * (exp(x) - result_type(1));
0313     }
0314
0315     static result_type chi(const RealType& x,
0316                            const RealType& s, const RealType& t,
0317                            const RealType& s_deriv, const RealType& t_deriv,
0318                            const RealType& eta, const RealType& zeta, const RealType& theta, const RealType& xi)
0319     {
0320 #ifndef BOOST_NO_STDC_NAMESPACE
0321         using std::exp;
0322 #endif
0323         if (x >= -s_deriv && x <= t_deriv) {
0324             return result_type(1);
0325         } else if (x > t_deriv) {
0326             return exp(-eta - zeta * (x - t));
0327         }
0328         return exp(-theta + xi * (x + s));
0329     }
0330     /// \endcond
0331 };
0332
0333 } // namespace random
0334
0335 using random::generalized_inverse_gaussian_distribution;
0336
0337 } // namespace boost
0338
0339 #endif // BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP