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

Question

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:

• Closed form for an infinite series of Bessel functions
• Closed form of $\sum_{n=1}^{\infty} \frac{J_0(2n)}{n^2}$

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.

Answer 1

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$