I think the problem has been that people saw the solution as being either too simple or not properly specified. You don't actually calculate the marginal probabilities ... you specify them. Then the rmvbin
function uses the specification of the marginal probabilities and the joint correlations to do the needed sampling to (on average) give joint distributions that match those specs.
library(bindata)
rmvbin<-rmvbin(100, margprob=rep(.5,50), bincorr=cor.mat)
> str(rmvbin)
num [1:100, 1:50] 0 0 0 1 0 0 0 1 0 0 ...
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : NULL
So to look at the sampling characteristics of this result, you can see what correlation there is with the first column:
Hmisc::describe(apply(rmvbin[,-1], 2, function(col) cor(col, rmvbin[,1]) ) )
apply(rmvbin[, -1], 2, function(col) cor(col, rmvbin[, 1]))
n missing unique Mean .05 .10 .25 .50 .75 .90
49 0 38 0.2009 0.05886 0.09874 0.13309 0.19372 0.25208 0.29723
.95
0.33772
lowest : 0.03508 0.04013 0.08696 0.09874 0.10889
highest: 0.29942 0.32450 0.34653 0.40902 0.46714
So the average correlation under sampling was pretty close to the nominal value 0.2. but it did vary fairly widely.