The function implements the algorithm proposed by Qaqish (2003) to
generate a random sample of d (=length(p)) correlated binary
variables. The sample is generated based on given marginal
probabilities p of the d variables and their correlation matrix
R. The algorithm starts by generating a data for the first
variable X_1 and generates succesively the data for X_2, ... based
on their conditional probabilities P(X_j|X_[i-1],...,X_1),
j=1,...,d.