log_gamma_q_dgamma.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP
#define STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/digamma.hpp>
#include <stan/math/prim/fun/exp.hpp>
#include <stan/math/prim/fun/fabs.hpp>
#include <stan/math/prim/fun/gamma_p.hpp>
#include <stan/math/prim/fun/gamma_q.hpp>
#include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
#include <stan/math/prim/fun/inv.hpp>
#include <stan/math/prim/fun/lgamma.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log1m.hpp>
#include <stan/math/prim/fun/tgamma.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
#include <cmath>
 
namespace stan {
namespace math {
 
namespace internal {
 
constexpr double LOG_Q_GAMMA_CF_PRECISION = 1.49012e-12;
 
template <typename T_a, typename T_z>
inline return_type_t<T_a, T_z> log_q_gamma_cf(const T_a& a, const T_z& z,
                                              double precision
                                              = LOG_Q_GAMMA_CF_PRECISION,
                                              int max_steps = 250) {
  using T_return = return_type_t<T_a, T_z>;
  const T_return log_prefactor = a * log(z) - z - lgamma(a);
 
  T_return b_init = z + 1.0 - a;
  T_return C = (fabs(value_of_rec(b_init)) >= EPSILON)
                   ? b_init
                   : std::decay_t<decltype(b_init)>(EPSILON);
  T_return D = 0.0;
  T_return f = C;
  for (int i = 1; i <= max_steps; ++i) {
    T_a an = -i * (i - a);
    const T_return b = b_init + 2.0 * i;
    D = b + an * D;
    D = (fabs(value_of_rec(D)) >= EPSILON) ? D
                                           : std::decay_t<decltype(D)>(EPSILON);
    C = b + an / C;
    C = (fabs(value_of_rec(C)) >= EPSILON) ? C
                                           : std::decay_t<decltype(C)>(EPSILON);
    D = inv(D);
    const T_return delta = C * D;
    f *= delta;
    const double delta_m1 = fabs(value_of_rec(delta) - 1.0);
    if (delta_m1 < precision) {
      break;
    }
  }
  return log_prefactor - log(f);
}
 
}  // namespace internal
 
template <typename T_a, typename T_z>
inline std::pair<return_type_t<T_a, T_z>, return_type_t<T_a, T_z>>
log_gamma_q_dgamma(const T_a& a, const T_z& z,
                   double precision = internal::LOG_Q_GAMMA_CF_PRECISION,
                   int max_steps = 250) {
  using T_return = return_type_t<T_a, T_z>;
  const double a_val = value_of(a);
  const double z_val = value_of(z);
  // For z > a + 1, use continued fraction for better numerical stability
  if (z_val > a_val + 1.0) {
    std::pair<T_return, T_return> result{
        internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0};
    // For gradient, use: d/da log(Q) = (1/Q) * dQ/da
    // grad_reg_inc_gamma computes dQ/da
    const T_return Q_val = exp(result.first);
    const double dQ_da
        = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val));
    result.second = dQ_da / Q_val;
    return result;
  } else {
    // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy
    const double P_val = gamma_p(a_val, z_val);
    std::pair<T_return, T_return> result{log1m(P_val), 0.0};
    // Gradient: d/da log(Q) = (1/Q) * dQ/da
    // grad_reg_inc_gamma computes dQ/da
    const T_return Q_val = exp(result.first);
    if (Q_val > 0) {
      const double dQ_da
          = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val));
      result.second = dQ_da / Q_val;
    } else {
      // Fallback if Q rounds to zero - use asymptotic approximation
      result.second = log(z_val) - digamma(a_val);
    }
    return result;
  }
}
 
}  // namespace math
}  // namespace stan
 
#endif