Log of modified_bessel_first_kind

Developers
1

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.

2

Returning pairs is OK or you can just pass in mutable references for the results.

Cutting over to approximations at some threshold should also be OK.

I donā€™t know how much the higher dimensional generalizations come up or whether the Besselā€™s used anywhere else.

3

Any follow-up ? I would be interested by using real orders in thus functionā€¦

4

It got implemented in Stan Math but is not exposed in the Stan language yet.

5

OK, thanks!
I guess I will be patient ;-)

6

@bgoodri, is there an issue in stan-dev/stan to add it with a signature? New functions are trivial to add to Stanā€”all we need is the doc for the manual and the function signature in stan/lang/function_signatures.h.

7

I canā€™t remember. Likely not. I can do it.

8

I see. Hereā€™s the GitHub issue. Iā€™ll assign myself and up-priortize it so I get to it.

@bgoodri: I see stan/math/prim/scal/fun/log_modified_bessel_first_kind.hpp. Can the Stan signature be:

real log_modified_bessel_first_kind(int v, real z);

Weā€™ll be adding a generic data qualifier as soon as the refactor gets moved, so we could have a qualfiied real variable data real v that would have to resolve to data (i.e., something that doesnā€™t need a derivative).

9

If we are going to expose it, then I should implement the rev version.

10

We can do the two things independently. Itā€™s still better as is than having the user write it, I suspect.

11

OK, although we should doc that the log_ version requires the first argument to be greater than -1; otherwise we would be taking the log of something that might not be positive.

12

What is the best way to implement the derivatives of log_modified_bessel_first_kind? We know what the numerators are:


but those should get calculated simultaneously with I_v\left(z\right) because they share basically all of the computation. Is the best thing we can do to make a copy of the prim version, stick into the rev/ directory and insert things like

if (!is_constant_struct<T1>::value) {
  // calculate digamma, etc.
} 

if (!is_constant_struct<T2>::value) {
  // pretend v is v + 1 and calculate part of log_bessel_first_kind(v + 1, z) 
}
13

That sort of code can live in prim. We do things with shared computations all over the probability libraries and all the univariate ones only live in prim.

Hopefully that helps. Want me to send a specific example?

14

You can use MathJax now, e.g., $\int e^x\,\mathrm{d}x$ becomes \int e^x\,\mathrm{d}x.

15

So do it with operands_and_partials?

16

Good day! Could you put descriptions of your approximating formulas?

1 Like
17

I am working on a custom pdf that uses {\displaystyle {}_{1}F_{1}(a,2a,x)} that can be represented in terms of Bessel Functions as follows

{\displaystyle {}_{1}F_{1}(a,2a,x)=e^{x/2}\,{}_{0}F_{1}\left(;a+{\tfrac {1}{2}};{\tfrac {x^{2}}{16}}\right)=e^{x/2}\left({\tfrac {x}{4}}\right)^{1/2-a}\Gamma \left(a+{\tfrac {1}{2}}\right)I_{a-1/2}\left({\tfrac {x}{2}}\right).}

The current Stan implementation of the Bessel function does not accept non-integer orders.
Ideally, I would have used the {\displaystyle {}_{1}F_{1}} directly, but the function is currently not available in stan either.

Is there a possibility of having this function for users in the near future?

1 Like
18

@andrjohns

19

I think the problem may have been calculating derivatives with respect to non-integer orders. There are a couple of relevant issues with commentary and potential workarounds:

2 Likes
20

For completenes, Iā€™ll also add that I spent some effort on implementing the log of modified bessel of the second time and getting the derivatives right (especially w.r.t. order) proved to be beyond my skill: