Multivariate Bernoulli Distribution

• The multivariate Bernoulli distribution is a probability model for binary vectors defined with fixed margins and potential dependence structures via pairwise correlations.
• Its joint density can be represented as a convex combination of extremal ray densities, forming a convex polytope subject to margin constraints.
• Algorithmic construction through linear programming ensures feasible correlation assignments, making it essential for simulation in risk aggregation and discrete multivariate modeling.

A multivariate Bernoulli distribution is a probability law for a vector $X = (X_1,\dots,X_m)$ of binary random variables, $X_i \in \{0,1\}$, defined on the discrete hypercube $\{0,1\}^m$, with prescribed marginal distributions (each $X_i \sim \text{Bernoulli}(p_i)$) and, possibly, specified dependence structure such as pairwise correlation matrix $(\rho_{ij})$. The Fréchet class of such distributions comprises all joint laws with fixed marginals $p = (p_1,\dots,p_m)$. Characterizing and constructing multivariate Bernoulli laws with given margins and correlations is central to discrete multivariate modeling, combinatorics, and dependent risk simulation (Fontana et al., 2017).

1. Fréchet Class: Definition and Structure

For $m \geq 1$ and margin vector $p = (p_1, \dots, p_m)$ with $0 < p_i < 1$, the Fréchet class is

$$F(p_1,\dots,p_m) = \{\, F \in F_m: X_i \sim \text{Bernoulli}(p_i),\, i=1,\dots,m \,\}$$

where $F_m$ denotes all distribution functions on $\{0,1\}^m$. This class is equivalently described in terms of joint densities $f: \{0,1\}^m \to [0,1]$ with the constraints: $$\sum_{x \in \{0,1\}^m} f(x) = 1, \qquad \sum_{x: x_i=1} f(x) = p_i, \;\; i=1,\dots,m$$ Thus, $F(p)$ is a convex polytope cut out by the $m$ margin constraints and the simplex condition. This setup encompasses all feasible dependence structures (including both extreme positive/negative association and independence) compatible with the specified margins (Fontana et al., 2017).

2. Convex Geometry and Extremal Representation

Theorem 3.2 in (Fontana et al., 2017) formally establishes that the set $\mathcal{D}(p)$ of all densities $f$ with margins $p$ is a convex polytope in $\mathbb{R}^{2^m}$, whose extreme points (vertices) can be determined algebraically. Explicitly, every $f \in \mathcal{D}(p)$ is a convex combination

$f(x) = \sum_{i=1}^{n_F} \lambda_i R^{(i)}(x)$

where $\{R^{(i)}\}_{i=1}^{n_F}$ are the ray densities (vertices or extremal distributions) of the polytope, $\lambda_i \geq 0$, and $\sum_{i} \lambda_i = 1$. These ray densities are computed as nonnegative solutions to the homogeneous margin-constraint system $H r = 0$ (where $H$ is the $m \times 2^m$ margin-constraint matrix), and generally possess support of minimal size dictated by the polytope's combinatorics (Fontana et al., 2017).

3. Polynomial and Copula-Type Expansions

Every $f \in \mathcal{D}(p)$ admits a unique polynomial (Farlie-Gumbel-Morgenstern-type) expansion indexed by subsets $A \subseteq \{1, \dots, m\}$: $$f(x) = \left( D^{\otimes m} F \right)(x) = \sum_{A \subseteq \{1,\dots,m\}} \theta_A \prod_{i \in A} (-1)^{1-x_i} p_i^{x_i} q_i^{1-x_i}$$ where $q_i = 1 - p_i$, $D$ is the finite difference operator, and the parameter vector $(\theta_A)$ encodes dependence at all interaction orders up to $m$ (Fontana et al., 2017). The cumulative function $F$ is expressed as a multilinear polynomial in $(1 - u_i)$, generalizing copula representations in the discrete setting: $$F(u) = 1 + \sum_{i} \theta_i (1-u_i) + \sum_{i<j} \theta_{ij}(1-u_i)(1-u_j) + \ldots + \theta_{12\cdots m} \prod_{i=1}^m (1-u_i)$$ This structure mirrors the FGM copula, but, in finite support and with matching margins, the parameters must satisfy linear constraints induced by the margin equations.

4. Compatibility and Realizability of Correlation Matrices

For a given pairwise correlation matrix $(\rho_{ij})$, the compatibility problem is: does there exist $f \in \mathcal{D}(p)$ with pairwise second moments matching

$$\mathbb{E}[X_i X_j] = p_i p_j + \rho_{ij} \sqrt{p_i q_i p_j q_j}$$

for all $1 \leq i < j \leq m$? Proposition 3.1 in (Fontana et al., 2017) shows that this is the case if and only if the target vector of second moments $\mu_2$ lies in the convex hull of the second-moment columns of the ray matrix $A_{2,p}$. Explicitly, the feasibility system is

$$A_{2,p} \lambda = \mu_2, \qquad \sum_i \lambda_i = 1, \quad \lambda_i \geq 0$$

where $A_{2,p}\in\mathbb{R}^{m(m-1)/2 \times n_F}$ collects the pairwise moments of the rays. If a nonnegative $\lambda$ exists, the desired law can be constructed as $f = \sum_i \lambda_i R^{(i)}$. If not, $\mu_2$ is infeasible, and one may project onto the feasible region to find the nearest compatible correlation structure (Fontana et al., 2017).

5. Bounding Achievable Correlations

Proposition 3.2 gives explicit bounds for each pairwise correlation: $L_{ij} \leq \mathbb{E}[X_i X_j] \leq U_{ij}$ where $L_{ij}$ and $U_{ij}$ are the minimal and maximal pairwise moments among all rays. In terms of correlations,

$$\rho_{ij} \in \left[ \frac{L_{ij} - p_i p_j}{\sqrt{p_i q_i p_j q_j}},\, \frac{U_{ij} - p_i p_j}{\sqrt{p_i q_i p_j q_j}} \right]$$

The bivariate case reduces to the classical Fréchet–Hoeffding bounds: $$\max\{0, p_1 + p_2 - 1\} \leq \mathbb{E}[X_1 X_2] \leq \min\{p_1, p_2\}$$ This establishes the sharp compatibility region for each pair of margins and attainable correlation (Fontana et al., 2017).

6. Algorithmic Construction and Numerical Illustration

The procedure for constructing a multivariate Bernoulli law with prescribed margins and correlations involves:

1. Building the constraint matrix $H$ for the margin system $H f = 0$;
2. Computing the extremal rays $R^{(i)}$ as generators of the kernel cone (e.g., via 4ti2);
3. Assembling the ray matrix $R_p = [R^{(1)},\dots,R^{(n_F)}]$ and the second-moment matrix $A_{2,p}$;
4. Solving the linear program $$A_{2,p} \lambda = \mu_2, \sum \lambda = 1, \lambda \geq 0$$ for feasibility;
5. Forming the joint probability vector $f = R_p \lambda$ if a solution exists.

In the case $m=3$, $p_i = 0.5$ for $i = 1,2,3$, and desired $(\rho_{12},\rho_{13},\rho_{23}) = (0.2, -0.3, 0.4)$, the algorithm yields the explicit weight vector $\lambda$ and the corresponding joint pmf on $\{0,1\}^3$ (Fontana et al., 2017).

7. Theoretical and Practical Implications

This convex-geometric characterization unifies the study of discrete dependence, simulation, and parameter compatibility in multivariate Bernoulli laws. Through explicit polytope structure, any marginal/correlation assignment can be checked for feasibility; if feasible, extremal and mixed laws can be synthesized efficiently. The theoretical framework extends naturally to related combinatorial problems, correlation polytopes, and computational methods for high-dimensional discrete data (Fontana et al., 2017). The approach is foundational for simulation in risk aggregation, network modeling, and dependence-extremal analysis in statistics and applied probability.