grad_reg_inc_gamma.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
#define STAN_MATH_PRIM_FUN_GRAD_REG_INC_GAMMA_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/gamma_p.hpp>
#include <stan/math/prim/fun/gamma_q.hpp>
#include <stan/math/prim/fun/is_any_nan.hpp>
#include <stan/math/prim/fun/is_inf.hpp>
#include <stan/math/prim/fun/lmultiply.hpp>
#include <cmath>
#include <limits>
 
namespace stan {
namespace math {
 
template <typename T1, typename T2>
inline return_type_t<T1, T2> grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig,
                                                double precision = 1e-6,
                                                int max_steps = 1e5) {
  using std::exp;
  using std::fabs;
  using std::log;
  using TP = return_type_t<T1, T2>;
 
  if (is_any_nan(a, z, g, dig)) {
    return std::numeric_limits<TP>::quiet_NaN();
  }
 
  T2 l = log(z);
  if (z >= a && z >= 8) {
    // asymptotic expansion http://dlmf.nist.gov/8.11#E2
    TP S = 0;
    T1 a_minus_one_minus_k = a - 1;
    T1 fac = a_minus_one_minus_k;  // falling_factorial(a-1, k)
    T1 dfac = 1;                   // d/da[falling_factorial(a-1, k)]
    T2 zpow = z;                   // z ** k
    TP delta = dfac / zpow;
 
    for (int k = 1; k < 10; ++k) {
      a_minus_one_minus_k -= 1;
 
      S += delta;
 
      zpow *= z;
      dfac = a_minus_one_minus_k * dfac + fac;
      fac *= a_minus_one_minus_k;
      delta = dfac / zpow;
 
      if (is_inf(delta)) {
        throw_domain_error("grad_reg_inc_gamma", "is not converging", "", "");
      }
    }
 
    return gamma_q(a, z) * (l - dig) + exp(-z + (a - 1) * l) * S / g;
  } else {
    // gradient of series expansion http://dlmf.nist.gov/8.7#E3
    TP S = 0;
    TP log_s = 0.0;
    double s_sign = 1.0;
    T2 log_z = log(z);
    TP log_delta = log_s - lmultiply(2, a);
    for (int k = 1; k <= max_steps; ++k) {
      S += s_sign >= 0.0 ? exp(log_delta) : -exp(log_delta);
      log_s += log_z - log(k);
      s_sign = -s_sign;
      log_delta = log_s - lmultiply(2, k + a);
      if (is_inf(log_delta)) {
        throw_domain_error("grad_reg_inc_gamma", "is not converging", "", "");
      }
      if (log_delta <= log(precision)) {
        return gamma_p(a, z) * (dig - l) + exp(a * l) * S / g;
      }
    }
    throw_domain_error(
        "grad_reg_inc_gamma", "k (internal counter)", max_steps, "exceeded ",
        " iterations, gamma function gradient did not converge.");
    return INFTY;
  }
}
 
}  // namespace math
}  // namespace stan
#endif