As I was saying this morning, I have been working on the logarithm of one of the Bessel functions and have some questions about the implementation. Currently our Bessel functions can only be evaluated at integer orders even though Boost can calculate them with any order, but Boost does not implement log-versions. Here I am only talking about the (logarithm of the) modified Bessel function of the first kind.

The immediate reason for doing so is that line 72 of

https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/prob/von_mises_lpdf.hpp#L72

overflows for kappa > 710. That can be fixed by using `log_sum_exp`

logic, although for large kappa it requires summing many terms. For sufficiently large kappa it is much cheaper to use an asymptotic expansion, which Boost does for kappa > 100 and very small orders. Tiny values of kappa are not really a problem numerically, although there is an asymptotic expansion for them.

The first question is what to do about the derivatives. The derivative with respect to kappa involves the same Bessel function but at order -1. It is somewhat inefficient to evaluate the Bessel function separately for one order and then an adjacent order, when both can be calculated simultaneously. But that is not something that we typically do under `stan/math/prim/scal/fun`

. Is it okay to return a `std::pair<double,double>`

? It would make perfect sense for the version under `stan/math/rev/scal/fun`

to calculate the value and the derivative simultaneously (although that is not what is being done now), but `von_mises_lpdf()`

does not call the `rev`

version.

A second reason for reimplementing this is that we cannot currently do the Von Mises-Fisher generalization

because you have to evaluate the (log of) the Bessel function at order p / 2 - 1 and that would currently only work for even p (the dimensionality of the hypersphere). My general implementation works for any real order greater than or equal to -1, although there are specializations that could be implemented for half-integer orders that only involve summation of a finite number of terms.

I am not sure if it is worth the effort.

The last question is whether this Bessel function comes up in any other statistical contexts? In other words, do we only need to get this right for the von Mises (- Fisher) distribution? The modified Bessel function of the second kind is more common in statistics and is defined in terms of the modified Bessel function of the first kind

but there are different computational issues there that would require a separate function.