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!