# How to generate 2 correlated Beta random variables

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:

• There is nothing preventing Y1 or Y2 from being outside $(0,1)$ so they cannot be beta distributed? Commented Jun 13, 2017 at 15:17
• Exactly how do you want the variables to be correlated? There are myriad ways to construct them, so we need further restrictions. A general method is to use a copula; another method would be to select a suitable Dirichlet distribution.
– whuber
Commented Jun 13, 2017 at 15:53
• The interesting paper at arxiv.org/pdf/1406.5881.pdf develops a bivariate distribution with beta marginals, that admits of all correlations in $[-1,1]$. They do not give an simulation algorithm but that shouldn't be to difficult. Commented Jun 13, 2017 at 16:08
• Proportions don't necessarily have Beta distributions. In fact, for a fixed denominator that would be a mathematical impossibility: the proportion is discrete while a Beta variable is continuous. You could immediately adopt the solution I posted yesterday about generating correlated Binomial variates, because those provide nice models of the counts that you really need for the t-test: you can't (validly) do a t-test on proportions alone.
– whuber
Commented Jun 13, 2017 at 16:54
• Indeed, the method I proposed for generating correlated binomials seems to fit your situation perfectly, because it constructs them as sums of correlated binary ("dichotomously scored") random variables. Thus, all you have to do is decide how much correlation you want among the pre- and post-responses, item by item, and you're on your way.
– whuber
Commented Jun 13, 2017 at 17:41

## 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.

• The {faux} package now also contains functions that implement this approach, allowing for the simulation of correlated variables drawn from several different distributions: debruine.github.io/faux/articles/norta.html Commented Jul 13 at 16:39

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)
}