Generating Random Variables from the Generalized Beta Distribution There is a very nice solution to generating random numbers from the Generalized Distribution of the second kind which can be found here.
There is a more general form of this which was developed by McDonald et. al. I have provided the link to my other question pertaining to this distribution where you can find further details as well as the functional form.
Is there a way to incorporate the parameter $c$ into the solution found in the aforementioned solution? I expect there to be a further transformation which would allow this but I have not made any progress on this front. I am curious to see if the community can assist.
 A: The Generalized Beta Distribution is based on a family of density functions parameterized by numbers $0\le c\le 1,$ $p \gt 0,$ and $q \gt 0,$ given by the formula
$$f(x;c,p,q) = \frac{x^{p-1}(1 - (1-c)x)^{q-1}}{B(p,q)(1 + cx)^{p+q}}$$
for $0 \lt x \lt 1/(1-c).$  The quantity $B(p,q)$ is the Beta function, normalizing $f$ to integrate to unity as required of any density.
How does this arise?  Let $Z$ have a Beta$(p,q)$ distribution, which means it is supported on the interval $[0,1]$ with density
$$f_Z(z; p,q) = \frac{1}{B(p,q)} z^{p-1} (1-z)^{q-1}.$$
Pick $0\le c \le 1$ and transform $Z$ into
$$X = \frac{Z}{1-cZ},$$
so that
$$Z = \frac{X}{1 + cX}$$
for $0 \lt X \lt 1/(1-c).$
(Thus, $X$ is a rational function of $Z$ with a simple zero at $0$ and simple pole at $1/c \ge 1.$  This helps us visualize what the transformation does: it stretches $Z$ out more and more near $1,$ with the amount of stretching controlled by $c.$)
The probability element of $X$ thereby becomes, as an easy calculation shows,
$$f_Z\left(\frac{x}{1+cx};p,q\right)\,\mathrm{d}\left(\frac{x}{1+cx}\right) = f(x;c,p,q)\,\mathrm{d}x.$$
When we further transform this random variable $X$ by taking its $a^\text{th}$ root (evidently for $a\ne 0$) and scaling that by a factor $b \gt 0,$ creating the variable
$$Y = b X^{1/a} = b\left(\frac{Z}{1-cZ}\right)^{1/a}, \tag{*}$$
it follows the generalized Beta distribution with parameters $a,b,c,p,q.$

Therefore, you can generate instances of $Y$ by generating instances of $Z$ using a Beta$(p,q)$ random number generator and transforming them according to $(*).$

A simple R implementation is
rbetagen <- function(n, a, b, c, p, q) {
  z <- rbeta(n, p, q)
  b * (z / (1 - c*z))^(1/a)
}

As usual, the first argument n is the number of independent values to obtain and the remaining arguments are the parameters.
(In practice you ought to check the arguments for validity: in particular, that $0\le c \le 1.$)
As a check, I wrote a function to compute the Generalized Beta density function; created some random sets of values of the five parameters; drew ten thousand values from the corresponding distributions; plotted their histograms; and overlaid (in red) graphs of the density for comparison.  Here is the resulting rogue's gallery:

Printed above each plot are the values of the parameters $a,b,c,p,q.$
The agreements between the simulated distributions and the computed densities across a range of parameters show this analysis is correct and that the code works.  (It does not fully test the code, which may have difficulties for some parameter values, especially for $a,$ $p,$ and $q$ close to zero and $c$ close to $1,$ as well as for large negative $a.$  These are the parameter values creating the highly skewed distributions portrayed as spikes in the figure.)
Appendix
Here is the R code to compute the Generalized Beta Density.
#
# The density of X.
#
dbeta3 <- function(x, c, p, q, log.p = FALSE) {
  # 0 <= c <= 1; p, q positive.
  x <- pmax(0, pmin(1 / (1-c), x))
  value <- (p-1) * log(x) + (q-1) * log(1 - (1-c)*x) - (p + q) * log(1 + c*x) -
    lgamma(p) - lgamma(q) + lgamma(p+q)
  if(!isTRUE(log.p)) value <- exp(value)
  value
}
#
# The full generalized Beta: the density of Y.
#
dbetagen <- function(y, a, b, c, p, q, log.p = FALSE) {
  # b > 0 and a != 0.
  x <- (y/b)^a
  value <- dbeta3(x, c, p, q, log.p)
  log.J <-  log(abs(a)) - log(b) + (1 - 1/a) * log(x)
  if (isTRUE(log.p)) value <- value + log.J else value <- value * exp(log.J)
  value
}

