classify functions for numerical instability

[bob-carpenter]

Flag naive implementations (like incomplete beta derivatives) in the manual.

Ideally, characterize at what points in the domain functions are stable for.

Make sure that iterative functions have finite bounds and throw exceptions if they don't converge.

[betanalpha]

Basically anything that used my fragile implementations of the hypergeometric functions, 3F2 and 2F1.

• src/stan/prob/distributions/univariate/continuous/beta.hpp
• src/stan/prob/distributions/univariate/continuous/student_t.hpp
• src/stan/prob/distributions/univariate/discrete/beta_binomial.hpp
• src/stan/prob/distributions/univariate/discrete/neg_binomial.hpp

Some effort should go into testing my implementation for the gradients of the incomplete gamma, too,

• src/stan/prob/distributions/univariate/continuous/inv_chi_square.hpp
• src/stan/prob/distributions/univariate/continuous/scaled_inv_chi_square.hpp
• src/stan/prob/distributions/univariate/continuous/inv_gamma.hpp

[bob-carpenter]

I'd add

• all the complementary cumulative distribution functions (ccdf) which are now implemented as 1-cdf, for values where cdf is near 1
• operator- when numbers are too close to each other (the above problem is due to this and the naive implementation of ccdfs as 1 - cdf); similarly for operator+ when one input is negative and one positive and they're close in absolute value
• exp when inputs are too large (> 800) or too small (< -800), which lead to overflow to +infinity and underflow to 0 respectively

It'd be nice if we could make some of the functions better behaved on the log scale, because what we really need is log_cdf. The problem is that the derivatives involve 1/cdf, so unless things cancel, we're stuck having no more precision than for cdf itself in terms of gradients.

• the basic cov_matrix and corr_matrix data types when there are too many dimension (maybe this is going to be better now that we're not checking symmetry but just using the lower triangular portions --- that did go in, didn't it?)

[aadler]

Michael, just in case you haven't seen these, they may prove helpful:

• DLMF: Hypergeometric §15.19 Methods of Computation
• Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Poschl-Teller-Ginocchio potential wave functions
• cephes / qfloat / qhy2f1.c

[betanalpha]

It'd be nice if we could make some of the functions better behaved on the log scale, because what we really need is log_cdf. The problem is that the derivatives involve 1/cdf, so unless things cancel, we're stuck having no more precision than for cdf itself in terms of gradients.

Technically the gradients require grad(cdf) / cdc which is another one of those functions that really should be
computed directly (a la the polygamma functions) if we wanted to be really robust.

The strategy would involve someone familiar with

a) Implementing recursion relations (most of these functions have nice analytic continuations)
b) Computing asymptotic expansions or Taylor series (asymptotic expansions seem to be more efficient)
c) Given a series for f(x) or 1 + f(x), finding a good series approximation to log(f(x)) or log(1 + f(x)).
d) Pade approximants and rational functions for developing efficient approximations to ratios of functions.