I was wondering if it might be possible to generate 2 correlated $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:

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 $Beta$ random variables go beyond 1:

enter image description here

P.S. There is apparently a principled way to do this as shown HERE.

  • 2
    $\begingroup$ There is nothing preventing Y1 or Y2 from being outside $(0,1)$ so they cannot be beta distributed? $\endgroup$ Jun 13, 2017 at 15:17
  • 5
    $\begingroup$ 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. $\endgroup$
    – whuber
    Jun 13, 2017 at 15:53
  • 4
    $\begingroup$ 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. $\endgroup$ Jun 13, 2017 at 16:08
  • 4
    $\begingroup$ 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. $\endgroup$
    – whuber
    Jun 13, 2017 at 16:54
  • 3
    $\begingroup$ 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. $\endgroup$
    – whuber
    Jun 13, 2017 at 17:41

1 Answer 1


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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