grad_reg_lower_inc_gamma.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_LOWER_REG_INC_GAMMA_HPP
#define STAN_MATH_PRIM_FUN_LOWER_REG_INC_GAMMA_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/digamma.hpp>
#include <stan/math/prim/fun/exp.hpp>
#include <stan/math/prim/fun/gamma_p.hpp>
#include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
#include <stan/math/prim/fun/is_any_nan.hpp>
#include <stan/math/prim/fun/is_inf.hpp>
#include <stan/math/prim/fun/lgamma.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log1p.hpp>
#include <stan/math/prim/fun/sqrt.hpp>
#include <stan/math/prim/fun/tgamma.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
#include <limits>
#include <cmath>
 
namespace stan {
namespace math {
 
template <typename T1, typename T2>
return_type_t<T1, T2> grad_reg_lower_inc_gamma(const T1& a, const T2& z,
                                               double precision = 1e-10,
                                               int max_steps = 1e5) {
  using std::exp;
  using std::log;
  using std::pow;
  using std::sqrt;
 
  using TP = return_type_t<T1, T2>;
 
  if (is_any_nan(a, z)) {
    return std::numeric_limits<TP>::quiet_NaN();
  }
 
  check_positive_finite("grad_reg_lower_inc_gamma", "a", a);
 
  if (z == 0.0) {
    return 0.0;
  }
  check_positive_finite("grad_reg_lower_inc_gamma", "z", z);
 
  if ((a < 0.8 && z > 15.0) || (a < 12.0 && z > 30.0)
      || a < sqrt(-756 - value_of_rec(z) * value_of_rec(z)
                  + 60 * value_of_rec(z))) {
    T1 tg = tgamma(a);
    T1 dig = digamma(a);
    return -grad_reg_inc_gamma(a, z, tg, dig, max_steps, precision);
  }
 
  T2 log_z = log(z);
  T2 emz = exp(-z);
 
  int n = 0;
  T1 a_plus_n = a;
  TP sum_a = 0.0;
  T1 lgamma_a_plus_1 = lgamma(a + 1);
  T1 lgamma_a_plus_n_plus_1 = lgamma_a_plus_1;
  TP term;
  while (true) {
    term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1);
    sum_a += term;
    if (term <= precision) {
      break;
    }
    if (n >= max_steps) {
      throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
                         max_steps, "exceeded ",
                         " iterations, gamma_p(a,z) gradient (a) "
                         "did not converge.");
    }
    ++n;
    lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
    ++a_plus_n;
  }
 
  n = 1;
  a_plus_n = a + 1;
  TP digammap1 = digamma(a_plus_n);
  TP sum_b = digammap1 * exp(a * log_z - lgamma_a_plus_1);
  lgamma_a_plus_n_plus_1 = lgamma_a_plus_1 + log(a_plus_n);
  while (true) {
    digammap1 += 1 / a_plus_n;
    term = exp(a_plus_n * log_z - lgamma_a_plus_n_plus_1) * digammap1;
    sum_b += term;
    if (term <= precision) {
      return emz * (log_z * sum_a - sum_b);
    }
    if (n >= max_steps) {
      throw_domain_error("grad_reg_lower_inc_gamma", "n (internal counter)",
                         max_steps, "exceeded ",
                         " iterations, gamma_p(a,z) gradient (a) "
                         "did not converge.");
    }
    ++n;
    lgamma_a_plus_n_plus_1 += log1p(a_plus_n);
    ++a_plus_n;
  }
}
 
}  // namespace math
}  // namespace stan
 
#endif