Details of Mittag-Leffler random variate generation
In MittagLeffleR: Mittag-Leffler Family of Distributions

First type Mittag-Leffler distribution

Random variate generation

For the efficient generation of random variates, we use the following useful fact (see e.g. Theorem 19.1 in @Haubold2011a): A standard $\alpha$-Mittag-Leffler random variable $Y$ has the representation:

$$Y \stackrel{d}{=} X^{1/\alpha} Z$$

where $X$ is standard exponentially distributed, $Z$ is $\alpha$-stable with Laplace Transform $$\mathbf E[\exp(-sZ)] = \exp(-s^\alpha),$$ $X$ and $Z$ are independent, and $\stackrel{d}{=}$ means equality in distribution.

Generating $X$

n <- 5
x <- rexp(n)

Generating $Z$

To generate such random variates $Z$, we use

a <- 0.8
sigma <- (cos(pi*a/2))^(1/a)
z <- stabledist::rstable(n = n, alpha = a, beta = 1, gamma = sigma, delta = 0, pm = 1)

Below are the details of the calculation. We use the parametrization of the stable distribution by @SamorodnitskyTaqqu as it has become standard. For $\alpha \in (0,1)$ and $\alpha \in (1,2)$,

$$\mathbf E[\exp(it Z)] = \exp\left\lbrace -\sigma^\alpha |t|^\alpha \left[1 - i \beta {\rm sgn}t \tan \frac{\pi \alpha}{2}\right] + i a t\right\rbrace$$

As in @thebook, Equation (7.28), set

$$\sigma^\alpha = C \Gamma(1-\alpha) \cos \frac{\pi\alpha}{2},$$

for some constant $C > 0$, set $\beta = 1$, set $a = 0$, and the log-characteristic function becomes

\begin{align} -C \frac{\Gamma(2-\alpha)}{1-\alpha} \cos \frac{\pi\alpha}{2} |t|^\alpha \left[1 - i\, {\rm sgn}(t) \tan \frac{\pi \alpha}{2}\right] \ = -C \Gamma(1-\alpha)|t|^\alpha \left[ \cos \frac{\pi \alpha}{2} - i\,{\rm sgn}(t) \sin \frac{\pi \alpha}{2}\right] \ = -C \Gamma(1-\alpha)|t|^\alpha\left(\exp(-i {\rm sgn}(t) \pi/2)\right)^\alpha \ = -C \Gamma(1-\alpha)(-i |t| {\rm sgn}(t))^\alpha \ = -C \Gamma(1-\alpha)(-it)^\alpha \end{align}

Setting $t = is$ recovers the Laplace transform, and to match the Laplace transform $\exp(-s^\alpha)$ of $Z$, it is necessary that $C \Gamma(1-\alpha) = 1$. But then $\sigma^\alpha = \cos(\pi \alpha/2)$, and we see that

$$Z \sim S(\alpha, \beta, \sigma, a) = S(\alpha, 1, \cos(\pi\alpha/2)^{1/\alpha}, 0)$$