grad_reg_lower_inc_gamma()

template<typename T1 , typename T2 >

return_type_t< T1, T2 > stan::math::grad_reg_lower_inc_gamma	(	const T1 &	a,
		const T2 &	z,
		double	precision = 1e-10,
		int	max_steps = 1e5
	)

Computes the gradient of the lower regularized incomplete gamma function.

The lower incomplete gamma function derivative w.r.t its first parameter (a) seems to have no standard source. It also appears to have no widely known approximate implementation. Gautschi (1979) has a thorough discussion of the calculation of the lower regularized incomplete gamma function itself and some stability issues.

Reference: Gautschi, Walter (1979) ACM Transactions on mathematical software. 5(4):466-481

We implemented calculations for d(gamma_p)/da by taking derivatives of formulas suggested by Gauschi and others and testing them against an outside source (Mathematica). We took three implementations which can cover the range {a:[0,20], z:[0,30]} with absolute error < 1e-10 with the exception of values near (0,0) where the error is near 1e-5. Relative error is also <<1e-6 except for regions where the gradient approaches zero.

Gautschi suggests calculating the lower incomplete gamma function for small to moderate values of $z$ using the approximation:

\[ \frac{\gamma(a,z)}{\Gamma(a)}=z^a e^-z \sum_n=0^\infty \frac{z^n}{\Gamma(a+n+1)} \]

We write the derivative in the form:

\[ \frac{d\gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z} \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)} - \frac{1}{e^z} \sum_n=0^\infty \frac{z^{a+n}}{\Gamma(a+n+1)}\psi^0(a+n+1) \]

This calculation is sufficiently accurate for small $a$ and small $z$. For larger values and $a$ and $z$ we use it in its log form:

\[ \frac{d \gamma(a,z)\Gamma(a)}{da} = \frac{\log z}{e^z} \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1)] - \sum_n=0^\infty \exp[(a+n)\log z - \log\Gamma(a+n+1) + \log\psi^0(a+n+1)] \]

For large $z$, Gauschi recommends using the upper incomplete Gamma instead and the negative of its derivative turns out to be more stable and accurate for larger $z$ and for some combinations of $a$ and $z$. This is a log-scale implementation of the derivative of the formulation suggested by Gauschi (1979). For some values it defers to the negative of the gradient for the gamma_q function. This is based on the suggestion by Gauschi (1979) that for large values of $z$ it is better to carry out calculations using the upper incomplete Gamma function.

Branching for choice of implementation for the lower incomplete regularized gamma function gradient. The derivative based on Gautschi's formulation appears to be sufficiently accurate everywhere except for large z and small to moderate a. The intersection between the two regions is a radius 12 quarter circle centered at a=0, z=30 although both implementations are satisfactory near the intersection.

Some limits that could be treated, e.g., infinite z should return tgamma(a) * digamma(a), throw instead to match the behavior of, e.g., boost::math::gamma_p

Template Parameters

T1	type of a
T2	type of z

Parameters

[in]	a	shared with complete Gamma
[in]	z	value to integrate up to
[in]	precision	series terminates when increment falls below this value.
[in]	max_steps	number of terms to sum before throwing

Exceptions

std::domain_error

if the series does not converge to requested precision before max_steps.

Definition at line 110 of file grad_reg_lower_inc_gamma.hpp.