mirror of
https://github.com/boostorg/random.git
synced 2026-07-21 13:33:46 +00:00
Feature: Generalized inverse Gaussian random generator (#126)
* use Devroye instead of Hormann for GIG * add gig header file * add gig test * move using declarations before member declaration * move using into each function * fix tab size and add init() * add test to jamfile * remove duplicate iostream operators * fix formatting again * add to docs * remove static because they use non-static member * add const in member functions * remove duplicate test * fix typo * need dist3 in three parameter test * specify test1 params * missing using std::cosh * sinh and cosh under c++11 * set arg3 different in test * set test arg value different from the default * add test3 to cover every line in gig * remove uninitialized variables compile warning in gig operator * did not add test3 mistakenly
This commit is contained in:
Young Geun Kim
@@ -72,6 +72,7 @@ doxygen_files =
|
||||
splitmix64
|
||||
xoshiro
|
||||
inverse_gaussian_distribution
|
||||
generalized_inverse_gaussian_distribution
|
||||
;
|
||||
|
||||
path-constant here : . ;
|
||||
|
||||
@@ -108,4 +108,5 @@ statistically to it are not acceptable.
|
||||
[choosing a random point on Earth (assumed to be a
|
||||
sphere) where to spend the next vacations]]
|
||||
[[__inverse_gaussian_distribution] [Inverse Gaussian Distribution] [The distribution is used to model non-negative, positively skewed data and has a wide variety of applications in business, survival analysis, finance, medicine, and even in labor dispute resolution.]]
|
||||
[[__generalized_inverse_gaussian_distribution][Generalized inverse Gaussian Distribution][-]]
|
||||
]
|
||||
|
||||
@@ -109,6 +109,7 @@
|
||||
[def __piecewise_constant_distribution [classref boost::random::piecewise_constant_distribution piecewise_constant_distribution]]
|
||||
[def __piecewise_linear_distribution [classref boost::random::piecewise_linear_distribution piecewise_linear_distribution]]
|
||||
[def __inverse_gaussian_distribution [classref boost::random::inverse_gaussian_distribution inverse_gaussian_distribution]]
|
||||
[def __generalized_inverse_gaussian_distribution [classref boost::random::generalized_inverse_gaussian_distribution generalized_inverse_gaussian_distribution]]
|
||||
|
||||
[include performance_data.qbk]
|
||||
|
||||
|
||||
@@ -89,6 +89,7 @@
|
||||
#include <boost/random/uniform_smallint.hpp>
|
||||
#include <boost/random/weibull_distribution.hpp>
|
||||
#include <boost/random/inverse_gaussian_distribution.hpp>
|
||||
#include <boost/random/generalized_inverse_gaussian_distribution.hpp>
|
||||
|
||||
#include <boost/random/generate_canonical.hpp>
|
||||
|
||||
|
||||
@@ -0,0 +1,339 @@
|
||||
/* boost random/generalized_inverse_gaussian_distribution.hpp header file
|
||||
*
|
||||
* Copyright Young Geun Kim 2025
|
||||
* Distributed under the Boost Software License, Version 1.0. (See
|
||||
* accompanying file LICENSE_1_0.txt or copy at
|
||||
* http://www.boost.org/LICENSE_1_0.txt)
|
||||
*
|
||||
* See http://www.boost.org for most recent version including documentation.
|
||||
*
|
||||
* $Id$
|
||||
*/
|
||||
|
||||
#ifndef BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP
|
||||
#define BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP
|
||||
|
||||
#include <boost/config/no_tr1/cmath.hpp>
|
||||
#include <istream>
|
||||
#include <iosfwd>
|
||||
#include <limits>
|
||||
#include <boost/assert.hpp>
|
||||
#include <boost/limits.hpp>
|
||||
#include <boost/random/detail/config.hpp>
|
||||
#include <boost/random/detail/operators.hpp>
|
||||
#include <boost/random/uniform_01.hpp>
|
||||
|
||||
namespace boost {
|
||||
namespace random {
|
||||
|
||||
/**
|
||||
* The generalized inverse gaussian distribution is a real-valued distribution with
|
||||
* three parameters p, a, and b. It produced values > 0.
|
||||
*
|
||||
* It has
|
||||
* \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$.
|
||||
* where \f$\displaystyle K_{p}\f$ is a modified Bessel function of the second kind.
|
||||
*
|
||||
* The algorithm used is from
|
||||
*
|
||||
* @blockquote
|
||||
* "Random variate generation for the generalized inverse Gaussian distribution",
|
||||
* Luc Devroye,
|
||||
* Statistics and Computing,
|
||||
* Volume 24, 2014, Pages 236 - 246
|
||||
* @endblockquote
|
||||
*/
|
||||
template<class RealType = double>
|
||||
class generalized_inverse_gaussian_distribution
|
||||
{
|
||||
public:
|
||||
typedef RealType result_type;
|
||||
typedef RealType input_type;
|
||||
|
||||
class param_type {
|
||||
public:
|
||||
typedef generalized_inverse_gaussian_distribution distribution_type;
|
||||
|
||||
/**
|
||||
* Constructs a @c param_type object from the "p", "a", and "b"
|
||||
* parameters.
|
||||
*
|
||||
* Requires:
|
||||
* a > 0 && b >= 0 if p > 0,
|
||||
* a > 0 && b > 0 if p == 0,
|
||||
* a >= 0 && b > 0 if p < 0
|
||||
*/
|
||||
explicit param_type(RealType p_arg = RealType(1.0),
|
||||
RealType a_arg = RealType(1.0),
|
||||
RealType b_arg = RealType(1.0))
|
||||
: _p(p_arg), _a(a_arg), _b(b_arg)
|
||||
{
|
||||
BOOST_ASSERT(
|
||||
(p_arg > RealType(0) && a_arg > RealType(0) && b_arg >= RealType(0)) ||
|
||||
(p_arg == RealType(0) && a_arg > RealType(0) && b_arg > RealType(0)) ||
|
||||
(p_arg < RealType(0) && a_arg >= RealType(0) && b_arg > RealType(0))
|
||||
);
|
||||
}
|
||||
|
||||
/** Returns the "p" parameter of the distribution. */
|
||||
RealType p() const { return _p; }
|
||||
/** Returns the "a" parameter of the distribution. */
|
||||
RealType a() const { return _a; }
|
||||
/** Returns the "b" parameter of the distribution. */
|
||||
RealType b() const { return _b; }
|
||||
|
||||
/** Writes a @c param_type to a @c std::ostream. */
|
||||
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
|
||||
{
|
||||
os << parm._p << ' ' << parm._a << ' ' << parm._b;
|
||||
return os;
|
||||
}
|
||||
|
||||
/** Reads a @c param_type from a @c std::istream. */
|
||||
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
|
||||
{
|
||||
is >> parm._p >> std::ws >> parm._a >> std::ws >> parm._b;
|
||||
return is;
|
||||
}
|
||||
|
||||
/** Returns true if the two sets of parameters are the same. */
|
||||
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
|
||||
{ return lhs._p == rhs._p && lhs._a == rhs._a && lhs._b == rhs._b; }
|
||||
|
||||
/** Returns true if the two sets of parameters are different. */
|
||||
BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
|
||||
|
||||
private:
|
||||
RealType _p;
|
||||
RealType _a;
|
||||
RealType _b;
|
||||
};
|
||||
|
||||
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
|
||||
BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
|
||||
#endif
|
||||
|
||||
/**
|
||||
* Constructs an @c generalized_inverse_gaussian_distribution from its "p", "a", and "b" parameters.
|
||||
*
|
||||
* Requires:
|
||||
* a > 0 && b >= 0 if p > 0,
|
||||
* a > 0 && b > 0 if p == 0,
|
||||
* a >= 0 && b > 0 if p < 0
|
||||
*/
|
||||
explicit generalized_inverse_gaussian_distribution(RealType p_arg = RealType(1.0),
|
||||
RealType a_arg = RealType(1.0),
|
||||
RealType b_arg = RealType(1.0))
|
||||
: _p(p_arg), _a(a_arg), _b(b_arg)
|
||||
{
|
||||
BOOST_ASSERT(
|
||||
(p_arg > RealType(0) && a_arg > RealType(0) && b_arg >= RealType(0)) ||
|
||||
(p_arg == RealType(0) && a_arg > RealType(0) && b_arg > RealType(0)) ||
|
||||
(p_arg < RealType(0) && a_arg >= RealType(0) && b_arg > RealType(0))
|
||||
);
|
||||
init();
|
||||
}
|
||||
/** Constructs an @c generalized_inverse_gaussian_distribution from its parameters. */
|
||||
explicit generalized_inverse_gaussian_distribution(const param_type& parm)
|
||||
: _p(parm.p()), _a(parm.a()), _b(parm.b())
|
||||
{
|
||||
init();
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a random variate distributed according to the
|
||||
* generalized inverse gaussian distribution.
|
||||
*/
|
||||
template<class URNG>
|
||||
RealType operator()(URNG& urng) const
|
||||
{
|
||||
#ifndef BOOST_NO_STDC_NAMESPACE
|
||||
using std::sqrt;
|
||||
using std::log;
|
||||
using std::min;
|
||||
using std::exp;
|
||||
using std::cosh;
|
||||
#endif
|
||||
RealType t = result_type(1);
|
||||
RealType s = result_type(1);
|
||||
RealType log_concave = -psi(result_type(1));
|
||||
if (log_concave >= result_type(.5) && log_concave <= result_type(2)) {
|
||||
t = result_type(1);
|
||||
} else if (log_concave > result_type(2)) {
|
||||
t = sqrt(result_type(2) / (_alpha + _abs_p));
|
||||
} else if (log_concave < result_type(.5)) {
|
||||
t = log(result_type(4) / (_alpha + result_type(2) * _abs_p));
|
||||
}
|
||||
log_concave = -psi(result_type(-1));
|
||||
if (log_concave >= result_type(.5) && log_concave <= result_type(2)) {
|
||||
s = result_type(1);
|
||||
} else if (log_concave > result_type(2)) {
|
||||
s = sqrt(result_type(4) / (_alpha * cosh(1) + _abs_p));
|
||||
} else if (log_concave < result_type(.5)) {
|
||||
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)));
|
||||
}
|
||||
RealType eta = -psi(t);
|
||||
RealType zeta = -psi_deriv(t);
|
||||
RealType theta = -psi(-s);
|
||||
RealType xi = psi_deriv(-s);
|
||||
RealType p = result_type(1) / xi;
|
||||
RealType r = result_type(1) / zeta;
|
||||
RealType t_deriv = t - r * eta;
|
||||
RealType s_deriv = s - p * theta;
|
||||
RealType q = t_deriv + s_deriv;
|
||||
RealType u = result_type(0);
|
||||
RealType v = result_type(0);
|
||||
RealType w = result_type(0);
|
||||
RealType cand = result_type(0);
|
||||
do
|
||||
{
|
||||
u = uniform_01<RealType>()(urng);
|
||||
v = uniform_01<RealType>()(urng);
|
||||
w = uniform_01<RealType>()(urng);
|
||||
if (u < q / (p + q + r)) {
|
||||
cand = -s_deriv + q * v;
|
||||
} else if (u < (q + r) / (p + q + r)) {
|
||||
cand = t_deriv - r * log(v);
|
||||
} else {
|
||||
cand = -s_deriv + p * log(v);
|
||||
}
|
||||
} while (w * chi(cand, s, t, s_deriv, t_deriv, eta, zeta, theta, xi) > exp(psi(cand)));
|
||||
cand = (_abs_p / _omega + sqrt(result_type(1) + _abs_p * _abs_p / (_omega * _omega))) * exp(cand);
|
||||
return _p > 0 ? cand * sqrt(_b / _a) : sqrt(_b / _a) / cand;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a random variate distributed accordint to the beta
|
||||
* distribution with parameters specified by @c param.
|
||||
*/
|
||||
template<class URNG>
|
||||
result_type operator()(URNG& urng, const param_type& parm) const
|
||||
{
|
||||
return generalized_inverse_gaussian_distribution(parm)(urng);
|
||||
}
|
||||
|
||||
/** Returns the "p" parameter of the distribution. */
|
||||
RealType p() const { return _p; }
|
||||
/** Returns the "a" parameter of the distribution. */
|
||||
RealType a() const { return _a; }
|
||||
/** Returns the "b" parameter of the distribution. */
|
||||
RealType b() const { return _b; }
|
||||
|
||||
/** Returns the smallest value that the distribution can produce. */
|
||||
RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
|
||||
{ return RealType(0.0); }
|
||||
/** Returns the largest value that the distribution can produce. */
|
||||
RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
|
||||
{ return (std::numeric_limits<RealType>::infinity)(); }
|
||||
|
||||
/** Returns the parameters of the distribution. */
|
||||
param_type param() const { return param_type(_p, _a, _b); }
|
||||
/** Sets the parameters of the distribution. */
|
||||
void param(const param_type& parm)
|
||||
{
|
||||
_p = parm.p();
|
||||
_a = parm.a();
|
||||
_b = parm.b();
|
||||
init();
|
||||
}
|
||||
|
||||
/**
|
||||
* Effects: Subsequent uses of the distribution do not depend
|
||||
* on values produced by any engine prior to invoking reset.
|
||||
*/
|
||||
void reset() { }
|
||||
|
||||
/** Writes an @c generalized_inverse_gaussian_distribution to a @c std::ostream. */
|
||||
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, generalized_inverse_gaussian_distribution, wd)
|
||||
{
|
||||
os << wd.param();
|
||||
return os;
|
||||
}
|
||||
|
||||
/** Reads an @c generalized_inverse_gaussian_distribution from a @c std::istream. */
|
||||
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, generalized_inverse_gaussian_distribution, wd)
|
||||
{
|
||||
param_type parm;
|
||||
if(is >> parm) {
|
||||
wd.param(parm);
|
||||
}
|
||||
return is;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns true if the two instances of @c generalized_inverse_gaussian_distribution will
|
||||
* return identical sequences of values given equal generators.
|
||||
*/
|
||||
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(generalized_inverse_gaussian_distribution, lhs, rhs)
|
||||
{ return lhs._p == rhs._p && lhs._a == rhs._a && lhs._b == rhs._b; }
|
||||
|
||||
/**
|
||||
* Returns true if the two instances of @c generalized_inverse_gaussian_distribution will
|
||||
* return different sequences of values given equal generators.
|
||||
*/
|
||||
BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(generalized_inverse_gaussian_distribution)
|
||||
|
||||
private:
|
||||
RealType _p;
|
||||
RealType _a;
|
||||
RealType _b;
|
||||
// some data precomputed from the parameters
|
||||
RealType _abs_p;
|
||||
RealType _omega;
|
||||
RealType _alpha;
|
||||
|
||||
/// \cond hide_private_members
|
||||
void init()
|
||||
{
|
||||
#ifndef BOOST_NO_STDC_NAMESPACE
|
||||
using std::abs;
|
||||
using std::sqrt;
|
||||
#endif
|
||||
_abs_p = abs(_p);
|
||||
_omega = sqrt(_a * _b); // two-parameter representation (p, omega)
|
||||
_alpha = sqrt(_omega * _omega + _abs_p * _abs_p) - _abs_p;
|
||||
}
|
||||
|
||||
result_type psi(const RealType& x) const
|
||||
{
|
||||
#ifndef BOOST_NO_STDC_NAMESPACE
|
||||
using std::cosh;
|
||||
using std::exp;
|
||||
#endif
|
||||
return -_alpha * (cosh(x) - result_type(1)) - _abs_p * (exp(x) - x - result_type(1));
|
||||
}
|
||||
|
||||
result_type psi_deriv(const RealType& x) const
|
||||
{
|
||||
#ifndef BOOST_NO_STDC_NAMESPACE
|
||||
using std::sinh;
|
||||
using std::exp;
|
||||
#endif
|
||||
return -_alpha * sinh(x) - _abs_p * (exp(x) - result_type(1));
|
||||
}
|
||||
|
||||
static result_type chi(const RealType& x,
|
||||
const RealType& s, const RealType& t,
|
||||
const RealType& s_deriv, const RealType& t_deriv,
|
||||
const RealType& eta, const RealType& zeta, const RealType& theta, const RealType& xi)
|
||||
{
|
||||
#ifndef BOOST_NO_STDC_NAMESPACE
|
||||
using std::exp;
|
||||
#endif
|
||||
if (x >= -s_deriv && x <= t_deriv) {
|
||||
return result_type(1);
|
||||
} else if (x > t_deriv) {
|
||||
return exp(-eta - zeta * (x - t));
|
||||
}
|
||||
return exp(-theta + xi * (x + s));
|
||||
}
|
||||
/// \endcond
|
||||
};
|
||||
|
||||
} // namespace random
|
||||
|
||||
using random::generalized_inverse_gaussian_distribution;
|
||||
|
||||
} // namespace boost
|
||||
|
||||
#endif // BOOST_GENERALIZED_RANDOM_INVERSE_GAUSSIAN_DISTRIBUTION_HPP
|
||||
@@ -165,6 +165,7 @@ run test_non_central_chi_squared_distribution.cpp /boost/test//boost_unit_test_f
|
||||
run test_hyperexponential.cpp math_test ;
|
||||
run test_hyperexponential_distribution.cpp math_test /boost/test//boost_unit_test_framework ;
|
||||
run test_inverse_gaussian_distribution.cpp /boost/test//boost_unit_test_framework ;
|
||||
run test_generalized_inverse_gaussian_distribution.cpp /boost/test//boost_unit_test_framework ;
|
||||
|
||||
# run nondet_random_speed.cpp ;
|
||||
# run random_device.cpp ;
|
||||
|
||||
@@ -210,6 +210,9 @@ BOOST_AUTO_TEST_CASE(test_generation) {
|
||||
boost::minstd_rand0 gen;
|
||||
BOOST_RANDOM_DISTRIBUTION dist BOOST_RANDOM_TEST1_PARAMS;
|
||||
BOOST_RANDOM_DISTRIBUTION dist_two BOOST_RANDOM_TEST2_PARAMS;
|
||||
#ifdef BOOST_RANDOM_TEST3_PARAMS
|
||||
BOOST_RANDOM_DISTRIBUTION dist_three BOOST_RANDOM_TEST3_PARAMS;
|
||||
#endif
|
||||
typedef BOOST_RANDOM_DISTRIBUTION::result_type result_type;
|
||||
for(int i = 0; i < 10; ++i) {
|
||||
result_type value = dist(gen);
|
||||
@@ -227,6 +230,16 @@ BOOST_AUTO_TEST_CASE(test_generation) {
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST2_MAX
|
||||
BOOST_CHECK_LE(value_two, BOOST_RANDOM_TEST2_MAX);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_PARAMS
|
||||
result_type value_three = dist_three(gen);
|
||||
use(value_three);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_MIN
|
||||
BOOST_CHECK_GE(value_three, BOOST_RANDOM_TEST3_MIN);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_MAX
|
||||
BOOST_CHECK_LE(value_three, BOOST_RANDOM_TEST3_MAX);
|
||||
#endif
|
||||
result_type value_param = dist_two(gen, dist.param());
|
||||
use(value_param);
|
||||
@@ -243,6 +256,16 @@ BOOST_AUTO_TEST_CASE(test_generation) {
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST2_MAX
|
||||
BOOST_CHECK_LE(value_two_param, BOOST_RANDOM_TEST2_MAX);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_PARAMS
|
||||
result_type value_three_param = dist(gen, dist_three.param());
|
||||
use(value_three_param);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_MIN
|
||||
BOOST_CHECK_GE(value_three_param, BOOST_RANDOM_TEST3_MIN);
|
||||
#endif
|
||||
#ifdef BOOST_RANDOM_TEST3_MAX
|
||||
BOOST_CHECK_LE(value_three_param, BOOST_RANDOM_TEST3_MAX);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
@@ -0,0 +1,43 @@
|
||||
/* test_generalized_inverse_gaussian_distribution.cpp
|
||||
*
|
||||
* Distributed under the Boost Software License, Version 1.0. (See
|
||||
* accompanying file LICENSE_1_0.txt or copy at
|
||||
* http://www.boost.org/LICENSE_1_0.txt)
|
||||
*
|
||||
* $Id$
|
||||
*
|
||||
*/
|
||||
|
||||
#include <boost/random/generalized_inverse_gaussian_distribution.hpp>
|
||||
#include <limits>
|
||||
|
||||
#define BOOST_RANDOM_DISTRIBUTION boost::random::generalized_inverse_gaussian_distribution<>
|
||||
#define BOOST_RANDOM_ARG1 p
|
||||
#define BOOST_RANDOM_ARG2 a
|
||||
#define BOOST_RANDOM_ARG3 b
|
||||
#define BOOST_RANDOM_ARG1_DEFAULT 1.0
|
||||
#define BOOST_RANDOM_ARG2_DEFAULT 1.0
|
||||
#define BOOST_RANDOM_ARG3_DEFAULT 1.0
|
||||
#define BOOST_RANDOM_ARG1_VALUE 0.5
|
||||
#define BOOST_RANDOM_ARG2_VALUE 1.5
|
||||
#define BOOST_RANDOM_ARG3_VALUE 2.0
|
||||
|
||||
#define BOOST_RANDOM_DIST0_MIN 0.0
|
||||
#define BOOST_RANDOM_DIST0_MAX (std::numeric_limits<double>::infinity)()
|
||||
#define BOOST_RANDOM_DIST1_MIN 0.0
|
||||
#define BOOST_RANDOM_DIST1_MAX (std::numeric_limits<double>::infinity)()
|
||||
#define BOOST_RANDOM_DIST2_MIN 0.0
|
||||
#define BOOST_RANDOM_DIST2_MAX (std::numeric_limits<double>::infinity)()
|
||||
#define BOOST_RANDOM_DIST3_MIN 0.0
|
||||
#define BOOST_RANDOM_DIST3_MAX (std::numeric_limits<double>::infinity)()
|
||||
|
||||
#define BOOST_RANDOM_TEST1_PARAMS (2.0, 1.0, 1.0)
|
||||
#define BOOST_RANDOM_TEST1_MIN 0.0
|
||||
|
||||
#define BOOST_RANDOM_TEST2_PARAMS (0.5, 0.5, 0.5)
|
||||
#define BOOST_RANDOM_TEST2_MIN 0.0
|
||||
|
||||
#define BOOST_RANDOM_TEST3_PARAMS (2.0, 5.0, 5.0)
|
||||
#define BOOST_RANDOM_TEST3_MIN 0.0
|
||||
|
||||
#include "test_distribution.ipp"
|
||||
Reference in New Issue
Block a user