How to generate 2 correlated Beta random variables

Question

I was wondering if it might be possible to generate 2 correlated $\mathsf{Beta}$ random variables?

In other words, I want to generate two Beta random variables which can be said to have come from two Beta distributions whose correlation is known to be "$\rho =$ some number"?

In R, I tried a simple linear transformation approach to achieve this:

set.seed(0)
X1 = rbeta(1e4, 5, 1) ; X2 = rbeta(1e4, 6, 1) ;  X3 = rbeta(1e4, 7, 1) ; a = -.5

Y1 = X1 + (a*X2) ; Y2 = X2 + (a*X3) ## Y1 and Y2 are meant to be correlated

cor(Y1, Y2)

plot(Y1~Y2, col = densCols(Y1, Y2) ) ; abline(lm(Y1~Y2), lty = 3, lwd = 2)

But the linear transformation approach above is particularly unprincipled. If you change a to .5, then the meant-to-be $\mathsf{Beta}$ random variables go beyond 1:

Answer 1

How can you simulate correlated beta distributed variables?

This question has attracted a lot of views, suggesting (future) readers might still be interested in the titular question. I would like to add an approach, with some examples demonstrating that it works. A general description for how to implement it, as well as an implementation in R, is included at the end.

The approach I will be highlighting here, is by using the probability integral transform, which states that samples from a known, continuous probability distribution can be converted to a uniformly distributed sample between $0$ and $1$, using that original distribution's inverse cumulative distribution function. The resulting (correlated) probabilities can then be used as input for the beta distribution's quantile function.

This method roughly preserves rank correlation between both variables. Linear correlation is not preserved, but this is not a reasonable demand, as for many choices of $\alpha$ and $\beta$, even perfect correlation does not yield a linear relationship between both variables.

---

Equal parameters for both distributions

Below are the results for various values of $\rho$, used as input to create two beta distributed variables $x$ and $y$ with the same shape parameters:

$$\begin{aligned}\\x &\sim\mathsf{Beta}(5, 1) \\ \\ y &\sim\mathsf{Beta}(5, 1) \\\\ \end{aligned}$$

Different parameters for both distributions

Below are the results for various values of $\rho$, used as input to create two beta distributed variables $x$ and $y$ with a different set of shape parameters:

$$\begin{aligned}\\x &\sim\mathsf{Beta}(1, 6) \\ \\ y &\sim\mathsf{Beta}(3, 3) \\\\ \end{aligned}$$

(Click on the images to see the full-resolution version.)

---

How to implement this approach

The process works as follows:

1. Draw samples from a bivariate normal distribution $$\begin{aligned}\begin{pmatrix}a \\ b\end{pmatrix} &\sim \mathcal{N}\!\left(\mathbf{0}, \mathbf{\Sigma} \right), \\\\ \mathbf{\Sigma} &= \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\end{aligned}$$ for some given correlation $\rho$ (e.g., using MASS::mvrnorm).
2. Convert these samples to probabilities, using the inverse cumulative distribution function (pnorm).$^\dagger$
3. Use the resulting (correlated) probabilities to compute quantiles of the target distribution (qbeta).$^\ddagger$

This works because simulating from a multivariate normal distribution allows us to directly specify a covariance matrix, and the probability integral transform allows us to convert those samples to (correlated) probabilities. We then simply ask: "Up till what point does the target distribution equal this probability?"

---

Implementation in R

require("MASS")
require("sfsmisc")

set.seed(2024)
n     <- 1e4
alpha <- c(1, 6)
beta  <- c(3, 3)
rho   <- 0.99
Sigma <- matrix(c(1, rho,
                  rho, 1), nrow = 2, byrow = TRUE)
XY    <- mvrnorm(n, c(0, 0), Sigma)
PIT   <- pnorm(XY)
x     <- qbeta(PIT[, 1], alpha[1], beta[1])
y     <- qbeta(PIT[, 2], alpha[2], beta[2])

plot(y ~ x, col = densCols(y, x), pch = 16, bty = "n", axes = FALSE,
     main = bquote(rho == .(rho)), cex.main = 4)
eaxis(1)
eaxis(2)

---

$^\dagger$: This part is the probability integral transform.
$^\ddagger$: This step can be easily adapted for other distributions than beta.

Answer 2

As others have said, the number of possible answers here is endless. We can ignore the issue of the actual distribution. Obviously, for shift and scale distributions, the regression approach is the easiest method to generate correlated data, but it's not an omnibus. Beta RVs are not shift/scale distributions, as you have clearly noticed, if you add 1 to a Beta random variable, it is no longer Beta distributed. The copula method is the most theoretically sound approach for which there are several papers, code examples, and R packages.

In R, a very interesting method is just to sort arrays. Back to your case, simulate independent Beta random variables stored in a matrix. and for a desired correlation, choose a number of rows proportional to that correlation and jointly sort those rows. To gain a closer match, you can iteratively sort more or less rows until convergence is met. Sorting the vectors maintains the univariate distributional properties.

The below example is by no means an efficient implementation, but rather an illustrative example.

## initial setup
set.seed(123)
n <- 100 
x <- matrix(rbeta(n*2, 2, 2), n, 2)

## user supplied desired correlation
p <- 0.5

## np is the number of rows to sort, npprev and npnext to evaluate convergence 
npprev <- n*p
np <- n*p
xnew <- x
xnew[1:np, ] <- apply(xnew[1:np, ], 2, sort)
pest <- cor(xnew[,1], xnew[ ,2])

## add 1 row to sort if under correlated, subtract 1 row if overcorrelated
npnext <- np + sign(p-pest)

## iterative until exact match OR alternating pattern reached
while (pest != p & npnext != npprev ) {
  npprev <- np
  np <- npnext
  xnew <- x
  xnew[1:np, ] <- apply(xnew[1:np, ], 2, sort)
  pest <- cor(xnew[,1], xnew[ ,2])
  npnext <- np + sign(p-pest)
}