Problem understanding suggestion to generate RVS I came across this interesting preprint which has the potential for solving a number of problems I'm having with asymmetric t-distributed data. So I tried implementing it.
The pdf and cdf work great, but I am having numerical stability issues when attempting to make draws from the resulting distribution. On pp. 6 the author states:

Random numbers may best be generated from the GAT distribution by generating them from the corresponding Beta distribution, and then transforming so that
$$
X = \mu + \frac{\varphi}{2} \left[c^{−1}\left(\frac{q}{1 − q}\right)^\delta − c\left(\frac{q}{1−q}\right)^{-\delta}\right]
$$

which has me confused. $\mu, \varphi, \delta$ are all well-defined, but I have no idea what $q$ is supposed to represent here, nor where the Beta distribution factors in. The closest I can think is that $q$ represents a sample from said Beta distribution, although I don't see what that is as the distribution makes use of the Beta function, and then the formula for $X$ represents a simple algebraic operation?
If anybody could shed a bit of light on this, I'd be extremely grateful!

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

