How to generate data from a generalized Dirichlet distribution?

Question

I need to generate data from a generalized Dirichlet distribution in Python to test my model, but I have no idea how can I proceed with that, can anyone guide me?

Answer 1

Because there are several conventions for describing this distribution, let's begin by examining the one used in the question. And since the formula looks a little confusing, it might be clearer to use a specific example.

Taking $n=3$ ought to reveal the pattern. In this case the variable is $(X_1,X_2,X_3).$ It is limited to the simplex where $X_i\ge 0$ and $X_1+X_2+X_3 \le 1.$ There are three alphas and three betas for the parameters. The density is

$$\left[\frac{x_1^{\alpha_1-1}\left(1 - x_1\right)^{\beta_1-(\alpha_2+\beta_2)}}{B(\alpha_1,\beta_1)}\right] \left[\frac{x_2^{\alpha_2-1}\left(1 - x_1-x_2\right)^{\beta_2-(\alpha_3+\beta_3)}}{B(\alpha_2,\beta_2)}\right] \left[\frac{x_3^{\alpha_3-1}\left(1 - x_1-x_2-x_3\right)^{\beta_3-1}}{B(\alpha_3,\beta_3)}\right]$$

Each expression in brackets looks much like a Beta density, suggesting we attempt a change of variables where

$$Y_1 = X_1;\ Y_2(1-X_1)=X_2;\ Y_3(1-(X_1+X_2)) = X_3.$$

It follows

$$\mathrm{d}x_1\,\mathrm{d}x_2\,\mathrm{d}x_3 = (1-y_1)^2(1-y_2)\,\mathrm{d}y_1\,\mathrm{d}y_2\,\mathrm{d}y_3$$

and thus, in terms of the $y_i,$ the density is

$$\left[\frac{y_1^{\alpha_1-1}(1-y_1)^{\beta_1-1}}{B(\alpha_1,\beta_1)}\right] \left[\frac{y_2^{\alpha_2-1}(1-y_2)^{\beta_2-1}}{B(\alpha_2,\beta_2)}\right] \left[\frac{y_3^{\alpha_3-1}(1-y_3)^{\beta_3-1}}{B(\alpha_3,\beta_3)}\right].$$

Because this factors into $n=3$ Beta densities,

$(Y_1,Y_2, \ldots, Y_n)$ is independent with marginal Beta distributions having parameters $(\alpha_1,\beta_1),$ $(\alpha_2,\beta_2), \ldots,$ and $(\alpha_n,\beta_n).$

From such a collection of $(Y_i)$ we can recover a Generalized Dirichlet variable $(X_i)$ recursively via

$$X_1 = Y_1;\ X_2 = Y_2(1-X_1);\ X_3 = Y_3(1-X_1-X_2);\text{ etc.}$$

Thus, the problem is solved by generating univariate Beta values and combining them as shown.

An R implementation is shown below. As a check, it compares several low-order multivariate moments of a large sample (with, in this example, $n=4$) to the theoretical moments. The figure shows good agreement, indicating the implementation is correct. The code will work efficiently with values of $n$ into the thousands (at least) if you suppress the moment calculations (the number of them increases as $n^m$ where $m$ is the order of the moment).

#
# Generation of `n` random variates.
#
rGDirichlet <- function(n, alpha, beta) {
  d <- max(length(alpha), length(beta))
  Y <- matrix(rbeta(d*n, alpha, beta), d) # Each sample is in a column of Y
  q <- rep(0, n)
  for (i in seq_len(d)) {                 # Convert `Y` to `X` (in place)
    Y[i, ] <- Y[i, ] * (1 - q)
    q <- q + Y[i, ]
  }
  return(Y)
}
#
# Multivariate moments.
#
mGDirichlet <- function(r, a, b, use.log=TRUE) {
  d <- sum(r) - cumsum(r)
  z <- sum(lgamma(a+b) + lgamma(a+r) + lgamma(b+d) - 
              (lgamma(a) + lgamma(b) + lgamma(a+b+r+d)))
  if(!isTRUE(use.log)) z <- exp(z) # The `r` moment of a GDirichlet(a, b) variate
  z 
}
#------------------------------------------------------------------------------#
#
# Random draws.
#
set.seed(17)
alpha <- 0.1 + round(rexp(4, 1/3), 1)            # alpha parameters
beta <- 0.1 + round(rexp(length(alpha), 1/3), 1) # beta parameters
X <- rGDirichlet(1e4, alpha, beta)
#
# Check some low moments for agreement with theory.
#
m <- 3 # Maximum total order of the moments
R <- as.matrix(do.call(expand.grid, lapply(seq_along(alpha), function(i) 0:m)))
R <- R[rowSums(R) <= m, ]
rownames(R) <- apply(R, 1, paste0, collapse="")
Log.Moments <- apply(R, 1, function(r) {
  c(Theoretical=mGDirichlet(r, alpha, beta),
    Empirical=log(mean(apply(X, 2, function(x) prod(x^r)))))
})
plot(t(Log.Moments), 
     pch=21, bg=hsv(rowSums(R)/(m+1), .8, 1, .25), cex=1.5,
     main=paste0("Logs of Moments Up to Total Order ", m))
abline(0:1, col="Gray", lwd=2)
mtext(paste0("alpha = (", paste(signif(alpha, 2), collapse= ", "),
            ")  beta = (", paste(signif(beta, 2), collapse=", "), ")"), side=3)

Edit: Here's what @whuber's answer looks like in Python:

def GDirichlet(n, alpha, beta):
    d = max(len(alpha), len(beta))
    my_generator = np.random.default_rng()
    Y = my_generator.beta(np.array(alpha), np.array(beta), (n, d)).reshape((d, n))
    q = 0
    for i in range(0, d):
        Y[i, :] = Y[i, :] * (1 - q)
        q = q + Y[i, :]
    return Y