Generating correlated binomial random variables

Question

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)

Answer 1

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)

Answer 2

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

Answer 3

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 correlation 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,]);
[1] 0.300122

mean(SAMPLE[2,]);
[1] 0.800145

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

Answer 4

Here is another R implementation that uses the linear transformation mentioned in the original post.

set.seed(64378)
corr_binomial <- function(n, k, p, rho) {
z <- rbinom(n, k , p)
x_raw <- rbinom(n, k , p)
x <-  z*rho + x_raw*(1 - rho^2)^0.5
observed_corr <- round(cor(x, z), 3)
scatter_z_x <- plot(z, x)
abline(lm(x ~ z))
return(list(observed_corr, summary(x))) 

}

corr_binomial(1000, 50, .5, .6) # Correlated Binomial trials with a correlation of approximately r = .60