ENH: Implement exponentially scaled modified Bessel function Op

[dehorsley]

Before

pt.iv(J, a) / pt.exp(a)

After

pt.ive(J, a)

Context for the issue:

Periodic HSGPs (now in PyMC 🥳 ) require iv(J, a) / exp(a) as normalizing factor. Using this directly could lead to overflows however (as iv grows like exp(x)). Instead ive(nu,x) = iv(nu,x) / exp(abs(x)) exists, and as a bonus cheaper to compute than the combined!

After #403 we already have an implementation of ive in JAX, and it's in scipy, so should be straightforward!

[dehorsley]

Happy to make the PR if it's agreed it's worth it

[maresb]

Woah, this is wild. Literally 5 minutes ago I was trying to normalize a periodicized Gaussian and ran into this exact expression.

I'll defer to Ricardo regarding the implementation, but would be very grateful for a reference for this normalization, either in code or literature. (I'm just screwing around in Mathematica.)

[ricardoV94]

Woah, this is wild. Literally 5 minutes ago I was trying to normalize a periodicized Gaussian and ran into this exact expression.

I'll defer to Ricardo regarding the implementation, but would be very grateful for a reference for this normalization, either in code or literature. (I'm just screwing around in Mathematica.)

Scipy docstrings are not too bad here

[maresb]

Sorry, I should have been more specific. I'm looking to see how else ive arises as a normalization factor, or if the HSGP stuff is normalizing the exact same function as me.

[dehorsley]

@maresb Solin & Särkkä 2014 is where I learnt it form

[maresb]

Awesome, thanks @dehorsley!!!

Indeed I was normalizing exactly (23), so it's great to connect this with the literature.