How to speed up calculation of Mittlag-Leffler's function in Stan?


I was coding the log of the Mittag-Leffler function in Stan. Everything fine, except its terrible slow.
Any ideas?

real mlf_log(real z, real a, real b) {
  real s = -lgamma(b);
  real si = s;
  int k = 1;
  real log_z = log(z);
  if(z < 10e-15)
    return s;
  while(si > -35) { // precision limit exp(-35) below 6.305117e-16 stop iteration
     si = k * log_z - lgamma(a * k + b);
     s = log_sum_exp(s, si);
     k += 1;
  return s;


Pretty sure the slow part here is going to be the lgamma (and derivatives) since each of those involves an iterative algorithm as well. Could you calculate the lgamma once and then just update it on each iteration?


Thanks for the quick reply. You’re right the iterative lgamma is slowing down the calculation!
If a is a positive integer, then the lgamma in the loop could replaced by a precalculation
of initial lgamma calculation and log updates in the loop.


If this were coded in c++ directly I think you could also roll over much of the work from one lgamma to the next (at least in derivatives).


Yes, you are right. Thank you!! I forgot about the derivatives:


That’s sophisticating “easy” then. Except I’ve to learn Stans C++ interface.

Added later: Not so easy though, if considering the derivatives of alpha and beta.
No closed form exists in this case.


Yeah, it’s a significant chunk of work. So first question is whether you need it to be faster!


None of this is closed form so that’s ok, it’s just a matter of getting representations that are good enough and carry over as much work as possible iteration-to-iteration.


My implementation is way too slow, about one million seconds per 1000 iterations.
I found out that log_sum_exp calls are not the fastest too.

The Mittag-Leffler function covers many cases whereas the shape parameters are
continuous. I coded already a lot of alternative functions already, but I’m always keen
on to improve my models and learn more.
I used it as a mixture with the Poisson distribution. It’s a very flexible type
of a Hyper-Poisson distribution.

I don’t expect the Stan Team to take any action on this. Maybe it’s a good starter to learn
the Stan C++ API.


Sure, if you want to do it I’m happy to give pointers.