Hello,

I read Ben Goodrich’s post from last summer on stan-dev about reimplementing the modified bessel function of the first kind (even more specific: its logarithm).

Ben was curious in what statistical contexts this Bessel function comes up.

Well …

Last week I started fitting Stan models with the Skellam distribution, and this needs …. The modified Bessel of the first kind!

The Skellam distribution is the distribution for the difference of two independent Poisson distributed variables.

See e.g. here

It appears well suited for modelling football score differences:

During fitting I ran into the overflow in boost::math::cyl_bessel_i problem.

This was reported on the Stan users list before, e.g. in august 2016 by Mike Lawrence for the von mises dist.

It appears only at the start of a new chain (but not always), the chains that get past this point behave as expected.

For now, I hard coded a ceiling at x = 700:

if (x < 700)

log_prob = total + log(modified_bessel_first_kind(k, x));

else

log_prob = total + log(modified_bessel_first_kind(k, 700));

This doesn’t feel right, but it seems to work.

I’m a bit unsure how to proceed, and would very much appreciate some advice.

Gertjan