# Generating correlated binomial random variables

I was wondering if it might be possible to generate correlated random binomial variables following a linear transformation approach?

Below, I tried something simple in R and it produces some correlation. But I was wondering if there is a principled way to do this?

X1 = rbinom(1e4, 6, .5) ; X2 = rbinom(1e4, 6, .5) ;  X3 = rbinom(1e4, 6, .5) ; a = .5

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

cor(Y1, Y2)

• $Y_1$ and $Y_2$ may be correlated, but they will no longer be Binomial. Example, $X_1 = X_2 = 1$ then $Y_1 = 1.5$ hence the $Y_i$ cannot be Binomial random variables. I would suggest you look into the Multinomial distribution. – knrumsey Jun 12 '17 at 19:37
• The short answer to the question is to seek the keyword copula, which helps in generating dependent variables with fixed margins. – Xi'an Jun 13 '17 at 12:34

Binomial variables are usually created by summing independent Bernoulli variables. Let's see whether we can start with a pair of correlated Bernoulli variables $(X,Y)$ and do the same thing.

Suppose $X$ is a Bernoulli$(p)$ variable (that is, $\Pr(X=1)=p$ and $\Pr(X=0)=1-p$) and $Y$ is a Bernoulli$(q)$ variable. To pin down their joint distribution we need to specify all four combinations of outcomes. Writing $$\Pr((X,Y)=(0,0))=a,$$ we can readily figure out the rest from the axioms of probability: $$\Pr((X,Y)=(1,0))=1-q-a, \\\Pr((X,Y)=(0,1))=1-p-a, \\\Pr((X,Y)=(1,1))=a+p+q-1.$$

Plugging this into the formula for the correlation coefficient $\rho$ and solving gives $$a = (1-p)(1-q) + \rho\sqrt{{pq}{(1-p)(1-q)}}.\tag{1}$$

Provided all four probabilities are non-negative, this will give a valid joint distribution--and this solution parameterizes all bivariate Bernoulli distributions. (When $p=q$, there is a solution for all mathematically meaningful correlations between $-1$ and $1$.) When we sum $n$ of these variables, the correlation remains the same--but now the marginal distributions are Binomial$(n,p)$ and Binomial$(n,q)$, as desired.

### Example

Let $n=10$, $p=1/3$, $q=3/4$, and we would like the correlation to be $\rho=-4/5$. The solution to $(1)$ is $a=0.00336735$ (and the other probabilities are around $0.247$, $0.663$, and $0.087$). Here is a plot of $1000$ realizations from the joint distribution: The red lines indicate the means of the sample and the dotted line is the regression line. They are all close to their intended values. The points have been randomly jittered in this image to resolve the overlaps: after all, Binomial distributions only produce integral values, so there will be a great amount of overplotting.

One way to generate these variables is to sample $n$ times from $\{1,2,3,4\}$ with the chosen probabilities and then convert each $1$ into $(0,0)$, each $2$ into $(1,0)$, each $3$ into $(0,1)$, and each $4$ into $(1,1)$. Sum the results (as vectors) to obtain one realization of $(X,Y)$.

### Code

Here is an R implementation.

#
# Compute Pr(0,0) from rho, p=Pr(X=1), and q=Pr(Y=1).
#
a <- function(rho, p, q) {
rho * sqrt(p*q*(1-p)*(1-q)) + (1-p)*(1-q)
}
#
# Specify the parameters.
#
n <- 10
p <- 1/3
q <- 3/4
rho <- -4/5
#
# Compute the four probabilities for the joint distribution.
#
a.0 <- a(rho, p, q)
prob <- c((0,0)=a.0, (1,0)=1-q-a.0, (0,1)=1-p-a.0, (1,1)=a.0+p+q-1)
if (min(prob) < 0) {
print(prob)
stop("Error: a probability is negative.")
}
#
# Illustrate generation of correlated Binomial variables.
#
set.seed(17)
n.sim <- 1000
u <- sample.int(4, n.sim * n, replace=TRUE, prob=prob)
y <- floor((u-1)/2)
x <- 1 - u %% 2
x <- colSums(matrix(x, nrow=n)) # Sum in groups of n
y <- colSums(matrix(y, nrow=n)) # Sum in groups of n
#
# Plot the empirical bivariate distribution.
#
plot(x+rnorm(length(x), sd=1/8), y+rnorm(length(y), sd=1/8),
pch=19, cex=1/2, col="#00000010",
xlab="X", ylab="Y",
main=paste("Correlation is", signif(cor(x,y), 3)))
abline(v=mean(x), h=mean(y), col="Red")
abline(lm(y ~ x), lwd=2, lty=3)

• Can this approach be extended for generating any number of binary variables? - to fit the given correlation matrix (or maximally close to fit it)? – ttnphns Jun 12 '17 at 21:02
• @ttnphns Yes, but the options explode: the $2^k$ probability table must be determined by $k$ marginal parameters, the sum-to-unity constraint, and (therefore) $2^k-k-1$ additional parameters. Evidently you would have a lot of freedom to select (or constrain) those parameters, according to the multivariate properties you wish to create. Also, a similar approach could be used to generate correlated Binomial variables with different values of their "$n$" parameters. Parvin: I believe "when we sum $n$ of these variables" unambiguously explains what $n$ represents. – whuber Jun 12 '17 at 21:08
• This is a nice result. Just to pick a little on your first sentence. To get a binomial from independent Bernoulli variables don't they need to have the same p? This has no effect on what you did since it is just a motivation for you approach. – Michael R. Chernick Jun 14 '17 at 18:38
• @Michael Thank you--you're quite right. Another reason it has no bearing on the method outlined here is that this method still involves summing Bernoulli variables with a common parameter: the parameter is $p$ for all the $X$ variables and $q$ for all the $Y$ variables. To keep the post reasonably short I did not present histograms of the marginal distributions to demonstrate they are indeed Binomial (but I actually did that in my original analysis just to make sure they were working!). – whuber Jun 14 '17 at 20:15
• @TNick I believe you are incorrect, but if you could show me a good reason--such as an example--I would be glad to consider it. Please note that the "full range of correlations" is not $[-1,1]$: there is a universal lower bound and the upper bound is less than $1$ whenever $p\ne q.$ – whuber Jan 28 '19 at 13:55

Using the method described by whuber in his excellent answer, I have programmed a function that generates pairs of correlated binomial random variables using the standard syntax for distributions in R. You can call this function to generate any desired number of correlated Bernoulli random variables, with specified probabilities prob1 and prob1 and specified corelation corr. Note that the correlation coefficient is the correlation of the individual Bernoulli values that sum to the binomial, not the correlation between the binomial values themselves.

rcorrbinom <- function(n, size = 1, prob1, prob2, corr = 0) {

#Check inputs
if (!is.numeric(n))             { stop('Error: n must be numeric') }
if (length(n) != 1)             { stop('Error: n must be a single number') }
if (as.integer(n) != n)         { stop('Error: n must be a positive integer') }
if (n < 1)                      { stop('Error: n must be a positive integer') }
if (!is.numeric(size))          { stop('Error: n must be numeric') }
if (length(size) != 1)          { stop('Error: n must be a single number') }
if (as.integer(size) != size)   { stop('Error: n must be a positive integer') }
if (size < 1)                   { stop('Error: n must be a positive integer') }
if (!is.numeric(prob1))         { stop('Error: prob1 must be numeric') }
if (length(prob1) != 1)         { stop('Error: prob1 must be a single number') }
if (prob1 < 0)                  { stop('Error: prob1 must be between 0 and 1') }
if (prob1 > 1)                  { stop('Error: prob1 must be between 0 and 1') }
if (!is.numeric(prob2))         { stop('Error: prob2 must be numeric') }
if (length(prob2) != 1)         { stop('Error: prob2 must be a single number') }
if (prob2 < 0)                  { stop('Error: prob2 must be between 0 and 1') }
if (prob2 > 1)                  { stop('Error: prob2 must be between 0 and 1') }
if (!is.numeric(corr))          { stop('Error: corr must be numeric') }
if (length(corr) != 1)          { stop('Error: corr must be a single number') }
if (corr < -1)                  { stop('Error: corr must be between -1 and 1') }
if (corr > 1)                   { stop('Error: corr must be between -1 and 1') }

#Compute probabilities
P00   <- (1-prob1)*(1-prob2) + corr*sqrt(prob1*prob2*(1-prob1)*(1-prob2));
P01   <- 1 - prob1 - P00;
P10   <- 1 - prob2 - P00;
P11   <- P00 + prob1 + prob2 - 1;
PROBS <- c(P00, P01, P10, P11)
if (min(PROBS) < 0)       { stop('Error: corr is not in the allowable range') }

#Generate the output
RAND <- array(sample.int(4, size = n*size, replace = TRUE, prob = PROBS),
dim = c(n, size));
VALS <- array(0, dim = c(2, n, size));
OUT  <- array(0, dim = c(2, n));
for (i in 1:n)    {
for (j in 1:size) {
VALS[1,i,j] <- (RAND[i,j] %in% c(3, 4));
VALS[2,i,j] <- (RAND[i,j] %in% c(2, 4)); }
OUT[1, i]   <- sum(VALS[1,i,]);
OUT[2, i]   <- sum(VALS[2,i,]); }

#Give output
OUT; }


Here is an example of using this function to produce a sample array containing a large number of correlated Bernoulli random variables. We can confirm that, for a large sample, the sampled values have sample means and sample correlation that is close to the specified parameters.

#Set parameters
n     <- 10^6;
PROB1 <- 0.3;
PROB2 <- 0.8;
CORR  <- 0.2;

#Generate sample of correlated Bernoulli random variables
set.seed(1);
SAMPLE <- rcorrbinom(n = n, prob1 = PROB1, prob2 = PROB2, corr = CORR);

#Check the properties of the sample
str(SAMPLE);
num [1:2, 1:10000] 0 1 0 1 1 1 0 0 0 1 ...

mean(SAMPLE[1,]);
 0.300122

mean(SAMPLE[2,]);
 0.800145

cor(SAMPLE[1,], SAMPLE[2,]);
 0.20018

• This version in my archives generates an $n\times 2$ array whose first column is $x$ and second column is $y.$ The last line may be of greatest interest. rbernoulli2 <- function(n, p=1/2, q=1/2, rho=0) { a <- (1-p) * (1-q) + rho * sqrt(p * (1-p) * q * (1-q)); prob <- c(01=1-p-a, 00=a, 11=a+p+q-1, 10=1-q-a); if (!all(prob >= 0)) stop("Invalid arguments ", p, ",", q, ",", rho); i <- sample.int(4, n, replace=TRUE, prob=prob); cbind(c(0,1)[ceiling(i/2)], c(0,1)[i %% 2 + 1]) } – whuber Jul 15 at 4:43
• Thanks whuber --- looks good. – Ben Jul 15 at 5:12

A Python (python3) implementation of @whuber's solution:

import numpy as np

def bernoulli_sample(n=100, p=0.5, q=0.5, rho=0):
p1 = rho * np.sqrt(p * q * (1 - p) * (1 - q)) + (1 - p) * (1 - q)
p2 = 1 - p - p1
p3 = 1 - q - p1
p4 = p1 + p + q - 1
samples = np.random.choice([0, 1, 2, 3], size=n, replace=True, p=[p1, p2, p3, p4])
samples = list(map(lambda x: np.array(tuple(np.binary_repr(x, width=2))).astype(np.int), samples))
return np.array(samples).sum(0)

def gen_correlated_bernoulli(size, n=100, p=0.5, q=0.5, rho=0):
return np.array([bernoulli_sample(n, p, q, rho) for _ in range(size)])