# Closed-form solution or approximation to an infinite series of modified Bessel Functions

## Context:

$$f(x)$$ is the pdf of $$x$$ on a circle with domain $$x \in$$[-$$\pi$$;$$\pi$$], with $$c$$ and $$k$$ being positive real numbers (usually in the range of $$1\div10$$):

$$f(x) = \frac{\exp\bigg(\dfrac{{c \exp\ (k \cos x)}}{2\pi I_0(k)}\bigg)}{Z}$$

This is an exponentiated and weighted (by $$c$$) pdf of the Von Mises distribution. I am trying to figure out a closed-form solution or approximation to $$Z$$, which is the normalizing factor to ensure that $$f(x)$$ integrates to $$1$$. I know I can do it numerically, but for my application speed is important (I want to fit hierarchical bayesian models with this likelihood).

## Work so far

I have been able to derive a power series representation of Z using the series of $$e$$:

\begin{align} Z &= \int_{-\pi}^{\pi} \exp\bigg(\frac{{c \exp\ (k \cos x)}}{2\pi I_0(k)}\bigg) dx \\ &=\int_{-\pi}^{\pi} \sum_{n=0}^{\infty}\bigg(\frac{c}{2\pi I_0(k)}\bigg)^n\frac{1}{n!}\exp(n\ k\cos x) \ dx\\ & = \sum_{n=0}^{\infty}\bigg(\frac{c}{2\pi I_0(k)}\bigg)^n\frac{1}{n!} \int_{-\pi}^{\pi}\exp(n\ k \cos x) \ dx\\ &= 2\pi\sum_{n=0}^{\infty}\bigg(\frac{c}{2\pi I_0(k)}\bigg)^n\frac{I_0(nk)}{n!} \end{align}

I confirmed with a simulation that this is accurate and depending on $$c$$ and $$k$$, I need between $$10$$ and $$30$$ terms of the series for an acceptable estimate of $$Z$$. Unfortunately that takes just as long as the direct numerical integration.

I can expand the modified Bessel functions ($$I_0(\dots))$$ into their series representations, and into their product representations using the zeros of the BesselJ function, but that got me nowhere after several days of trying different things.

We can represent $$Z$$ as a power series of the form:

$$Z = 2\pi\sum_{n=0}^{\infty} \frac{a_n c^n}{n!}$$

where $$a_n = I_0(nk)/{2\pi I_0(k)^n}$$ for all $$n\in\Bbb N$$ and this makes $$Z$$ an exponential generating function of $$a_n$$. Also, note that for the $$n^{th}$$ derivative of $$Z$$ we have$$Z^{(n)}(0) = a_n$$.

## Simplified problem

I also tried to work with a simplified version of f(x):

$$\hat f(x) = \frac{\exp({c \exp\ (k \cos x)})}{\hat Z}$$

This removes one of the modified bessel functions from the series representation of $$\hat Z$$:

$$Z = 2\pi\sum_{n=0}^{\infty}\ \frac{I_0(nk) c^n}{n!} = 2\pi\sum_{n=0}^{\infty}\ \frac{c^n}{n!} \sum_{m=0}^{\infty} \frac{(nk)^{2m}}{m! \ m! \ 2^{2m}}$$

But again I'm stuck here. I also tried expressing the $$I_0$$ functions in terms of sums of $$I_n$$ functions, but also no help.

## Related questions that might help

The following two questions that ask about infinite series besselJ functions seem relevant, but I have not been able to follow them since I don't know anything about contour integration yet:

Specifically, the first question asks for a very related sum:

$$\sum_{n=1}^{\infty} \frac{J_k(nz)}{n^k}$$

The solution there is intriguing, but I don't understand it, and can't borrow it directly since the simplified sum above includes $$c^n/n!$$ instead of $$n^k$$

## Goal:

I want to get to a closed-form solution to either $$Z$$ or $$\hat Z$$. I know this might not exist, but I don't know that it doesn't. Alternatively, a good approximation to either that covers the range of $$c$$ and $$k$$ from $$0$$ to $$\sim20$$ and is fast to compute would be great too. Ultimately, even if a solution does not exist, I'd like to understand how to approach the problem beyond what I already tried, so that I learn something new and useful.

• I learned today about approximation with Chebyshev polynomials <embeddedrelated.com/showarticle/152.php> and it works great. It's not a closed-form solution, but it gives a fixed set of coefficients to evaluate Z for any c and k. Now the question is if I can use that efficiently.
– Ven
Nov 11, 2023 at 16:44
• The honest answer is I don't know - big errors cause the parameter likelihood estimation to fail, but I'm not sure what an acceptable level would be. I could test that using parameter recovery simulations, but I suspect 1% is not too bad (currently using 0.01% percent from the chebyshev approximation with 20 terms.)
– Ven
Nov 14, 2023 at 7:20
• I think you might attract more potential answer writers to this question if one of the tags were [approximation]. While [probability-distributions] is appropriate given the context, it probably does not help much with attracting answers to the question at hand. If you have come up with a good solution yourself in the meantime, please do not hesitate to self-answer your question so we all can learn. Thanks! Nov 15, 2023 at 1:39

If moderate accuracy is sufficient, one can retain the series representation of $$Z$$ proposed by the asker and focus on lowering the computational cost.

$$2 \pi \sum_{n=0}^{\infty}\left( \frac{c}{2 \pi I_{0}(k)} \right)^{n} \frac{I_{0} (nk)}{n!}$$

Note that the above requires IEEE-754 double-precision computation if $$n \in [0, 30]$$ as envisioned by the asker, as the computation of $$I_{0}$$ will cause overflow in intermediate computation otherwise. Clearly, we would want to compute $$\frac{c}{2 \pi I_{0}(k)}$$ once and then just multiply by it for each term of the sum to affect exponentiation. Also, we would want to pre-compute the $$\frac{1}{n!}$$ and just multiply with the appropriate value for each term.

This leaves the computation of $$I_{0}(nk)$$ as the most expensive part of the computation. We would want to compute this cheaply from a simpler function, with moderate accuracy. General shape and series expansions suggest that $$\cosh$$ would be a good candidate, e.g. $$I_{0}(x) \approx \mathrm{R}_{m,n}(x^{2}) \cosh x$$, where $$\mathrm{R}_{m,n}$$ is the quotient of a polynomial of degree $$m$$ and a polynomial of degree $$n$$.

Using $$\cosh$$ as the basic building block has the further advantage that most of the values $$\cosh (nk)$$ can be computed cheaply from previous values via the multi-angle formulas for $$\cosh (2x)$$, $$\cosh (3x)$$, $$\cosh (4x)$$, and $$\cosh (5x)$$. To maintain accuracy and for best performance, fused multiply-add (FMA) operations shall be used throughout. These are available on all major processor architectures, in particular x86-64 and ARM64, and are exposed in C and C++ via a standard math function fma().

J. Olivares, P. Martin, and E. Valero, "A simple approximation for the modified Bessel function of zero order $$I_{0}(x)$$," Journal of Physics: Conference Series 1043 (2018), proposes the following elegant approximation

$$I_{0}(x) \approx \frac{\cosh x}{\sqrt[4]{\left(1+\left(\frac{1}{2}x\right)^2\right)}} \frac{1 + 0.24273 \ x^2}{1 + 0.43023 \ x^2}$$

Using this approximation to compute the $$I_{0}(nk)$$, I observe a maximum relative error in $$Z$$ of $$4.4 \cdot 10^{-3}$$ for $$c, k$$ in $$[\frac{2}{3}, \frac{32}{3}]$$. As an alternative I tried a straightforward computation $$I_{0}(x) \approx \mathrm{R}_{6,6} \cosh x$$ where $$\mathrm{R}_{6,6}$$ is a minimax approximation over $$[0, 1000]$$. With this I observe a maximum relative error in $$Z$$ of $$6.6 \cdot 10^{-4}$$ for $$c, k$$ in $$[\frac{2}{3}, \frac{32}{3}]$$. The numerator of this particular $$\mathrm{R}_{6,6}$$ is

$$3.069864 \cdot 10^{-21} \ x^{12} + 4.711973 \cdot 10^{-15} \ x^{10} + 4.246135 \cdot 10^{-10} \ x^{8} + 3.362724 \cdot 10^{-6} \ x^{6} + 2.486299 \cdot 10^{-3} \ x^{4} + 1.735820 \cdot 10^{-1} \ x^{2} + 1.000000$$

while the denominator is

$$1.891267 \cdot 10^{-19} \ x^{12} + 1.296936 \cdot 10^{-13} \ x^{10} + 6.274928 \cdot 10^{-9} \ x^{8} + 2.735696 \cdot 10^{-5} \ x^{6} + 1.114439 \cdot 10^{-2} \ x^{4} + 4.206092 \cdot 10^{-1} \ x^{2} + 1.000673$$

• This seems very promising! I actually came up across the same paper a few days ago. So far I have one approximation with Chebyshev polynomials, and one with the formula from the Olivares paper, plugging it in a linear regression with several instances of I0 with c and k as x. Both have similarly low error as well. I will compare the efficiency of my two methods and the one you provided!
– Ven
Nov 15, 2023 at 11:24