How to generate random samples from GB1? Suppose that $X$ follows the Generalized Beta of the first kind (GB1), so its probability density function is given by
$$
f(x \mid a, b, p, q) = \dfrac{a\,x^{ap - 1}\,(1 - (x/b)^a)^{q-1}}{b^{ap}B(p,q)} \hspace{1cm} 0 < x < b
$$
where $ a > 0 , p > 0$ and $q > 0$ are the shape parameters.
My question is how can I generate pseudo-random samples from $X$?
 A: Look at the formula and compare it to the core part of the PDF for a Beta distribution,
$$g(y\mid \alpha,\beta) \propto y^{\alpha-1}(1-y)^{\beta-1}.$$
Evidently the relationship is
$$y=\left(\frac{x}{b}\right)^a.\tag{1}$$
We may verify that by computing
$$\mathrm{d}y = a\frac{1}{b}\left(\frac{x}{b}\right)^{a-1}\mathrm{d}x,$$
whence
$$\eqalign{
g(y\mid \alpha,\beta) \mathrm{d}y &\propto \left(\left(\frac{x}{b}\right)^a\right)^{\alpha-1}\left(1-\left(\frac{x}{b}\right)^a\right)^{\beta-1}\left(\frac{x}{b}\right)^{a-1}\mathrm{d}x \\
&\propto x^{a(\alpha-1) + a-1}\left(1 - \left(\frac{x}{b}\right)^a\right)^{\beta-1}\mathrm{d}x \\
&\propto f(x\mid a,b,\alpha,\beta).
}$$
This identifies $\alpha$ with $p$ and $\beta$ with $q$.
Solving $(1)$ for $x$ yields
$$x = b (y)^{1/a}.\tag{2}$$
Therefore you may generate realizations of a random variable $X$ with distribution given by $f$ simply by generating realizations of a Beta$(p,q)$ random variable $Y$ and computing $X$ with $(2)$.

As a demonstration, here is a plot of 10,000 realizations of $X$ for $(a,b,p,q)=(1,2,3,4)$ generated according to this recipe.  On it is superimposed the graph of $f$.  Their agreement is excellent.

The code (in R) that produced this plot is readily modified to examine other values of the parameters.
#
# Specify the parameters.
#
a <- 1
b <- 2
p <- 3
q <- 4
#
# Define the PDF of a generalized Beta distribution.
#
f <- function(x, a, b, p, q) {
  a * x^(a*p-1) * (1 - (x/b)^a)^(q-1) / (b^(a*p) * beta(p, q))
}
#
# Generate realizations from a generalized Beta.
#
n <- 1e4
set.seed(17)
y <- rbeta(n, p, q)
x <- b*y^(1/a)
#
# Demonstrate the generation is correct by comparing the histogram
# to the theoretical curve.
#
hist(x, freq=FALSE)
curve(f(x, a, b, p, q), add=TRUE, lwd=2, col="#c00000")

