Simulate correlated Bernoulli variables where the sample distribution is exactly the same as the population distribution

Question

Is there a way to generate two correlated variables with a Bernoulli distribution where the sample distribution is exactly the same as the population distribution?

I can easily generate a pair of correlated variables with a Bernoulli distribution. Here is an example where the two variables have a correlation of 0.5, and the probability is 0.5.

# simulate normally distributed data

data = MASS::mvrnorm(n=1000, 
                     mu=c(0,0), 
                     Sigma=matrix(data = c(1 , 0.5,
                                          0.5, 1),
                                  nrow = 2,
                                  byrow = TRUE), 
                     empirical=TRUE)

# transform to a Bernoulli distribution

data2 = apply(data,2,function(X){
p <- stats::pnorm(X, 0, 1)
p2<-stats::qbinom(p, 1, 0.5 )
return(p2)
})

cor(data2)
##          [,1]      [,2]
##[1,] 1.0000000 0.3561253
##[2,] 0.3561253 1.0000000

Note though that in this sample, the probability is close, but not exactly 0.5.

sum(data2[,1]/dim(data2)[1])
## [1] 0.507
sum(data2[,2]/dim(data2)[1])
## [1] 0.497

Normally this isn't an issue, and is probably even a feature! Is there a way though to generate correlated variables with a Bernoulli distribution where the sample probability is exactly the same as the population probability? Where the number in each group is exactly what you'd expect? Obviously this is trivial with a single variable (at least when the probability is 0.5).

x1 = sample(rep(c(0,1),1000/2),replace = F,1000)
sum(x1)/1000
##[1] 0.5

But I don't know how I would generate a second variable with the same distribution that has the desired correlation with the first.

Thanks in advance!

Edit While prior questions have asked about simulating Bernoulli data with a known population correlation, this question is about attaining a specific correlation within a sample. I have already accepted an answer, but I think the question is different enough that it is inaccurate to mark it as a duplicate.

Answer 1

As describe in the answer by @jbowman,
you may think of the joint distribution of these two potentially dependent identically distributed Bernoulli variables as a 2-by-2 contingency table:

id	0	1	margin
0	$1-2p+q$	$p-q$	$1-p$
1	$p-q$	$q$	$p$
margin	$1-p$	$p$	$1$

The correlation of those two variables will equal $$ \rho = \frac{ q \cdot (1-p)^2 + (1-2p+q) \cdot p^2 - 2\cdot (p-q) \cdot (1-p)\cdot p}{p \cdot (1-p)} = \frac{q-p^2}{p-p^2}$$ This makes $$ q = \rho \cdot (p - p^2) + p^2 $$ Again, as mentioned by @jbowman, one needs to be very cautious in selection of the sample size. But in some cases, this code may work:

p <- 0.5
ro <- 0.8
q <- ro*(p - p^2) + p^2
joint_distribution <- matrix(c(1-2*p+q, p-q, p-q, q), nrow=2)
n <- 20
sample <- rbind(matrix(rep(c(0,0), round(n * joint_distribution[1,1])), ncol=2, byrow=T),
                matrix(rep(c(1,0), round(n * joint_distribution[2,1])), ncol=2, byrow=T),
                matrix(rep(c(0,1), round(n * joint_distribution[1,2])), ncol=2, byrow=T),
                matrix(rep(c(1,1), round(n * joint_distribution[2,2])), ncol=2, byrow=T))
print(cor(sample[,1], sample[,2]))
print(mean(sample[,1]))
print(mean(sample[,2]))

Output:

[1] 0.8
[1] 0.5
[1] 0.5

Answer 2

If all you care about is generating pairs of Bernoulli variables with a given correlation, note that this can be modeled as a 2x2 table with row labels $(0,1)$ for the first variable and column labels $(0,1)$ for the second variable, indicating whether the two variables equaled $0$ or $1$. You can calculate the exact probabilities for a random draw falling in each of the four cells of the table from the underlying distribution.

However, to find a sample with exactly the same observed frequencies as the underlying probabilities requires a sample size of exactly the right size so that the calculated frequencies can equal the probabilities. If your probabilities are irrational numbers, tough luck! If they are numbers with (hopefully short) finite expansions in some base, e.g., $(0.25, 0.25, 0.25, 0.25)$ in base $10$, then you would need a sample size that is, in this case, a multiple of $4$, and the cell counts would be the probabilities times the sample size. No randomness is involved, so it's really a misnomer to refer to the process as "generating a random sample."