Add Boost's `tgamma_lower` and `tgamma(a,z)` functions

[andrjohns]

Description

The Math library currently has the normalised incomplete gamma functions (gamma_p and gamma_q), which return values within [0, 1]. The log of the gamma_p function is used for the gamma_lcdf function, but given that gamma_p normalises the result to [0,1], we lose resolution around the boundaries of the LCDF.

If we add the tgamma(a,z) function (incomplete gamma, not normalised), then we can define a log_gamma_p function which should be more numerically stable.

Where the normalised lower incomplete gamma (gamma_p) is defined:

$$ P(a,z) = \frac{\Gamma(a) - \Gamma(a,z)}{\Gamma(a)} $$

With the tgamma_lower function, we can keep this on the log scale as:

double lgamma_a = lgamma(a);

log_diff_exp(lgamma_a, log(tgamma(a,z)) - lgamma_a;

With an analogous approach for gamma_q and log_gamma_q

Current Version:

v4.4.0

EDIT: The definition should have used tgamma, not tgamma_lower!

[hhau]

Chiming in to say accurately computing the right tail for the gamma LCDF is challenging, and even R's implementation tends to overflow to 1 sooner than you'd hope (particularly noticeable when working with mixtures of gammas). There's a recent paper that I came across when trying to work around this behaviour, and it might be interesting / relevant. Unsure what boost does.