---
abstract: 'Consider the problem of drawing random variates $(X_1,\ldots,X_n)$ from a distribution where the marginal of each $X_i$ is specified, as well as the correlation between every pair $X_i$ and $X_j$. For given marginals, the Fréchet-Hoeffding bounds put a lower and upper bound on the correlation between $X_i$ and $X_j$. Hence any achievable correlation can be uniquely represented by a convexity parameter $\lambda_{ij} \in [0,1]$ where 1 gives the maximum correlation and 0 the minimum correlation. We show that for a given convexity parameter matrix, the worst case is when the marginal distribution are all Bernoulli random variables with parameter 1/2 (fair 0-1 coins). It is worst case in the sense that given a convexity parameter matrix that is obtainable when the marginals are all fair 0-1 coins, it is possible to simulate from any marginals with the same convexity parameter matrix. In addition, we characterize completely the set of convexity parameter matrices for symmetric Bernoulli marginals in two, three and four dimensions.'
author:
- 'Mark Huber[^1]  and Nevena Marić[^2]'
title: Simulation of multivariate distributions with fixed marginals and correlations
---

Introduction
============

Consider the problem of simulating a random vector $(X_1,\ldots,X_n)$ where for all $i$ the cumulative distribution function (cdf) of $X_i$ is $F_i$, and for all $i$ and $j$ the correlation between $X_i$ and $X_j$ should be $\rho_{ij} \in [0,1]$. The correlation here is the usual notion $${\text{Corr}}(X,Y) = \frac{{\mathbb{E}}[(X - {\mathbb{E}}(X))(Y - {\mathbb{E}}(Y))]}{{\text{SD}}(X){\text{SD}}(Y)}
 = \frac{{\mathbb{E}}[XY] - {\mathbb{E}}[X]{\mathbb{E}}[Y]}{{\text{SD}}(X){\text{SD}}(Y)},$$ for standard deviations ${\text{SD}}(X)$ and ${\text{SD}}(Y)$ that are finite.

Let $\Omega$ denote the set of matrices with entries in $[-1,1]$, all the diagonal entries equal $1$, and are nonnegative definite. Then it is well known that any correlation matrix $(\rho_{ij})$ must lie in $\Omega$.

This problem, in different guises, appears in numerous fields: physics [@smith1981gust], engineering [@lampard1968stochastic], ecology [@Dias2008], and finance [@Lawrance1981], to name just a few. Due to its applicability in the generation of synthetic optimization problems, it has also received special attention by the simulation community  [@hill1994], [@henderson2000].

It is a very well studied problem with a variety of approaches. When the marginals are normal and the distribution is continuous with respect to Lebesgue measure, this is just the problem of generating a multivariate normal with specified correlation matrix. It is well known how to accomplish this (see, for instance [@fishman1996], p. 223) for any matrix in $\Omega$.

For marginals that are not normal, the question is very much harder. A common method is to employ families of copulas (see for instance [@nelson1999]), but there are very few techniques that apply to general marginals. Instead, different families of copulas typically focus on different marginal distributions.

Devroye and Letac [@devroyel2010] showed that if the marginals are beta distributed with equal parameters at least $1/2$, then when the dimension is three it is possible to simulate such a vector where the correlation is any matrix in $\Omega$. This set of beta distributions includes the important case of uniform $[0,1]$ marginals, but they have not been able to extend their technique to higher dimensions.

Chagnanty and Joe [@chagnantyj2006] characterized the achievable correlation matrices when the marginals are Bernoulli. When the dimension is 3 their characterization is easily checkable, in higher dimensions they give a number of inequalities that grows exponentially in the dimension.

For the case of general marginals, in statistics there is a tradition of using transformations of mutivariate normal vectors dating back to Mardia [@mardia1970] and Li and Hammond [@li1975]. This approach relies heavily on developing usable numerical methods. In this paper we approach the same problem using exclusively probabilistic techniques.

We show that for many correlation matrices the problem of simulating from a multivariate distribution with fixed marginals and specified correlation can be reduced to showing the existence of a multivariate distribution whose marginals are Bernoulli with mean $1/2$, and for each pair of marginals, there is a specified probability that the pair takes on the same value. For $n=2,3,4$ we are able to give necessary and sufficient conditions on those agreement probabilities in order for such a distribution to exist.

#### The convexity matrix approach.

Any two random variables $X$ and $Y$ have correlation in $[-1,1]$, but if the marginal distributions of $X$ and $Y$ are fixed, it is generally not possible to build a bivariate distribution for any correlation in $[-1,1]$. For instance, for $X$ and $Y$ both exponentially distributed, the correlation must lie in $[1 - \pi^2/6,1]$. The range of achievable correlations is always a closed interval.

For two dimensions it is well known how to find the minimum and maximum correlation. These come from the inverse transform method, which works as follows. First, given a cdf $F$, define the pseudoinverse of the cdf as $$F^{-1} = \inf\{x:F(x) \geq u\}.$$ When $U$ is uniform over the interval $[0,1]$ (write $U \sim {\textsf{Unif}}([0,1])$), $F^{-1}(U)$ is a random variable with cdf $F$ (see for instance p. 28 of [@devroye1986]). Since $U$ and $1 - U$ have the same distribution, both can be used in the inverse transform method. The random variables $U$ and $1 - U$ are [*antithetic*]{} random variables. Of course ${\text{Corr}}(U,U) = 1$ and ${\text{Corr}}(U,1-U) = -1$, so these represent an easy way to get minimum and maximum correlation when the marginals are uniform random variables.

The following theorem comes from work of Fréchet [@frechet1951] and Hoeffding [@hoeffding1940].

\[THM:fhbound\] For $X_1$ with cdf $F_1$ and $X_2$ with cdf $F_2$, and $U \sim {\textsf{Unif}}([0,1])$: $${\text{Corr}}(F_1^{-1}(U),F_2^{-1}(1 - U)) \leq {\text{Corr}}(X_1,X_2) \leq 
  {\text{Corr}}(F_1^{-1}(U),F_2^{-1}(U)).$$

In other words, the maximum correlation between $X_1$ and $X_2$ is achieved when the same uniform is used in the inverse transform method to generate both. The minimum correlation between $X_1$ and $X_2$ is achieved when antithetic random variates are used in the inverse transform method.

Now consider $(X_1,\ldots,X_n)$, where each $X_i$ has cdf $F_i$, and the correlation between $X_i$ and $X_j$ is $\rho_{ij}$. Then let $$\rho^-_{ij} = {\text{Corr}}(F_i^{-1}(U)F_j^{-1}(1 - U)) \text{ and }
\rho^+_{ij} = {\text{Corr}}(F_i^{-1}(U)F_j^{-1}(U)).$$ For the $\rho_{ij}$ to be achievable, it must be true that $\rho_{ij} \in [\rho_{ij}^-,\rho_{ij}^+]$, so there exists $\lambda_{ij} \in [0,1]$ so that $$\rho_{ij} = \lambda_{ij} \rho_{ij}^+ + (1 - \lambda_{ij}) \rho_{ij}^-.$$ Given a correlation matrix $(\rho_{ij})$ we will refer to $(\lambda_{ij})$ as the convexity parameter matrix. Our main result is as follows.

\[THM:main\] Suppose that for a convexity parameter matrix $\Lambda$, it is possible to simulate $(B_1,\ldots,B_n)$ where for all $i$, ${\mathbb{P}}(B_i = 1) = {\mathbb{P}}(B_i = 0) = 1/2$ and $$(\forall i,j)({\mathbb{P}}(B_i = B_j) = \lambda_{ij}).$$ Then for any specified set of marginals $F_1,\ldots,F_n$ it is possible to simulate in linear time from a multivariate distribution on $(X_1,\ldots,X_n)$ such that $$(\forall i)(X_i \sim F_i) \text{ and }
(\forall i,j)({\text{Corr}}(X_i,X_j) = 
  \lambda_{ij} \rho^+_{ij} + (1 - \lambda_{ij}) \rho^-_{ij}).$$

The next result characterizes when such a multivariate Bernoulli exists in two, three, and four dimensions, and gives necessary conditions for higher dimensions.

\[THM:nec\] Suppose $(B_1,B_2,\ldots,B_n)$ are random variables with ${\mathbb{P}}(B_i = 1) = {\mathbb{P}}(B_i = 0) = 1/2$ for all $i$. When $n = 2$, it is possible to simulate $(B_1,B_2)$ for any $\lambda_{12} \in [0,1]$. When $n = 3$, it is possible to simulate $(B_1,B_2,B_3)$ if and only if $$1 + 2\min\{\lambda_{23},\lambda_{12},\lambda_{13}\}
 \geq \lambda_{23} + \lambda_{12} + \lambda_{13} \geq 1.$$ When $n = 4$, it is possible to simulate $(B_1,B_2,B_3,B_4)$ if and only if $$1 - (1/2)\ell \leq (1/2)(u - 1)$$ where $$\begin{aligned}
\ell &= \min (\lambda_{14} + \lambda_{14} + \lambda_{13} + \lambda_{23},
 \lambda_{14} + \lambda_{34} + \lambda_{12} + \lambda_{23},
 \lambda_{24} + \lambda_{34} + \lambda_{12} + \lambda_{13}) \\
u &= \min_{\{i,j,k\}} \lambda_{ij} + \lambda_{jk} + \lambda_{ik}.\end{aligned}$$

The rest of the paper is organized as follows. In the next section the notion of a concurrence matrix is introduced, and it is shown that when the marginals are all Bernoulli with mean $1/2$, the concurrence matrix and convexity matrix are the same. Furthermore, given a multivariate Bernoulli $1/2$ with given convexity matrix, a straightforward algorithm gives the multivariate distribution with arbitrary marginals that has the same convexity matrix, proving Theorem \[THM:main\]. In Section \[SEC:nec\], the $n = 2$ and $n = 3$ cases are shown. In Section \[SEC:asymmetric\], a link is made between symmetric Bernoulli marginals and asymmetric, which in turn gives the proof of the $n = 4$ case.

The algorithm
=============

To present our algorithm, we begin by defining a matrix that measures the chance that any two components are equal.

For a random vector $(X_1,\ldots,X_n)$ the [*concurrence matrix*]{} $A = (a_{ij})$ is defined as $a_{ij} = {\mathbb{P}}(X_i = X_j)$.

Let ${\textsf{Bern}}(p)$ denote the Bernoulli distribution with parameter $p$, so for $X \sim {\textsf{Bern}}(p)$, ${\mathbb{P}}(X = 1) = p$ and ${\mathbb{P}}(X = 0) = 1-p$.

\[THM:concurrence\] For a random vector $(B_1,\ldots,B_n)$ with $B_i \sim {\textsf{Bern}}(1/2)$ for all $i$, the concurrence matrix and convexity matrix are the same matrix.

Let $(B_1,\ldots,B_n)$ be a random vector with marginals ${\textsf{Bern}}(1/2)$, and $i \neq j$ be elements of $\{1,\ldots,n\}$.

Let ${\textbf{1}}(\text{expression})$ denote the indicator function that is 1 if the expression is true and 0 otherwise. The inverse transform method for generating ${\textsf{Bern}}(1/2)$ random variables is $B = {\textbf{1}}(U > 1/2)$ where $U \sim {\textsf{Unif}}([0,1])$. From the Fréchet-Hoeffding bound (Theorem \[THM:fhbound\]), the correlation between $B_i$ and $B_j$ can range anywhere from -1 to 1. Hence $${\text{Corr}}(B_i,B_j) = \lambda_{ij} (1) + (1 - \lambda_{ij})(-1) \Rightarrow
 \lambda_{ij} = (1/2)({\text{Corr}}(B_i,B_j) + 1).$$

Now denote ${\mathbb{P}}(B_i = a,B_j = b)$ by $p_{ab}$. Then $$\begin{aligned}
{\mathbb{P}}(B_i = B_j) &= p_{11} + p_{00} = a_{ij} \\
{\mathbb{P}}(B_i = 1) &= p_{10} + p_{11} = 1/2 \\
{\mathbb{P}}(B_j = 1) &= p_{01} + p_{11} = 1/2 \\
\sum_{a,b} {\mathbb{P}}(B_i = a,B_j = b) &= p_{00} + p_{01} + p_{10} + p_{11} = 1\end{aligned}$$ Adding the top three equations and then using the fourth gives $$3p_{11} + p_{10} + p_{01} + p_{00} = a_{ij} + 1 = 2p_{11} + 1.$$ Hence $p_{11} = a_{ij} / 2.$

That means $${\text{Corr}}(B_i,B_j) = \frac{{\mathbb{E}}[B_i B_j] - {\mathbb{E}}[B_i]{\mathbb{E}}[B_j]}{{\text{SD}}(B_1){\text{SD}}(B_2)} 
  = \frac{(1/2)a_{ij} - (1/2)(1/2)}{(1/2)(1/2)} = 2a_{ij} - 1,$$ so $a_{ij} = (1/2)({\text{Corr}}(B_i,B_j) + 1)$, and $a_{ij} = \lambda_{ij}$.

Given the ability to generate from a multivariate Bernoulli with given convexity/concurrence matrix, generating from a multivariate with general marginals with the same convexity matrix turns out to be easy. Suppose the $n$ by $n$ convexity matrix $\Lambda=(\lambda_{ij})$ and cdf’s $F_1,\ldots,F_n$ are given.

[rl]{}\
1) & Draw $U$ uniformly from $[0,1]$.\
2) & Draw $(B_1,\ldots,B_n)$ $(B_i \sim {\textsf{Bern}}(1/2)$) with concurrence matrix $\Lambda = (\lambda_{ij})$.\
3) & For all $i \in \{1,\ldots,n\}$, let $X_i \leftarrow F_i^{-1}(UB_i + (1 - U)(1 - B_i)).$\

That is, if $B_i = 1$ then choose $X_i$ with the inverse transform method and $U$, otherwise use the inverse transform method and $1 - U$. Theorem \[THM:main\] can now be restated in terms of this algorithm.

\[THM:convexity\] The above algorithm generates $(X_1,\ldots,X_n)$ with $X_i \sim F_i$ for all $i$ and convexity matrix $\Lambda$.

Each $X_i$ can be written as $B_i F_i^{-1}(U) + (1 - B_i) F_i^{-1}(1 - U)$. That is, it is a mixture of $F_i^{-1}(U)$ and $F_i^{-1}(1 - U)$ with weights equal to $1/2$. Since $F_i^{-1}(U)$ and $F_i^{-1}(1 - U)$ have the same distribution with cdf $F_i$, the random variable $X_i$ also shares this distribution. Hence they all have the same mean (call it $\mu_i$) and standard deviation (call it $\sigma_i$.)

That means the correlations satisfy $$\begin{aligned}
{\text{Corr}}(X_i,X_j) \sigma_i \sigma_j &= {\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)] \\
&= {\mathbb{E}}[{\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)]|B_i,B_j] \\
&= \sum_{{a \in \{0,1\}} \atop {b \in \{0,1\}}} {\mathbb{P}}(B_i = a,B_j = b)
  {\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)|B_i = a,B_j = b].\end{aligned}$$

Note when $B_i = B_j = 0$, then $X_i = F_i^{-1}(U)$ and $X_j = F_j^{-1}(U)$ which are maximally correlated. So $${\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)|B_i = 0,B_j = 0]
 = \rho^+_{ij} \sigma_i \sigma_j.$$ Following this same logic further gives $$\begin{aligned}
{\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)|B_i = 1,B_j = 0]
 &= \rho^-_{ij} \sigma_i \sigma_j \\
{\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)|B_i = 0,B_j = 1]
 &= \rho^-_{ij} \sigma_i \sigma_j \\
{\mathbb{E}}[(X_i - \mu_i)(X_j - \mu_j)|B_i = 1,B_j = 1]
 &= \rho^+_{ij} \sigma_i \sigma_j \\\end{aligned}$$ Hence $${\text{Corr}}(X_i,X_j) = \rho^+_{ij} {\mathbb{P}}(B_i = B_j) + \rho^-_{ij} {\mathbb{P}}(B_i \neq B_j)
  = \rho^+_{ij} a_{ij} + \rho^-_{ij}(1 - a_{ij}).$$

Theorem \[THM:main\] is an immediate consequence.\

The complete algorithm for generating multivariate distributions where component $i$ has cdf $F_i$ and ${\text{Corr}}(X_i,X_j) = \rho_{ij}$ is as follows.

[rl]{}\
1) & For every $\{i,j\} \subseteq \{1,\ldots,n\}$ do (where $A \sim {\textsf{Unif}}([0,1])$)\
2) & $\rho^+_{ij} \leftarrow {\text{Corr}}(F_i^{-1}(A),F_i^{-1}(1 - A)),$ $\rho^+_{ij} \leftarrow {\text{Corr}}(F_i^{-1}(A),F_i^{-1}(1 - A))$.\
3) & $\lambda_{ij} \leftarrow (\rho_{ij} - \rho^-_{ij})/(\rho^+_{ij} - \rho^-_{ij})$\
4) & Draw $U \leftarrow {\textsf{Unif}}([0,1])$\
5) & Draw $(B_1,\ldots,B_n)$ $(B_i \sim {\textsf{Bern}}(1/2)$) with concurrence matrix $\Lambda = (\lambda_{ij})$.\
6) & For all $i \in \{1,\ldots,n\}$, let $X_i \leftarrow F_i^{-1}(U B_i + (1 - U)(1 - B_i)).$\

Line 2 just generates the maximum and minimum correlation bounds from Theorem \[THM:fhbound\]. Line 3 is set up so that $\rho_{ij} = \lambda_{ij} \rho^+_{ij} + (1 - \lambda_{ij})\rho^-_{ij}.$

Note that when the marginals for $i$ and $j$ are the same distribution, $\rho^+_{ij} = 1$. However, even for equal marginals, the minimum correlation is often strictly greater than $-1$. See [@dukicm2013] and the references therein for the minimum correlation for a variety of common distributions.

This algorithm is only applicable when it is actually possible to simulate from the multivariate symmetric Bernoulli distribution with specified concurrence matrix. The next section gives necessary and sufficient conditions for this to be possible for two, three, and four dimensions.

Attainable convexity matrices {#SEC:nec}
=============================

In this section we prove Theorem \[THM:nec\] which gives necessary and sufficient conditions for dimensions two, three, and four for when it is possible to create a multivariate Bernoulli 1/2 with a given convexity matrix. We break the Theorem into three parts.

\[LEM:two\] For any $\lambda_{12} \in [0,1]$, there exists a unique joint distribution on $\{0,1\}^2$ such that $(B_1,B_2)$ with this distribution has $B_1,B_2 \sim {\textsf{Bern}}(1/2)$ and ${\mathbb{P}}(B_1 = B_2) = \lambda_{12}$.

Let $p_{ij} = {\mathbb{P}}(B_1 = i,B_2 = j)$. Then the equations that are necessary and sufficient to meet the distribution and convexity conditions are: $$p_{10} + p_{11} = 0.5,\ p_{01} + p_{11} = 0.5,\ p_{11} + p_{00} = \lambda_{12},
 \text{ and } p_{00} + p_{01} + p_{10} + p_{11} = 1.$$ This system of linear equations has full rank, so there exists a unique solution. Given there exists a unique solution, it is easy to verify that $$p_{00} = (1/2)\lambda_{12},\ p_{01} = (1/2)[1-\lambda_{12}],
  \ p_{10} = (1/2)[1-\lambda_{12}],\ p_{11} = (1/2)\lambda_{12}$$ satisfies the equations.

This provides an alternate algorithm to that found in [@dukicm2013] for simulating from bivariate distributions with correlation between $\rho^-_{1,2}$ and $\rho^+_{1,2}$.

\[LEM:3B\] A random vector $(B_1,B_2,B_3)$ with $B_i \sim {\textsf{Bern}}(1/2)$ exists (and is possible to simulate from in a constant number of steps) if and only if the concurrence matrix satisfies $$1 \leq \lambda_{23} + \lambda_{12} + \lambda_{13} \leq 
  1 + 2 \min\{\lambda_{12},\lambda_{13},\lambda_{23}\}$$

Letting $p_{ijk} = {\mathbb{P}}(B_1 = i,B_2 = j,B_3 = k)$, there are three conditions from the marginals. $$\sum_{j,k \in \{0,1\}} p_{1jk} = 0.5,\ 
\sum_{i,k \in \{0,1\}} p_{i1k} = 0.5,\ 
\sum_{ij \in \{0,1\}} p_{ij1} = 0.5.$$ Then there are three conditions from the correlations $$\sum_{k \in \{0,1\}} p_{00k} + p_{11k} = \lambda_{12},\ 
\sum_{j \in \{0,1\}} p_{0j0} + p_{1j1} = \lambda_{13},\ 
\sum_{i \in \{0,1\}} p_{i00} + p_{i11} = \lambda_{23}.$$ A seventh condition is $\sum_{i,j,k} p_{ijk} = 1$. Since eight equations are needed, suppose that $p_{111} = \alpha$.

This 8 by 8 system of equations has full rank, so there is a unique solution. It is easy to verify that the solution is $$\begin{aligned}
p_{000} &= (1/2)(\lambda_{12} + \lambda_{13} + \lambda_{23} - 1) - \alpha \\
p_{001} &= (1/2)(1-(\lambda_{13} + \lambda_{23})) + \alpha \\
p_{010} &= (1/2)(1-(\lambda_{12} + \lambda_{23})) + \alpha \\
p_{011} &= (1/2)\lambda_{23} - \alpha \\
p_{100} &= (1/2)(1-(\lambda_{12} + \lambda_{13})) + \alpha \\
p_{101} &= (1/2)\lambda_{13} - \alpha \\
p_{110} &= (1/2)\lambda_{12} - \alpha \\
p_{111} &= \alpha \end{aligned}$$

In order for this solution to yield probabilities, all must lie in $[0,1]$. The $p_{011},p_{101},$ and $p_{110}$ equations give the restriction $$\label{EQN:restriction1}
\alpha \leq (1/2)\min\{\lambda_{12},\lambda_{23},\lambda_{13}\}.$$ Of course, $\alpha \geq 0$ so the $p_{111}$ equation can be satisfied.

The $p_{000}$ equation requires that $$\label{EQN:restriction2}
\alpha \leq (1/2)(\lambda_{13} + \lambda_{12} + \lambda_{23}  - 1).$$ With these two conditions, equations $p_{001}$, $p_{010}$, and $p_{100}$ give the constraint $$\label{EQN:restriction3}
\alpha \geq (1/2)(\lambda_{13} + \lambda_{12} + \lambda_{23} - 
  \min\{\lambda_{13},\lambda_{12},\lambda_{23}\} - 1).$$

The lower bound in  is below the upper bound in equation , but not necessarily below the upper bound in equation . So this gives another restriction: $$\label{EQN:restriction4}
(1/2)(\lambda_{13} + \lambda_{12} + \lambda_{23} - 
  \min\{\lambda_{13},\lambda_{12},\lambda_{23}\} - 1) \leq 
   (1/2) \min\{\lambda_{13},\lambda_{12},\lambda_{23}\}.$$

As long as  and  are satisfied, their will exist a solution. If satisfied with strict inequality, there will be an infinite number of solutions.

Note that not all positive definite correlation matrices are attainable with ${\textsf{Bern}}(1/2)$ marginals. For instance, if $\lambda_{12} = \lambda_{13} = \lambda_{23} = 0.3$, then $\rho_{12} = \rho_{13} = \rho_{23} = -0.4.$ With diagonal entries $1$, the $\rho$ values form a positive definite matrix, but it is impossible to build a multivariate distribution with ${\textsf{Bern}}(1/2)$ marginals with these correlations.

Asymmetric Bernoulli distributions {#SEC:asymmetric}
==================================

Now consider the problem of drawing a multivariate Bernoulli $(X_1,\ldots,X_n)$ where $X_i \sim {\textsf{Bern}}(p_i)$ where $i$ is not necessarily $1/2$, and the concurrence matrix is given.

\[LEM:asymmetric\] An $n$ dimensional multivariate Bernoulli distribution where the marginal of component $i$ is ${\textsf{Bern}}(p_i)$ and concurrence matrix $\Lambda$ exists if and only if an $n + 1$ dimensional multivariate Bernoulli distribution exists with ${\textsf{Bern}}(1/2)$ marginals and concurrence matrix $$\left(
\begin{array}{c|c}
\raisebox{-15pt}{{\Large\mbox{{$\Lambda$}}}} & p_1 \\[-4ex]
 & \vdots \\
 & p_n \\
\hline
p_1 \cdots p_n & 1 \\ 
\end{array}
\right)$$

Suppose such an $n + 1$ dimensional distribution exists with ${\textsf{Bern}}(1/2)$ marginals and specified concurrence matrix. Let $(B_1,\ldots,B_{n+1})$ be a draw from this distribution. Then set $X_i = {\textbf{1}}(B_{i} = B_{n+1}).$ The concurrence matrix gives ${\mathbb{P}}(X_i = 1) = p_i$, and for $i \neq j$, ${\mathbb{P}}(X_i = X_j) = {\mathbb{P}}(B_{i} = B_{j}) = \lambda_{ij}$.

Conversely, suppose such an $n$ dimensional distribution with ${\textsf{Bern}}(p_i)$ marginals exists. Let $B_{n+1} \sim {\textsf{Bern}}(1/2)$ independent of the $X_i$, and set $B_{i} = B_{n+1} X_i + (1 - B_{n+1})(1 - X_i).$ Then ${\mathbb{P}}(B_{i} = 1) = (1/2)p_i + (1/2)(1 - p_i) = 1/2$, and ${\mathbb{P}}(B_{i} = B_{n+1}) = p_i$, the correct concurrence parameter. Finally, for $i \neq j$, $${\mathbb{P}}(B_{i} = B_{j}) = {\mathbb{P}}(X_i = X_j) = \lambda_{ij}.$$

This result gives an alternate way of deriving Lemma \[LEM:3B\] using the Fréchet-Hoeffding bounds. Consider a bivariate asymmetric Bernoulli $(X,Y)$ where $X \sim {\textsf{Bern}}(p)$, and $Y \sim {\textsf{Bern}}(q)$. The pseudoinverse cdf of (p) is given by $F_p^{-1}(U) = {\textbf{1}}( U > 1-p)$. Hence the minimum and maximum correlation of $X$ and $Y$ can be directly found using Theorem \[THM:fhbound\]: $$\begin{aligned}
 \rho^{-}(X,Y)&= \frac{(p+q-1) {\textbf{1}}(p+q >1) - pq}{\sqrt{pq(1-p)(1-q)}}\\
   \rho^{+}(X,Y)&= \frac{\min(p,q) - pq}{\sqrt{pq(1-p)(1-q)}}.
 \end{aligned}$$

Suppose that we want to simulate $X$ and $Y$ so that ${\mathbb{P}}(X = Y) = r$. For what $r$ is this possible?

Let $p_{ij} = {\mathbb{P}}(X = i,Y = j)$. Then $$\begin{aligned}
 {\mathbb{P}}(X=Y)&=r = p_{11}+p_{00}\\
 {\mathbb{P}}(X=1)&=p = p_{10}+p_{11}\\
 {\mathbb{P}}(Y=1)&=q = p_{01}+p_{11}
 \end{aligned}$$ Adding the three equations above we have $$\begin{aligned}
 p+q+r =  p_{11}+p_{00}+p_{10}+p_{11}+p_{01}+p_{11} = 1+ 2 p_{11}
 \end{aligned}$$ and therefore $p_{11}= (p+q+r -1)/2$. On the other hand $${\mathbb{E}}(XY) = {\mathbb{E}}({\textbf{1}}(X=1, Y=1)) = {\mathbb{P}}(X=1,Y=1)= p_{11}.$$ The Fréchet-Hoeffding Theorem provides the condition $${\mathbb{E}}(F_p^{-1}(U)F_q^{-1}(1 - U)) \leq {\mathbb{E}}(XY) \leq {\mathbb{E}}(F_p^{-1}(U)F_q^{-1}(U))$$ which is equivalent to $$(p+q-1) {\textbf{1}}(p+q >1) \leq (p+q+r -1)/2 \leq  \min(p,q).$$ Solving the inequalities for $r$ we get the following condition $$|1-(p+q)| \leq r \leq 2 \min(p,q) + 1- (p+q).$$

Letting $p=\lambda_{13}$, $q=\lambda_{23}$, and $r=\lambda_{12}$ we obtain an equivalent condition to the one in Lemma \[LEM:3B\].

As a further consequence of this new proof technique, we obtain a simple method for simulating from the 3 dimensional symmetric Bernoulli distribution for given concurrence probabilities when possible.

[rl]{}\
\
\
1) & Draw $B_3 \leftarrow {\textsf{Bern}}(1/2)$ and $U \leftarrow {\textsf{Unif}}([0,1])$.\
2) & $X_1 \leftarrow {\textbf{1}}(U \in [0,\lambda_{13}])$\
3) & $X_2 \leftarrow {\textbf{1}}(U \in 
  [(1 + \lambda_{13} - \lambda_{23} - \lambda_{12})/2,
   (1 + \lambda_{13} + \lambda_{23} - \lambda_{12})/2])$\
4) & $B_1 \leftarrow B_3 X_1 + (1 - B_3)(1 - X_1)$\
5) & $B_1 \leftarrow B_3 X_2 + (1 - B_3)(1 - X_2)$\

\[PROP:threevaralg\] The above algorithm simulates $(B_1,B_2,B_3)$ where $B_i \sim {\textsf{Bern}}(1/2)$ for all $i$ and ${\mathbb{P}}(B_i = B_j) = \lambda_{ij}$ for all $i$ and $j$.

Note $1 \leq \lambda_{12} + \lambda_{13} + \lambda_{23} \leq 
 1 + 2\min\{\lambda_{12},\lambda_{13},\lambda_{23}\}$ if and only if $$0 \leq \frac{1 + \lambda_{13} - \lambda_{23} - \lambda_{12}}{2} \leq \lambda_{13}
 \leq \frac{1 + \lambda_{13} + \lambda_{23} - \lambda_{12}}{2} \leq 1.$$ Hence $$\begin{aligned}
{\mathbb{P}}(X_1 = 1) &= \lambda_{13} - 0 = \lambda_{13}. \\
{\mathbb{P}}(X_2 = 1) &= \frac{1 + \lambda_{13} + \lambda_{23} - \lambda_{12}}{2}
 - \frac{1 + \lambda_{13} - \lambda_{23} - \lambda_{12}}{2} = \lambda_{23}. \\
{\mathbb{P}}(X_1 = X_2) &= \left(\lambda_{13} -  
  \frac{1 + \lambda_{13} - \lambda_{23} - \lambda_{12}}{2}\right) + 
  \left(1 - \frac{1 + \lambda_{13} + \lambda_{23} - \lambda_{12}}{2}\right) 
  = \lambda_{12}.\end{aligned}$$ Then as in the proof of Lemma \[LEM:asymmetric\], the output has the desired trivariate symmetric Bernoulli distribution.

Lemma \[LEM:asymmetric\] also allows a characterization of the 4 dimensional symmetric Bernoulli concurrence (and so also convexity) matrices.

A random vector $(B_1,B_2,B_3,B_4)$ with $B_i \sim {\textsf{Bern}}(1/2)$ exists (and is possible to simulate in a constant number of steps) if and only if for $$\begin{aligned}
\ell &= \min (\lambda_{14} + \lambda_{14} + \lambda_{13} + \lambda_{23},
 \lambda_{14} + \lambda_{34} + \lambda_{12} + \lambda_{23},
 \lambda_{24} + \lambda_{34} + \lambda_{12} + \lambda_{13}) \\
u &= \min_{\{i,j,k\}} \lambda_{ij} + \lambda_{jk} + \lambda_{ik},\end{aligned}$$ it is true that $$1 - (1/2)\ell \leq (1/2)(u - 1)$$

By using Lemma \[LEM:asymmetric\], the problem is reduced to finding a distribution for $(X_1,X_2,X_3)$ where $X_i \sim {\textsf{Bern}}(\lambda_{i4})$ and the upper 3 by 3 minor of $\Lambda$ is the new concurrence matrix. Just as in Lemma \[LEM:3B\], this gives eight equations of full rank with a single parameter $\alpha$. Letting $q_{ijk} = {\mathbb{P}}(X_1 = i,X_2 = j,X_3 = k)$, the unique solution is $$\begin{aligned}
q_{000} &= (1/2)(\lambda_{12} + \lambda_{13} + \lambda_{23}) - (1/2) - \alpha \\
q_{001} &= -(1/2)(\lambda_{14} + \lambda_{24} + \lambda_{13} + \lambda_{23}) + 1 
   + \alpha \\
q_{010} &= -(1/2)(\lambda_{14} + \lambda_{34} + \lambda_{12} + \lambda_{23}) + 1 
   + \alpha \\
q_{011} &= (1/2)(\lambda_{24} + \lambda_{34} + \lambda_{23}) - (1/2) - \alpha \\
q_{100} &= -(1/2)(\lambda_{24} + \lambda_{34} + \lambda_{12} + \lambda_{13}) + 1 
   + \alpha \\
q_{101} &= (1/2)(\lambda_{14} + \lambda_{34} + \lambda_{13}) - (1/2) - \alpha \\
q_{110} &= (1/2)(\lambda_{34} + \lambda_{24} + \lambda_{12}) - (1/2) - \alpha \\
q_{111} &= \alpha.\end{aligned}$$

All of these right hand sides lie in $[0,1]$ if and only if $1 - (1/2)\ell \leq (1/2)(u - 1)$, and $\alpha$ is chosen to lie in $[1 - (1/2)\ell,(1/2)(u - 1)]$.

As with the 3 dimensional case, this proof can be used to simulate a 4 dimensional multivariate symmetric Bernoulli: generate $(X_1,X_2,X_3)$ using any $\alpha \in [1 - (1/2)\ell,(1/2)(u - 1)]$ and the $q$ distribution, then generate $B_4 \sim {\textsf{Bern}}(1/2)$, and then set $B_i$ to be $B_4 X_i + (1 - B_4)(1 - X_i)$ for $i \in \{1,2,3,4\}$.

Conclusions
===========

The Fréchet-Hoeffding bounds give a lower and upper bound on the pairwise correlation between two random variables with given marginals. Hence for higher dimensions the correlation matrix can be written as a convexity matrix whose parameters indicated where on the line from the lower to the upper bound the correlation lies. With a simple algorithm, we have shown that when viewed as a convexity matrix problem, the worse case for marginals is when they are all symmetric Bernoulli random variables. It is worst case in the sense that if it is possible to build a symmetric Bernoulli multivariate distribution with the given convexity matrix, then it is possible to build a multivariate distribution with that convexity matrix [*for any marginal distributions.*]{} For symmetric Bernoulli marginals, the convexity matrix is also the concurrence matrix that gives the probabilities that any pair of random variables are the same.

For three and four dimensions, the set of convexity matrices that yield a symmetric Bernoulli distribution is characterized completely. For five or higher dimensions, every subset of three and four have these characterizations as necessary conditions.

[**Acknowledgement**]{} Support from the National Science Foundation (grant DMS - 1007823) is gratefully acknowledged.

[10]{}

N. R. Chagnanty and H. Joe. Range of correlation matrices for dependent Bernoulli random variables. , 93:197–206, 2006.

L. Devroye. . Springer, 1986.

L. Devroye and G. Letac. Copulas in three dimensions with prescribed correlations. 2010. arXiv:1004.3146v1.

C. T. S. Dias, A. Samaranayaka, and B. Manly. . , 215(4): 293–300, 2008.

V. M Dukic and N. Mari[ć]{}. Minimum correlation in construction of multivariate distributions. , 87(3):032114, 2013.

G. S. Fishman. . Springer-Verlag, 1996.

M. Freéchet. Sur les tableaux de corrélation dont les marges sont donnés. , 4, 1951.

S. G. Henderson, B. A. Chiera, and R. M. Cooke. Generating ÒdependentÓ quasi-random numbers. , Vol. 1. IEEE, 2000.

R. R. Hill and C. H. Reilly. Composition for multivariate random variables. , 332–339, IEEE, 1994.

W. Hoeffding. Masstabinvariante korrelatiostheorie. , 5:179–233, 1940.

D. G. Lampard. A stochastic process whose successive intervals between events form a first order Markov chain: I. , pages 648–668, 1968.

A. J. Lawrance and P. A. W. Lewis. . , 13(4):826–845, 1981.

S. T. Li and J. L. Hammond. Generation of pseudorandom numbers with specified univariate distributions and correlation coefficients. (1975): 557-561.

K. V. Mardia. A translation family of bivariate distributions and Fréchet’s bounds. , Series A (1970): 119-122.

R. B. Nelson. . Springer, 1999.

O. E. Smith and S. I. Adelfang. Gust model based on the bivariate gamma probability distribution. , 18(6):545–549, 1981.

[^1]: Claremont McKenna College (email: [mhuber@cmc.edu]{})

[^2]: University of Missouri-St. Louis (email: [maric@math.umsl.edu]{})
