How to code Mittag-Leffler in Stan?

Modeling
1

I want to write Mittag-Leffler in Stan. I found the following code.

However, the code will not work if Z<0, because log(Z) does not exist for negative Z. Actually, for my data, Z is negative. I also found the Matlab code about Mittag-Leffler from the following website.

However, in function E = LTInversion(t,lambda,alpha,beta,gama,log_epsilon)
…
s_star = abs(lambda)^(1/alpha) * exp(1i*(theta+2k_vettpi)/alpha) ;

I noticed that there is a complex number 1i, I don’t know how to deal with a complex number in Stan.
Similarly, for an R package from the following website, it also has 1i.

Does anyone know how to code Mittag-Leffler in Stan? Or can I use existing Mittag-Leffler R package in Stan code? Thanks!

2

The Stan language does not support complex numbers. You need a numerical algorithm that is specialized for real, but possibly negative, values.

3

Thanks for telling me that Stan does not support complex number. I searched on the website, it seems that the algorithm of Mittag-Leffler function in both Matlab and R packages have 1i (complex number). I cannot follow their algorithm and reproduce that in Stan. Do you know how can I code Mittag-Leffler in Stan for real negative number? Thanks!

4

From Wikipedia:

In the case α and β are real and positive, the series converges for all values of the argument z

If it not converges for all values, eg. in the complex case, then the simulation is not guaranteed to terminate.

5

Thanks for your information that if α and β are real and positive, the series converges for all values of the argument z. I coded a small program in R, because I can compare the results generated from my program with results generated from MittagLefflerR package. If I make my program works in R, and I can convert it to Stan language. Here is my code.

library (MittagLeffleR)
mlf_my<- function (z,a,b)
{
N=10000
tolerance=1e-15
temp = 0
sum=1/gamma(b)
last=sum
for (i in 1:N){
now= z^i/(gamma(a*i+b))
sum=sum+now

if (abs(now)< tolerance && abs(now/last)< 1)
  break

last=now

}

return (sum)
}

my_array<- c(-1,-2,-3,-4,-5,-6,-7,-8,-9,-10)
a=0.5
b=0.5

result_from_mlf_my<- rep(0,10)
result_from_R_library ← rep(0,10)

for(i in 1:length(my_array)){
result_from_R_library[i]=mlf(my_array[i],a,b)
}

for(i in 1:length(my_array)){
result_from_mlf_my[i] = mlf_my(my_array[i],a,b)
}

The results are shown below
Error in if (abs(now) < tolerance && abs(now/last) < 1) break :
missing value where TRUE/FALSE needed
In addition: Warning message:
In mlf_my(my_array[i], a, b) : value out of range in ‘gammafn’

result_from_R_library
[1] 0.136606007 0.053398231 0.027186130 0.016191753 0.010666395 0.007530177 0.005589203 0.004308254 0.003420067 0.002779656
result_from_mlf_my
[1] 1.366060e-01 5.339823e-02 2.718613e-02 1.619172e-02 1.132374e-02 4.453422e+01 -3.891892e+06 0.000000e+00 0.000000e+00 0.000000e+00

It seems that my results are the same as those generated from MittagLeffleR when z = -1, -2, -3, -4. My function also generates some numbers when z= -5, -6, -7. But my results are different from those generated by MittagLeffleR. For z=-8, my program stop working because "value out of range in ‘gammafn’ ". Obviously, my program does not work for z = -5 to -10. Can anyone help me fix my program to make it work for all negative real z, with 0<= a<=1 and 0<= b<=1. Thank you very much!

6

If a and b are known, you could re-implement MittagLeffleR:::LTinversion in C++ using complex numbers and the complex-step method of obtaining the derivative. If you are estimating a and b, then you have the additional problem of getting the partial derivatives with respect to them.