Generating Random Variables from the Generalized Beta Distribution

Question

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.

See https://stats.stackexchange.com/questions/140496/how-to-draw-a-random-sample-from-a-generalized-beta-distribution-of-the-second-k — kjetil b halvorsen, Dec 15 '21 at 22:22
Part of your problem is that you don't get the PDF right: it must be a function of $y^a,$ but you omit the exponent in one of the terms. See https://en.wikipedia.org/wiki/Generalized_beta_distribution#Definition. — whuber, Dec 15 '21 at 23:26

whuber · Answer 1 · 2021-12-16T01:26:19.017

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
}

Generating Random Variables from the Generalized Beta Distribution

1 Answers1

Appendix