3
$\begingroup$

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

$\endgroup$
3
$\begingroup$

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.

Figure

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")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.