# 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$?

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")