The distribution function of this generalized asymmetric t-distribution is apparently related to the incomplete Beta function with a transformation $q(x)$.
$$F(x) = B\left(\frac{\nu}{\alpha(1+r^2)},\frac{\nu r}{\alpha(1+r^2)} ; q(x)\right) = \int_0^{q(x)} t^{{\nu}/{(\alpha(1+r^2))}} (1-t)^{{(\nu r)}/{(\alpha(1+r^2))}} dt$$
Since the incomplete Beta function is the distribution function of a beta distributed variable you can sample $q(x)$ from the beta distribution and then compute $x$ from that using the inverse of the formula for $q(x)$.
Below is a derivation of the invers computing $x$ in terms of $q(x)$.
If we simplify the expression by using $\mu = 0$ and $\phi = 1$, and abbreviate some expressions with $d = c^{-\eta}$ and $\eta = \alpha (1+r^2)/r$ then the function $q(x)$ becomes
$$q(r) = \frac{1}{1+d\left(x + \sqrt{1+x^2}\right)^{-\eta}}$$
and the expression $q(r)/(1-q(r))$ that occur in the reverse function are
$$\frac{q(r)}{1-q(r)}= {d\left(x + \sqrt{1+x^2}\right)^{\eta}}$$
or
$$\frac{1}{d}\left(\frac{q(r)}{1-q(r)}\right)^{-\eta}= x + \sqrt{1+x^2}$$
This I can inverse and it becomes
$$x = \frac{1}{2} \left( \frac{1}{d}\left(\frac{q(r)}{1-q(r)}\right)^{-\eta} - \frac{d}{1}\left(\frac{q(r)}{1-q(r)}\right)^{\eta} \right)$$
Note that in the non-generalized case there are also relationships with the beta distribution. If $X = \text{Beta}(1/2,\nu/2)$ and $Y$ is a Bernoulli variable that is either -1 or +1 with equal probability, then $$Y \sqrt{\nu \frac{X}{1-X}} \sim t_\nu$$
And the case $X = \text{Beta}(\nu/2,\nu/2)$ such that
$$\frac{\sqrt{\nu}}{2}\left(\frac{X}{1-X} + \frac{1-X}{X}\right) \sim t_\nu$$
relates to the transformation in Showing $\frac{2X}{1-X^2}$ is standard Cauchy when $X$ is standard Cauchy
This looks similar to the expression used in the article but I cannot get it to works. When I test it for $r=1,\alpha=1,c=1$ then I do not get a t-distribution. It seems to me that they made some typo somewhere. I can make it work when I change a constant
In the code below I use Y = 0.5*sqrt(nu)*(odds^0.5-1/odds^0.5) I stead of Y = 0.5*(odds^0.5-1/odds^0.5) which seems to be what I should use if I follow the formula
$$X = \mu + \frac{\varphi}{2} \left[c^{−1}\left(\frac{q}{1 − q}\right)^\delta − c\left(\frac{q}{1−q}\right)^{-\delta}\right] = \frac{1}{2}\left[\left(\frac{q}{1 − q}\right)^1 − \left(\frac{q}{1−q}\right)^{-1}\right] $$
set.seed(1)
nu = 3
n = 10^4
X = rbeta(n,nu/2,nu/2)
odds = X/(1-X)
Y = 0.5*sqrt(nu)*(odds^0.5-1/odds^0.5)
hist(Y, freq=0, breaks = seq(-1000,1000,0.2), xlim = c(-4,4))
xs=seq(-5,5,0.1)
lines(xs,dt(xs,nu))
This code gives a histogram of the sampled variable and it seems to follow the curve of a t-distribution
