3
$\begingroup$

Suppose I have $M$ Bernoulli distributions with parameters $p_i$, pairwise correlation $\rho_{ij}$ for $i\neq j$. I would like to generate $N$ samples from the joint distribution. The case of $M=2$ has been solved here. A similar question in R was asked here. The answers use bindata package from R.

I know that this distribution is not unique and is not always solvable. This paper describes one approach to this problem.

My question is about the general algorithm of the paper, not a specific package.

I am implementing the 3 x 3 example described on page 7. I know this specific case is viable via the algorithm above because the paper has this example. Here is my approach to implement the algorithm:

  1. Compute the quantiles for $\mu_1 = Q(p_1), \mu_2=Q(p_2)$ and $\mu_3=Q(p_3)$.
  2. For $1\leq i\neq j\leq 3$, compute $p_{ij}$ values from the $L(-\mu_i, -\mu_j, \rho_{ij})$. From Figure 2, it seems like my $p_{ij}$ values were in line with the paper's approach.
  3. For $1\leq i\neq j\leq 3$, compute covariance $c_{ij} = p_{ij} - p_i \cdot p_j$. This gives us the following covariance matrix $$ \begin{bmatrix} 1.00 & -0.03 & -0.01 \\ -0.03 & 1.00 & 0.03 \\ -0.01 & 0.03 & 1.00 \end{bmatrix} $$ However, the paper reports a different covariance matrix: $$ \begin{bmatrix} 1.0000 & -0.4464 & -0.1196 \\ -0.4464 & 1.0000 & 0.4442 \\ -0.1196 & 0.4442 & 1.0000 \end{bmatrix} $$

I would like to understand what would be the correct approach.

$\endgroup$
7
  • 1
    $\begingroup$ I’m not sure this is on-topic here, as it looks like code debugging or an inquiry about a particular function to use. The statistical theory to get dependent data, though, would be a “copula” such as the BBC copula package in R that might have something analogous in Python (perhaps in statsmodels). A hack might be to generate correlated $N(0,1)$ Gaussians in numpy and the write a list comprehension like [0 if < 0 else 1 for item in margin_x] (that syntax needs work but is just supposed to give the idea of what to do). $//$ Focusing on copulas would be totally on-topic here, however! $\endgroup$
    – Dave
    Commented Apr 4 at 13:38
  • $\begingroup$ @Dave I appreciate the gaussian copula suggestion. My question is about the algorithm described in the paper, and not a specific package. I have edited the question for clarification. $\endgroup$
    – zigs211567
    Commented Apr 4 at 15:05
  • 1
    $\begingroup$ There's a basic problem here: the pairwise correlations do not uniquely determine the joint distribution. Maybe you have a particular kind of joint Bernoulli distribution in mind? What is it? BTW, your post still emphasizes the issue of figuring out what's wrong with your code, which is a debugging task. It would be better to describe what you think it's doing. $\endgroup$
    – whuber
    Commented Apr 4 at 15:15
  • 1
    $\begingroup$ Thank you and I appreciate your patience. I have edited the question, removed the code, and provided the steps of my approach. $\endgroup$
    – zigs211567
    Commented Apr 4 at 15:58
  • 1
    $\begingroup$ @whuber, the problem seems well enough defined to me. The post asks for 8 probabilities satisfying 7 equations: summing to unity, and having the three given variances and three given correlations. This indeed leaves 1 degree of freedom, but inequalities for non-negative probabilities will constrain possible answers to a narrow range. So a nice algorithm that checks whether answers exist, and provides an interior answer if they do (perhaps in combination with two extreme answers) would be a feasible and worthwhile answer to the post. $\endgroup$
    – user225256
    Commented Apr 6 at 11:34

1 Answer 1

1
$\begingroup$

I have figured out how the algorithm works in the paper. I'm posting this answer just in case someone else is looking for the same problem.

First, note that as people have mentioned, this problem is underspecified for $M>2$. There are $2^M$ possible masses for the joint distribution of $M$ Bernoullis and there are $1 + M +\binom{M}{2}$ constraints. (The 1 is from the sum of the probability masses to be 1, $M$ is from $M$ variances, and $\binom{M}{2}$ are from the pairwise correlations.)

However, sampling from such distributions is often necessary, for instance in simulating clinical trials, the default risk of corporate bonds, or slashing risk in L2 blockchains.

The following algorithm works for practical purposes but note that the distribution it finds is not unique and it will fail silently if there is no such distribution. Hence, always verify the sample correlations.

Intuition

Let's start with the $M=2$ case. Suppose we want to sample from the joint distribution of two Bernoulli distributions with parameters $p_1$ and $p_2$ and their (Pearson) correlation $\rho$.

The first insight is observing that we can sample efficiently from a correlated bivariate Normal distribution and then convert the samples to Bernoulli like below:

converting Normal to Bernoulli

The normal samples on the first quadrant become $(1, 1)$, the samples on the 2nd quadrant become $(0, 1)$, and so on. The resulting Bernoulli samples are indeed correlated, but how can we make sure the correlation coefficient is $\rho$? We need to select the right normal distribution for that. Hence, to obtain a target correlation between the Bernoulli distributions, we find a normal distribution that satisfies our requirements.

Now, I'll introduce some math. Suppose $-1<\alpha<1$ and assume that the bivariate normal distribution $\mathcal N(\mu, \Sigma(\alpha))$ with mean $[\mu_x, \mu_y]$ and covariance matrix $\Sigma(\alpha)=[[1, \alpha],[\alpha, 1]]$ satisfies our requirements.

Suppose $(X_n, Y_n)\sim\mathcal N(\mu, \Sigma(\alpha))$ and $(X_b, Y_b)$ is the corresponding Bernoulli sample.

We want $$P(X_b=1) = p_1\implies P(X_n>0)=p_1\implies P(X_n-\mu_x>-\mu_x)=p_1$$ Now, $X_n-\mu_x$ is a standard normal since the variance of $X_n$ was 1. Thus, using the quantile function, $Q(\cdot)$, we have $\mu_x=Q(p_1)$.

Hence, $\mu=[Q(p_1), Q(p_2)]$.

Next step is to find $\Sigma(\alpha)$. Note that from the Pearson correlation formula

$$ \mathrm{Corr}(X_b, Y_b)=\rho\implies P(X_b=1, Y_b=1) = p_1p_2 + \rho\sqrt{p_1(1-p_1)p_2(1-p_2)} \quad(1) $$

On the other hand

$$ P(X_b=1, Y_b=1) = P(X_n > 0, Y_n>0) = \int_{0}^{\infty}\int_{0}^{\infty}f(x,y, \mu, \Sigma(\alpha))dydx \quad(2) $$

where $f(x,y, \mu, \Sigma(\alpha))$ is the pdf of the $\mathcal N(\mu, \Sigma(\alpha))$. For a given $\alpha$, this can be solved numerically (using scipy.stats for instance). We want an $\alpha$ so that the values in (1) and (2) are equal. Since $\alpha$ has a finite range, we can just search for such an $\alpha$.

The Algorithm $M=2$ Case

  1. Given $p_1$, $p_2$ and $\rho$, compute $P(X_b=1, Y_b=1)$ using (1).
  2. Compute $\mu_x=Q(p_1)$ and $\mu_y=Q(p_2)$, where $Q$ is the quantile function. This gives us the mean $\mu=[\mu_x, \mu_y]$ of the normal distribution.
  3. Now, for $-1<\alpha<1$s, define $\Sigma(\alpha)=[[1, \alpha],[\alpha, 1]]$. For various $\alpha$-s, compute $P(X_n > 0, Y_n>0)$ using (2). We want to find an $\alpha$ so that $P(X_n > 0, Y_n>0)$ matches with the step 1 value as closely as possible. A binary search on the range of $\alpha$ can be beneficial here. Suppose the $\tilde{\alpha}$ is the result of this step.
  4. Thus we know both $\mu$ and $\Sigma(\tilde{\alpha})$ of the normal distribution. Generate samples $(X_n, Y_n)$ from the distribution.
  5. Define $X_b=(X_n > 0)$ and similarly $Y_b = (Y_n > 0)$ (using numpy notations). These are the correlated Bernoulli samples.

$M > 2$ Case

In this case, the Normal distribution will be $M$-dimensional. The $\mu$ can be found using step 2. The diagonal entries of the covariance matrix will be $1$. And for entry $e_{ij}$ ($i$-th row and $j$-th column with $i\neq j$), we will use the algorithm for $M=2$ case for $p_i$, $p_j$ and $\rho_{ij}$ to find $e_{ij}=\tilde{\alpha_{ij}}$.

Finally, once we have the full covariance matrix, we will sample from the normal distribution, and using step 5, convert them to Bernoulli.

Conclusion

I have tested this algorithm for large $M$, and a wide range of $p_i$s and it seems to work very well unless either $p_i$ or $\rho_{ij}$ are near their boundaries (like $p_i\to 0$ or $\rho_{ij}\to1$ etc).

$\endgroup$
1
  • $\begingroup$ +1, especially for the illustration, which clearly shows what's going on. $\endgroup$
    – whuber
    Commented Apr 12 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.