Simulate correlated Bernoulli variables where the sample distribution is exactly the same as the population distribution 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.

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

