# How to code Mittag-Leffler in Stan?

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!

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

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!

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.

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
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.