Continuous modified second kind Bessel function

Hello,
I'm a bit unsure whether to open a new issue or comment in one of the already closed issues (see https://github.com/stan-dev/math/issues/23 ), so please redirect me if needed.

I'm interested in using the modified second-kind Bessel function inside Stan and taking gradients with respect to the order parameter nu. My understanding is that it would require changing the signature of the function so nu is real and define efficient derivatives (see #1112).

It turns out that such BesselK function comes up often in population genetic predictions, and I’m not aware of ways to avoid this if, for example, the order nu is a function of your parameter of interest.

I’ve ported my Stan models to Turing / Julia so I could use:

https://github.com/cgeoga/BesselK.jl

Which implements exactly this feature.

I’m now thinking that perhaps one could port the Julia implementation (described in https://arxiv.org/abs/2201.00090) and implement it in the Stan math library.

I’m not very mathy, and I’m unsure how this feature relates to the rest of Bessel functions. From the API perspective, it doesn’t feel elegant to only implement the derivatives to one of the functions.

If that’s not an issue, I’m happy to make a PR provided you give me a bit of help on how to start. It seems like URL links in this wiki-page are broken

https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

such

https://mc-stan.org/math/

Best,
Curro

Start with the Stan Math BesselK implementation and issue #1112, then compare the requested behavior with the Julia BesselK.jl implementation and its linked paper. Determine how gradients with respect to a real order parameter should fit the existing Bessel-function API; done means the feature is implemented with the requested automatic differentiation support and validated by appropriate tests.