Suppose I want to generate $X_1, \ldots, X_n$ Bernoulli random variables such that:
$$ Corr(X_i, X_j) = \rho $$
for all $i \neq j$.
I am wondering what method I might be able to use (I will be trying to do this in R). I am thinking of using Couplas, but could anyone guide me how to implement this? Thanks.
Because this correlation matrix is so symmetric, we might try to solve the problem with a symmetric distribution.
One of the simplest that gives sufficient flexibility in varying the correlation is the following. Given $d\ge 2$ variables, define a distribution on the set of $d$-dimensional binary vectors $X$ by assigning probability $q$ to $X=(1,1,\ldots, 1)$, probability $q$ to $X=(0,0,\ldots, 0)$, and distributing the remaining probability $1-2q$ equally among the $d$ vectors having exactly one $1$; thus, each of those gets probability $(1-2q)/d$. Note that this family of distributions depends on just one parameter $0\le q\le 1/2$.
It's easy to simulate from one of these distributions: output a vector of zeros with probability $q$, output a vector of ones with probability $q$, and otherwise select uniformly at random from the columns of the $d\times d$ identity matrix.
All the components of $X$ are identically distributed Bernoulli variables. They all have common parameter
$$p = \Pr(X_1 = 1) = q + \frac{1-2q}{d}.$$
Compute the covariance of $X_i$ and $X_j$ by observing they can both equal $1$ only when all the components are $1$, whence
$$\Pr(X_i=1=X_j) = \Pr(X=(1,1,\ldots,1)) = q.$$
This determines the mutual correlation as
$$\rho = \frac{d^2q - ((d-2)q + 1)^2}{(1 + (d-2)q)(d-1 - (d-2)q)}.$$
Given $d \ge 2$ and $-1/(d-1)\le \rho \le 1$ (which is the range of all possible correlations of any $d$-variate random variable), there is a unique solution $q(\rho)$ between $0$ and $1/2$.
Simulations bear this out. Beginning with a set of $21$ equally-spaced values of $\rho$, the corresponding values of $q$ were computed (for the case $d=8$) and used to generate $10,000$ independent values of $X$. The $\binom{8}{2}=28$ correlation coefficients were computed and plotted on the vertical axis. The agreement is good.
I carried out a range of such simulations for values of $d$ between $2$ and $99$, with comparably good results.
A generalization of this approach (namely, allowing for two, or three, or ... values of the $X_i$ simultaneously to equal $1$) would give greater flexibility in varying $E[X_i]$, which in this solution is determined by $\rho$. That combines the ideas related here with the ones in the fully general $d=2$ solution described at https://stats.stackexchange.com/a/285008/919.
The following R code features a function p to compute $q$ from $\rho$ and $d$ and exhibits a fairly efficient simulation mechanism within its main loop.
#
# Determine p(All zeros) = p(All ones) from rho and d.
#
p <- function(rho, d) {
if (rho==1) return(1/2)
if (rho <= -1/(d-1)) return(0)
if (d==2) return((1+rho)/4)
b <- d-2
(4 + 2*b + b^2*(1-rho) - (b+2)*sqrt(4 + b^2 * (1-rho)^2)) / (2 * b^2 * (1-rho))
}
#
# Simulate a range of correlations `rho`.
#
d <- 8 # The number of variables.
n.sim <- 1e4 # The number of draws of X in the simulation.
rho.limits <- c(-1/(d-1), 1)
rho <- seq(rho.limits[1], rho.limits[2], length.out=21)
rho.hat <- sapply(rho, function(rho) {
#
# Compute the probabilities from rho.
#
qd <- q0 <- p(rho, d)
q1 <- (1 - q0 - qd)
#
# First randomly select three kinds of events: all zero, one 1, all ones.
#
u <- sample.int(3, n.sim, prob=c(q0,q1,qd), replace=TRUE)
#
# Conditionally, when there is to be one 1, uniformly select which
# component will equal 1.
#
k <- diag(d)[, sample.int(d, n.sim, replace=TRUE)]
#
# When there are to be all zeros or all ones, make it so.
#
k[, u==1] <- 0
k[, u==3] <- 1
#
# The simulated values of X are the columns of `k`. Return all d*(d-1)/2 correlations.
#
cor(t(k))[lower.tri(diag(d))]
})
#
# Display the simulation results.
#
plot(rho, rho, type="n",
xlab="Intended Correlation",
ylab="Sample Correlation",
xlim=rho.limits, ylim=rho.limits,
main=paste(d, "Variables,", n.sim, "Iterations"))
abline(0, 1, col="Red", lwd=2)
invisible(apply(rho.hat, 1, function(y)
points(rho, y, pch=21, col="#00000010", bg="#00000004")))
