try older version
Varun Agrawal
parent
6fb3f0f209
commit
d614fda81f
|
|
@ -18,7 +18,6 @@
|
||||||
#pragma once
|
#pragma once
|
||||||
|
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
#include <iostream>
|
|
||||||
|
|
||||||
namespace gtsam {
|
namespace gtsam {
|
||||||
|
|
||||||
|
|
@ -103,29 +102,46 @@ static const long double gauss_legendre_50_weights[50] = {
|
||||||
|
|
||||||
namespace internal {
|
namespace internal {
|
||||||
|
|
||||||
|
/// 50 point Gauss-Legendre quadrature
|
||||||
template <typename T>
|
template <typename T>
|
||||||
class IncompleteGamma {
|
constexpr T incomplete_gamma_quad_inp_vals(const T lb, const T ub,
|
||||||
/// 50 point Gauss-Legendre quadrature values
|
|
||||||
static constexpr T quadrature_inp_vals(const T lb, const T ub,
|
|
||||||
const int counter) noexcept {
|
const int counter) noexcept {
|
||||||
return (ub - lb) * gauss_legendre_50_points[counter] / T(2) +
|
return (ub - lb) * gauss_legendre_50_points[counter] / T(2) +
|
||||||
(ub + lb) / T(2);
|
(ub + lb) / T(2);
|
||||||
}
|
}
|
||||||
|
|
||||||
/// 50 point Gauss-Legendre quadrature weights
|
template <typename T>
|
||||||
static constexpr T quadrature_weight_vals(const T lb, const T ub,
|
constexpr T incomplete_gamma_quad_weight_vals(const T lb, const T ub,
|
||||||
const int counter) noexcept {
|
const int counter) noexcept {
|
||||||
return (ub - lb) * gauss_legendre_50_weights[counter] / T(2);
|
return (ub - lb) * gauss_legendre_50_weights[counter] / T(2);
|
||||||
}
|
}
|
||||||
|
|
||||||
static constexpr T quadrature_fn(const T x, const T a,
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_quad_fn(const T x, const T a,
|
||||||
const T lg_term) noexcept {
|
const T lg_term) noexcept {
|
||||||
return exp(-x + (a - T(1)) * log(x) - lg_term);
|
return exp(-x + (a - T(1)) * log(x) - lg_term);
|
||||||
}
|
}
|
||||||
|
|
||||||
static constexpr T quadrature_lb(const T a, const T z) noexcept {
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_quad_recur(const T lb, const T ub, const T a,
|
||||||
|
const T lg_term,
|
||||||
|
const int counter) noexcept {
|
||||||
|
return (counter < 49 ? // if
|
||||||
|
incomplete_gamma_quad_fn(
|
||||||
|
incomplete_gamma_quad_inp_vals(lb, ub, counter), a, lg_term) *
|
||||||
|
incomplete_gamma_quad_weight_vals(lb, ub, counter) +
|
||||||
|
incomplete_gamma_quad_recur(lb, ub, a, lg_term, counter + 1)
|
||||||
|
:
|
||||||
|
// else
|
||||||
|
incomplete_gamma_quad_fn(
|
||||||
|
incomplete_gamma_quad_inp_vals(lb, ub, counter), a, lg_term) *
|
||||||
|
incomplete_gamma_quad_weight_vals(lb, ub, counter));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_quad_lb(const T a, const T z) noexcept {
|
||||||
// break integration into ranges
|
// break integration into ranges
|
||||||
return a > T(1000) ? std::max(T(0), std::min(z, a) - 11 * sqrt(a))
|
return (a > T(1000) ? std::max(T(0), std::min(z, a) - 11 * sqrt(a))
|
||||||
: a > T(800) ? std::max(T(0), std::min(z, a) - 11 * sqrt(a))
|
: a > T(800) ? std::max(T(0), std::min(z, a) - 11 * sqrt(a))
|
||||||
: a > T(500) ? std::max(T(0), std::min(z, a) - 10 * sqrt(a))
|
: a > T(500) ? std::max(T(0), std::min(z, a) - 10 * sqrt(a))
|
||||||
: a > T(300) ? std::max(T(0), std::min(z, a) - 10 * sqrt(a))
|
: a > T(300) ? std::max(T(0), std::min(z, a) - 10 * sqrt(a))
|
||||||
|
|
@ -135,11 +151,12 @@ class IncompleteGamma {
|
||||||
: a > T(50) ? std::max(T(0), std::min(z, a) - 7 * sqrt(a))
|
: a > T(50) ? std::max(T(0), std::min(z, a) - 7 * sqrt(a))
|
||||||
: a > T(40) ? std::max(T(0), std::min(z, a) - 6 * sqrt(a))
|
: a > T(40) ? std::max(T(0), std::min(z, a) - 6 * sqrt(a))
|
||||||
: a > T(30) ? std::max(T(0), std::min(z, a) - 5 * sqrt(a))
|
: a > T(30) ? std::max(T(0), std::min(z, a) - 5 * sqrt(a))
|
||||||
: std::max(T(0), std::min(z, a) - 4 * sqrt(a));
|
: std::max(T(0), std::min(z, a) - 4 * sqrt(a)));
|
||||||
}
|
}
|
||||||
|
|
||||||
static constexpr T quadrature_ub(const T a, const T z) noexcept {
|
template <typename T>
|
||||||
return a > T(1000) ? std::min(z, a + 10 * sqrt(a))
|
constexpr T incomplete_gamma_quad_ub(const T a, const T z) noexcept {
|
||||||
|
return (a > T(1000) ? std::min(z, a + 10 * sqrt(a))
|
||||||
: a > T(800) ? std::min(z, a + 10 * sqrt(a))
|
: a > T(800) ? std::min(z, a + 10 * sqrt(a))
|
||||||
: a > T(500) ? std::min(z, a + 9 * sqrt(a))
|
: a > T(500) ? std::min(z, a + 9 * sqrt(a))
|
||||||
: a > T(300) ? std::min(z, a + 9 * sqrt(a))
|
: a > T(300) ? std::min(z, a + 9 * sqrt(a))
|
||||||
|
|
@ -147,109 +164,92 @@ class IncompleteGamma {
|
||||||
: a > T(90) ? std::min(z, a + 8 * sqrt(a))
|
: a > T(90) ? std::min(z, a + 8 * sqrt(a))
|
||||||
: a > T(70) ? std::min(z, a + 7 * sqrt(a))
|
: a > T(70) ? std::min(z, a + 7 * sqrt(a))
|
||||||
: a > T(50) ? std::min(z, a + 6 * sqrt(a))
|
: a > T(50) ? std::min(z, a + 6 * sqrt(a))
|
||||||
: std::min(z, a + 5 * sqrt(a));
|
: std::min(z, a + 5 * sqrt(a)));
|
||||||
}
|
}
|
||||||
|
|
||||||
static constexpr T quadrature(const T a, const T z) noexcept {
|
template <typename T>
|
||||||
T lb = quadrature_lb(a, z);
|
constexpr T incomplete_gamma_quad(const T a, const T z) noexcept {
|
||||||
T ub = quadrature_ub(a, z);
|
return incomplete_gamma_quad_recur(incomplete_gamma_quad_lb(a, z),
|
||||||
T lg_term = lgamma(a);
|
incomplete_gamma_quad_ub(a, z), a,
|
||||||
T value = quadrature_fn(quadrature_inp_vals(lb, ub, 49), a, lg_term) *
|
lgamma(a), 0);
|
||||||
quadrature_weight_vals(lb, ub, 49);
|
|
||||||
for (size_t counter = 48; counter >= 0; counter--) {
|
|
||||||
value += quadrature_fn(quadrature_inp_vals(lb, ub, counter), a, lg_term) *
|
|
||||||
quadrature_weight_vals(lb, ub, counter);
|
|
||||||
}
|
|
||||||
return value;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
// reverse cf expansion
|
||||||
* @brief Reverse continued fraction expansion
|
// see: https://functions.wolfram.com/GammaBetaErf/Gamma2/10/0003/
|
||||||
* See: https://functions.wolfram.com/GammaBetaErf/Gamma2/10/0003/
|
|
||||||
*
|
template <typename T>
|
||||||
* @param a Degrees of freedom
|
constexpr T incomplete_gamma_cf_2_recur(const T a, const T z,
|
||||||
* @param z
|
|
||||||
* @param depth
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T cf_2_recur(const T a, const T z,
|
|
||||||
const int depth) noexcept {
|
const int depth) noexcept {
|
||||||
if (depth < 100) {
|
return (depth < 100 ? (1 + (depth - 1) * 2 - a + z) +
|
||||||
return (1 + (depth - 1) * 2 - a + z) +
|
depth * (a - depth) /
|
||||||
depth * (a - depth) / cf_2_recur(a, z, depth + 1);
|
incomplete_gamma_cf_2_recur(a, z, depth + 1)
|
||||||
} else {
|
: (1 + (depth - 1) * 2 - a + z));
|
||||||
return 1 + (depth - 1) * 2 - a + z;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Lower (regularized) incomplete gamma function
|
constexpr T incomplete_gamma_cf_2(
|
||||||
*
|
const T a,
|
||||||
* @param a Degrees of freedom
|
const T z) noexcept { // lower (regularized) incomplete gamma function
|
||||||
* @param z
|
return (T(1.0) - exp(a * log(z) - z - lgamma(a)) /
|
||||||
* @return constexpr T
|
incomplete_gamma_cf_2_recur(a, z, 1));
|
||||||
*/
|
|
||||||
static constexpr T continued_fraction_2(const T a, const T z) noexcept {
|
|
||||||
return T(1.0) - exp(a * log(z) - z - lgamma(a)) / cf_2_recur(a, z, 1);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
// cf expansion
|
||||||
* @brief Continued fraction expansion of Gamma function
|
// see: http://functions.wolfram.com/GammaBetaErf/Gamma2/10/0009/
|
||||||
* See: http://functions.wolfram.com/GammaBetaErf/Gamma2/10/0009/
|
|
||||||
*
|
template <typename T>
|
||||||
* @param a Degrees of freedom
|
constexpr T incomplete_gamma_cf_1_coef(const T a, const T z,
|
||||||
* @param z
|
|
||||||
* @param depth
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T cf_1_recur(const T a, const T z,
|
|
||||||
const int depth) noexcept {
|
const int depth) noexcept {
|
||||||
if (depth < 55) {
|
return (is_odd(depth) ? -(a - 1 + T(depth + 1) / T(2)) * z
|
||||||
T cf_coef = is_odd(depth) ? -(a - 1 + T(depth + 1) / T(2)) * z
|
: T(depth) / T(2) * z);
|
||||||
: T(depth) / T(2) * z;
|
|
||||||
return (a + depth - 1) + cf_coef / cf_1_recur(a, z, depth + 1);
|
|
||||||
} else {
|
|
||||||
return (a + depth - 1);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Lower (regularized) incomplete gamma function
|
constexpr T incomplete_gamma_cf_1_recur(const T a, const T z,
|
||||||
*
|
const int depth) noexcept {
|
||||||
* @param a Degrees of freedom
|
return (depth < 55 ? // if
|
||||||
* @param z
|
(a + depth - 1) + incomplete_gamma_cf_1_coef(a, z, depth) /
|
||||||
* @return constexpr T
|
incomplete_gamma_cf_1_recur(a, z, depth + 1)
|
||||||
*/
|
:
|
||||||
static constexpr T continued_fraction_1(const T a, const T z) noexcept {
|
// else
|
||||||
return exp(a * log(z) - z - lgamma(a)) / cf_1_recur(a, z, 1);
|
(a + depth - 1));
|
||||||
}
|
}
|
||||||
|
|
||||||
public:
|
template <typename T>
|
||||||
/**
|
constexpr T incomplete_gamma_cf_1(
|
||||||
* @brief Compute the CDF for the Gamma distribution.
|
const T a,
|
||||||
*
|
const T z) noexcept { // lower (regularized) incomplete gamma function
|
||||||
* @param a Degrees of freedom
|
return (exp(a * log(z) - z - lgamma(a)) /
|
||||||
* @param z
|
incomplete_gamma_cf_1_recur(a, z, 1));
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T compute(const T a, const T z) noexcept {
|
|
||||||
if (is_nan(a) || is_nan(z)) { // NaN check
|
|
||||||
return LIM<T>::quiet_NaN();
|
|
||||||
} else if (a < T(0)) {
|
|
||||||
return LIM<T>::quiet_NaN();
|
|
||||||
} else if (LIM<T>::min() > z) {
|
|
||||||
return T(0);
|
|
||||||
} else if (LIM<T>::min() > a) {
|
|
||||||
return T(1);
|
|
||||||
} else if (a < T(10) && z - a < T(10)) {
|
|
||||||
return continued_fraction_1(a, z);
|
|
||||||
} else if (a < T(10) || z / a > T(3)) {
|
|
||||||
return continued_fraction_2(a, z);
|
|
||||||
} else {
|
|
||||||
return quadrature(a, z);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
//
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_check(const T a, const T z) noexcept {
|
||||||
|
return ( // NaN check
|
||||||
|
(is_nan(a) || is_nan(z)) ? LIM<T>::quiet_NaN() :
|
||||||
|
//
|
||||||
|
a < T(0) ? LIM<T>::quiet_NaN()
|
||||||
|
:
|
||||||
|
//
|
||||||
|
LIM<T>::min() > z ? T(0)
|
||||||
|
:
|
||||||
|
//
|
||||||
|
LIM<T>::min() > a ? T(1)
|
||||||
|
:
|
||||||
|
// cf or quadrature
|
||||||
|
(a < T(10)) && (z - a < T(10)) ? incomplete_gamma_cf_1(a, z)
|
||||||
|
: (a < T(10)) || (z / a > T(3)) ? incomplete_gamma_cf_2(a, z)
|
||||||
|
:
|
||||||
|
// else
|
||||||
|
incomplete_gamma_quad(a, z));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T1, typename T2, typename TC = common_return_t<T1, T2>>
|
||||||
|
constexpr TC incomplete_gamma_type_check(const T1 a, const T2 p) noexcept {
|
||||||
|
return incomplete_gamma_check(static_cast<TC>(a), static_cast<TC>(p));
|
||||||
}
|
}
|
||||||
};
|
|
||||||
|
|
||||||
} // namespace internal
|
} // namespace internal
|
||||||
|
|
||||||
|
|
@ -270,212 +270,208 @@ class IncompleteGamma {
|
||||||
* \frac{\Gamma(a,x)}{\Gamma(a)} = 1 \f] When \f$ a > 10 \f$, a 50-point
|
* \frac{\Gamma(a,x)}{\Gamma(a)} = 1 \f] When \f$ a > 10 \f$, a 50-point
|
||||||
* Gauss-Legendre quadrature scheme is employed.
|
* Gauss-Legendre quadrature scheme is employed.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
template <typename T1, typename T2>
|
template <typename T1, typename T2>
|
||||||
constexpr common_return_t<T1, T2> incomplete_gamma(const T1 a,
|
constexpr common_return_t<T1, T2> incomplete_gamma(const T1 a,
|
||||||
const T2 x) noexcept {
|
const T2 x) noexcept {
|
||||||
using TC = common_return_t<T1, T2>;
|
return internal::incomplete_gamma_type_check(a, x);
|
||||||
return internal::IncompleteGamma<TC>::compute(static_cast<TC>(a),
|
|
||||||
static_cast<TC>(x));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
namespace internal {
|
namespace internal {
|
||||||
|
|
||||||
template <typename T>
|
template <typename T>
|
||||||
class IncompleteGammaInverse {
|
constexpr T incomplete_gamma_inv_decision(const T value, const T a, const T p,
|
||||||
/**
|
const T direc, const T lg_val,
|
||||||
* @brief Compute an initial value for the inverse gamma function which is
|
const int iter_count) noexcept;
|
||||||
* then iteratively updated.
|
|
||||||
*
|
//
|
||||||
* @param a Degrees of freedom
|
// initial value for Halley's method
|
||||||
* @param p
|
template <typename T>
|
||||||
* @return constexpr T
|
constexpr T incomplete_gamma_inv_t_val_1(const T p) noexcept { // a > 1.0
|
||||||
*/
|
return (p > T(0.5) ? sqrt(-2 * log(T(1) - p)) : sqrt(-2 * log(p)));
|
||||||
static constexpr T initial_val(const T a, const T p) noexcept {
|
}
|
||||||
if (a > T(1)) {
|
|
||||||
// Inverse gamma function initial value when a > 1.0
|
template <typename T>
|
||||||
const T t_val = p > T(0.5) ? sqrt(-2 * log(T(1) - p)) : sqrt(-2 * log(p));
|
constexpr T incomplete_gamma_inv_t_val_2(const T a) noexcept { // a <= 1.0
|
||||||
const T sgn_term = p > T(0.5) ? T(-1) : T(1);
|
return (T(1) - T(0.253) * a - T(0.12) * a * a);
|
||||||
const T initial_val_1 =
|
}
|
||||||
t_val -
|
|
||||||
|
//
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_initial_val_1_int_begin(
|
||||||
|
const T t_val) noexcept { // internal for a > 1.0
|
||||||
|
return (t_val -
|
||||||
(T(2.515517L) + T(0.802853L) * t_val + T(0.010328L) * t_val * t_val) /
|
(T(2.515517L) + T(0.802853L) * t_val + T(0.010328L) * t_val * t_val) /
|
||||||
(T(1) + T(1.432788L) * t_val + T(0.189269L) * t_val * t_val +
|
(T(1) + T(1.432788L) * t_val + T(0.189269L) * t_val * t_val +
|
||||||
T(0.001308L) * t_val * t_val * t_val);
|
T(0.001308L) * t_val * t_val * t_val));
|
||||||
const T signed_initial_val_1 = sgn_term * initial_val_1;
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_initial_val_1_int_end(
|
||||||
|
const T value_inp, const T a) noexcept { // internal for a > 1.0
|
||||||
return std::max(
|
return std::max(
|
||||||
T(1e-04),
|
T(1E-04), a * pow(T(1) - T(1) / (9 * a) - value_inp / (3 * sqrt(a)), 3));
|
||||||
a * pow(T(1) - T(1) / (9 * a) - signed_initial_val_1 / (3 * sqrt(a)),
|
|
||||||
3));
|
|
||||||
} else {
|
|
||||||
// Inverse gamma function initial value when a <= 1.0
|
|
||||||
T t_val = T(1) - T(0.253) * a - T(0.12) * a * a;
|
|
||||||
if (p < t_val) {
|
|
||||||
return pow(p / t_val, T(1) / a);
|
|
||||||
} else {
|
|
||||||
return T(1) - log(T(1) - (p - t_val) / (T(1) - t_val));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Compute the error value `f(x)` which we can use for root-finding.
|
constexpr T incomplete_gamma_inv_initial_val_1(
|
||||||
*
|
const T a, const T t_val, const T sgn_term) noexcept { // a > 1.0
|
||||||
* @param value
|
return incomplete_gamma_inv_initial_val_1_int_end(
|
||||||
* @param a Degrees of freedom
|
sgn_term * incomplete_gamma_inv_initial_val_1_int_begin(t_val), a);
|
||||||
* @param p
|
}
|
||||||
* @return constexpr T
|
|
||||||
*/
|
template <typename T>
|
||||||
static constexpr T err_val(const T value, const T a, const T p) noexcept {
|
constexpr T incomplete_gamma_inv_initial_val_2(
|
||||||
|
const T a, const T p, const T t_val) noexcept { // a <= 1.0
|
||||||
|
return (p < t_val ? // if
|
||||||
|
pow(p / t_val, T(1) / a)
|
||||||
|
:
|
||||||
|
// else
|
||||||
|
T(1) - log(T(1) - (p - t_val) / (T(1) - t_val)));
|
||||||
|
}
|
||||||
|
|
||||||
|
// initial value
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_initial_val(const T a, const T p) noexcept {
|
||||||
|
return (a > T(1) ? // if
|
||||||
|
incomplete_gamma_inv_initial_val_1(
|
||||||
|
a, incomplete_gamma_inv_t_val_1(p), p > T(0.5) ? T(-1) : T(1))
|
||||||
|
:
|
||||||
|
// else
|
||||||
|
incomplete_gamma_inv_initial_val_2(
|
||||||
|
a, p, incomplete_gamma_inv_t_val_2(a)));
|
||||||
|
}
|
||||||
|
|
||||||
|
//
|
||||||
|
// Halley recursion
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_err_val(
|
||||||
|
const T value, const T a, const T p) noexcept { // err_val = f(x)
|
||||||
return (incomplete_gamma(a, value) - p);
|
return (incomplete_gamma(a, value) - p);
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Derivative of the incomplete gamma function w.r.t. value
|
constexpr T incomplete_gamma_inv_deriv_1(
|
||||||
*
|
const T value, const T a,
|
||||||
* @param value
|
const T lg_val) noexcept { // derivative of the incomplete gamma function
|
||||||
* @param a Degrees of freedom
|
// w.r.t. x
|
||||||
* @param log_val
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T first_derivative(const T value, const T a,
|
|
||||||
const T lg_val) noexcept {
|
|
||||||
return (exp(-value + (a - T(1)) * log(value) - lg_val));
|
return (exp(-value + (a - T(1)) * log(value) - lg_val));
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Second derivative of the incomplete gamma function w.r.t. value
|
constexpr T incomplete_gamma_inv_deriv_2(
|
||||||
*
|
const T value, const T a,
|
||||||
* @param value
|
const T deriv_1) noexcept { // second derivative of the incomplete gamma
|
||||||
* @param a Degrees of freedom
|
// function w.r.t. x
|
||||||
* @param derivative
|
return (deriv_1 * ((a - T(1)) / value - T(1)));
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T second_derivative(const T value, const T a,
|
|
||||||
const T derivative) noexcept {
|
|
||||||
return (derivative * ((a - T(1)) / value - T(1)));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Compute \f[ \frac{f(x_n)}{f'(x_n)} \f] as part
|
constexpr T incomplete_gamma_inv_ratio_val_1(const T value, const T a,
|
||||||
* of the update denominator.
|
const T p,
|
||||||
*
|
const T deriv_1) noexcept {
|
||||||
* @param value
|
return (incomplete_gamma_inv_err_val(value, a, p) / deriv_1);
|
||||||
* @param a Degrees of freedom
|
|
||||||
* @param p
|
|
||||||
* @param derivative
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T ratio_val_1(const T value, const T a, const T p,
|
|
||||||
const T derivative) noexcept {
|
|
||||||
return (err_val(value, a, p) / derivative);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Compute \f[ \frac{f''(x_n)}{f'(x_n)} \f] as part
|
constexpr T incomplete_gamma_inv_ratio_val_2(const T value, const T a,
|
||||||
* of the update denominator.
|
const T deriv_1) noexcept {
|
||||||
*
|
return (incomplete_gamma_inv_deriv_2(value, a, deriv_1) / deriv_1);
|
||||||
* @param value
|
|
||||||
* @param a Degrees of freedom
|
|
||||||
* @param derivative
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T ratio_val_2(const T value, const T a,
|
|
||||||
const T derivative) noexcept {
|
|
||||||
return (second_derivative(value, a, derivative) / derivative);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Halley's method update delta
|
constexpr T incomplete_gamma_inv_halley(const T ratio_val_1,
|
||||||
*
|
const T ratio_val_2) noexcept {
|
||||||
* @param ratio_val_1
|
|
||||||
* @param ratio_val_2
|
|
||||||
* @return constexpr T
|
|
||||||
*/
|
|
||||||
static constexpr T halley(const T ratio_val_1, const T ratio_val_2) noexcept {
|
|
||||||
return (ratio_val_1 /
|
return (ratio_val_1 /
|
||||||
std::max(T(0.8), std::min(T(1.2), T(1) - T(0.5) * ratio_val_1 *
|
std::max(T(0.8), std::min(T(1.2), T(1) - T(0.5) * ratio_val_1 *
|
||||||
ratio_val_2)));
|
ratio_val_2)));
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
template <typename T>
|
||||||
* @brief Compute the iterative solution for the
|
constexpr T incomplete_gamma_inv_recur(const T value, const T a, const T p,
|
||||||
* incomplete inverse gamma function.
|
const T deriv_1, const T lg_val,
|
||||||
*
|
const int iter_count) noexcept {
|
||||||
* @param initial_val Initial value guess
|
return incomplete_gamma_inv_decision(
|
||||||
* @param a Degrees of freedom
|
value, a, p,
|
||||||
* @param p
|
incomplete_gamma_inv_halley(
|
||||||
* @param lg_val
|
incomplete_gamma_inv_ratio_val_1(value, a, p, deriv_1),
|
||||||
* @return constexpr T
|
incomplete_gamma_inv_ratio_val_2(value, a, deriv_1)),
|
||||||
*/
|
lg_val, iter_count);
|
||||||
static constexpr T find_root(const T initial_val, const T a, const T p,
|
|
||||||
const T lg_val) noexcept {
|
|
||||||
const int GAMMA_INV_MAX_ITER = 35;
|
|
||||||
T x = initial_val;
|
|
||||||
T derivative = first_derivative(initial_val, a, lg_val);
|
|
||||||
for (size_t counter = 1; counter <= GAMMA_INV_MAX_ITER; counter++) {
|
|
||||||
T direc = halley(ratio_val_1(x, a, p, derivative),
|
|
||||||
ratio_val_2(x, a, derivative));
|
|
||||||
derivative = first_derivative(x, a, lg_val);
|
|
||||||
x = x - direc;
|
|
||||||
}
|
|
||||||
return x - halley(ratio_val_1(x, a, p, derivative),
|
|
||||||
ratio_val_2(x, a, derivative));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
public:
|
template <typename T>
|
||||||
/**
|
constexpr T incomplete_gamma_inv_decision(const T value, const T a, const T p,
|
||||||
* @brief Compute the percent point function for the Gamma distribution.
|
const T direc, const T lg_val,
|
||||||
*
|
const int iter_count) noexcept {
|
||||||
* @param a Degrees of freedom
|
// return( abs(direc) > GCEM_INCML_GAMMA_INV_TOL ?
|
||||||
* @param p
|
// incomplete_gamma_inv_recur(value - direc, a, p,
|
||||||
* @return constexpr T
|
// incomplete_gamma_inv_deriv_1(value,a,lg_val), lg_val) : value - direc );
|
||||||
*/
|
#define INCML_GAMMA_INV_MAX_ITER 35
|
||||||
static constexpr T compute(const T a, const T p) noexcept {
|
return (iter_count <= INCML_GAMMA_INV_MAX_ITER ? // if
|
||||||
// Perform checks on the input and return the corresponding best answer
|
incomplete_gamma_inv_recur(
|
||||||
if (is_nan(a) || is_nan(p)) { // NaN check
|
value - direc, a, p,
|
||||||
return LIM<T>::quiet_NaN();
|
incomplete_gamma_inv_deriv_1(value, a, lg_val), lg_val,
|
||||||
} else if (LIM<T>::min() > p) { // Check lower bound
|
iter_count + 1)
|
||||||
return T(0);
|
:
|
||||||
} else if (p > T(1)) { // Check upper bound
|
// else
|
||||||
return LIM<T>::quiet_NaN();
|
value - direc);
|
||||||
} else if (LIM<T>::min() > abs(T(1) - p)) {
|
|
||||||
return LIM<T>::infinity();
|
|
||||||
} else if (LIM<T>::min() > a) { // Check lower bound for degrees of freedom
|
|
||||||
return T(0);
|
|
||||||
} else {
|
|
||||||
return find_root(initial_val(a, p), a, p, lgamma(a));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_begin(const T initial_val, const T a,
|
||||||
|
const T p, const T lg_val) noexcept {
|
||||||
|
return incomplete_gamma_inv_recur(
|
||||||
|
initial_val, a, p, incomplete_gamma_inv_deriv_1(initial_val, a, lg_val),
|
||||||
|
lg_val, 1);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
constexpr T incomplete_gamma_inv_check(const T a, const T p) noexcept {
|
||||||
|
return ( // NaN check
|
||||||
|
(is_nan(a) || is_nan(p)) ? LIM<T>::quiet_NaN() :
|
||||||
|
//
|
||||||
|
LIM<T>::min() > p ? T(0)
|
||||||
|
: p > T(1) ? LIM<T>::quiet_NaN()
|
||||||
|
: LIM<T>::min() > abs(T(1) - p) ? LIM<T>::infinity()
|
||||||
|
:
|
||||||
|
//
|
||||||
|
LIM<T>::min() > a ? T(0)
|
||||||
|
:
|
||||||
|
// else
|
||||||
|
incomplete_gamma_inv_begin(incomplete_gamma_inv_initial_val(a, p), a,
|
||||||
|
p, lgamma(a)));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T1, typename T2, typename TC = common_return_t<T1, T2>>
|
||||||
|
constexpr TC incomplete_gamma_inv_type_check(const T1 a, const T2 p) noexcept {
|
||||||
|
return incomplete_gamma_inv_check(static_cast<TC>(a), static_cast<TC>(p));
|
||||||
}
|
}
|
||||||
};
|
|
||||||
|
|
||||||
} // namespace internal
|
} // namespace internal
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Compile-time inverse incomplete gamma function
|
* Compile-time inverse incomplete gamma function
|
||||||
*
|
*
|
||||||
* Compute the value \f$ x \f$
|
* @param a a real-valued, non-negative input.
|
||||||
* such that \f[ f(x) := \frac{\gamma(a,x)}{\Gamma(a)} - p \f] equal to zero,
|
* @param p a real-valued input with values in the unit-interval.
|
||||||
* for a given \c p.
|
|
||||||
*
|
*
|
||||||
* We find this root using Halley's method:
|
* @return Computes the inverse incomplete gamma function, a value \f$ x \f$
|
||||||
* \f[ x_{n+1} = x_n - \frac{f(x_n)/f'(x_n)}{1 - 0.5 \frac{f(x_n)}{f'(x_n)}
|
* such that \f[ f(x) := \frac{\gamma(a,x)}{\Gamma(a)} - p \f] equal to zero,
|
||||||
* \frac{f''(x_n)}{f'(x_n)} } \f] where
|
* for a given \c p. GCE-Math finds this root using Halley's method: \f[ x_{n+1}
|
||||||
* \f[ \frac{\partial}{\partial x} \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) =
|
* = x_n - \frac{f(x_n)/f'(x_n)}{1 - 0.5 \frac{f(x_n)}{f'(x_n)}
|
||||||
* \frac{1}{\Gamma(a)} x^{a-1} \exp(-x) \f] \f[ \frac{\partial^2}{\partial x^2}
|
* \frac{f''(x_n)}{f'(x_n)} } \f] where \f[ \frac{\partial}{\partial x}
|
||||||
|
* \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1}
|
||||||
|
* \exp(-x) \f] \f[ \frac{\partial^2}{\partial x^2}
|
||||||
* \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1}
|
* \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1}
|
||||||
* \exp(-x) \left( \frac{a-1}{x} - 1 \right) \f]
|
* \exp(-x) \left( \frac{a-1}{x} - 1 \right) \f]
|
||||||
*
|
|
||||||
* @param a The degrees of freedom for the gamma distribution.
|
|
||||||
* @param p The quantile value for computing the percent point function.
|
|
||||||
*
|
|
||||||
* @return Computes the inverse incomplete gamma function.
|
|
||||||
*/
|
*/
|
||||||
|
|
||||||
template <typename T1, typename T2>
|
template <typename T1, typename T2>
|
||||||
constexpr common_return_t<T1, T2> incomplete_gamma_inv(const T1 a,
|
constexpr common_return_t<T1, T2> incomplete_gamma_inv(const T1 a,
|
||||||
const T2 p) noexcept {
|
const T2 p) noexcept {
|
||||||
using TC = common_return_t<T1, T2>;
|
return internal::incomplete_gamma_inv_type_check(a, p);
|
||||||
return internal::IncompleteGammaInverse<TC>::compute(static_cast<TC>(a),
|
|
||||||
static_cast<TC>(p));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue