A Practical Guide to Simulating Correlated Binary Outcomes

Chi Heem Wong; Zied Ben Chaouch

REVIEW 2 major objections 5 minor 29 references

The paper claims that simulating Bernoulli outcomes with prescribed means and pairwise correlations is exactly decidable by a linear program over the 2^N atomic probabilities, which returns an exact law when feasible and an infeasibility ce

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · deepseek-v4-flash

2026-08-01 19:53 UTC pith:QPDZRHIW

load-bearing objection A useful and mathematically sound clarification of when Gaussian-threshold methods are exact, packaged with an exact LP feasibility test for moderate N; the issues are reproducibility details, not conceptual flaws. the 2 major comments →

arxiv 2607.16801 v1 pith:QPDZRHIW submitted 2026-07-18 cs.CE q-fin.CPq-fin.RMstat.CO

A Practical Guide to Simulating Correlated Binary Outcomes

Chi Heem Wong , Zied Ben Chaouch This is my paper

classification cs.CE q-fin.CPq-fin.RMstat.CO MSC 62H2090C0565C10

keywords correlated Bernoullibinary outcomeslinear programmingfeasibility certificateGaussian thresholdtetrachoric calibrationcross-moment matchingsimulation

verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper establishes an exact feasibility test for the common risk-modeling task of generating dependent binary events with given means and pairwise correlations. It shows that the usual Gaussian-threshold shortcut — setting latent normal correlations to the target Bernoulli correlations, or calibrating them pairwise — is generally incorrect, because thresholding distorts correlations and the calibrated latent matrix may not be a valid covariance. Instead, the paper formulates the problem directly over the joint Bernoulli distribution: a linear program whose variables are the 2^N atomic probabilities, with linear constraints for normalization, nonnegativity, means, and pairwise cross-moments. The LP either recovers an exact law or certifies infeasibility, and the paper proves that every feasible target has a compatible law on at most 1+N+C(N,2) states while every infeasible target has a separating quadratic witness. This matters because it separates the feasibility question from the choice of model family and gives practitioners a principled way to know when their target moments are simulable at all.

Core claim

The central discovery is that Bernoulli feasibility with prescribed first and second moments is characterizable as membership in the convex hull of the binary moment feature vectors, and that this characterization is computationally actionable as a linear program. Specifically, the PMF-LP minimizes the zero function subject to 1ᵀα=1, α≥0, Mᵀα=μ, and κᵢⱼᵀα=Σᵢⱼ+μᵢμⱼ; a feasible solution is an exact law and an infeasible problem yields a separating quadratic certificate via the separating hyperplane theorem. The paper further shows that Gaussian-threshold constructions (CoMBO) are exact only when the tetrachorically calibrated latent matrix is positive semidefinite, and that even feasible Berno

What carries the argument

The central object is the joint Bernoulli probability mass function vector α indexed by the 2^N outcomes, with feature map φ(x)=(1, x_i, x_i x_j) for each state. The workhorse identity is that means and pairwise cross-moments are linear functions of α (Mᵀα and κᵢⱼᵀα), which turns the moment-matching problem into a linear feasibility program. The supporting machinery is the convex-hull/separating-hyperplane characterization: feasibility is membership of the target vector b in conv{φ(x)}, which yields both the sparse-support bound (Carathéodory) and the infeasibility certificate (a quadratic function nonnegative on all states but negative at the target). For the truncated-moment extension, Möb

Load-bearing premise

The load-bearing premise is that the linear program is solved in exact real arithmetic; in floating-point implementations, solver tolerances can produce tiny negative probabilities that are clipped to zero, which could mis-certify feasibility or return a law that violates the moment equations on boundary or ill-conditioned targets.

What would settle it

A concrete test: for the N=3 target in Appendix C (μ=(0.3,0.5,0.6), R with ρ12=0.2, ρ13=0.3, ρ23=0.6), solve the PMF-LP in rational or high-precision arithmetic and confirm the returned α recovers the target moments exactly. Then repeat with a boundary target such as μ_i=0.5 for all i and ρ_ij=-1 for all pairs; the exact LP must certify infeasibility and the separating quadratic must be negative at the target. If any solver returns a 'feasible' law for the second case, or an infeasibility certificate for a genuinely feasible target, the practical exactness claim fails.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

If this is right

If the PMF-LP is solved exactly, a practitioner obtains a definitive yes/no answer to whether the requested moments are simulable, without committing to a latent or parametric family.
For any feasible target, there exists an exact simulator supported on at most 1+N+C(N,2) states, so sparse-support methods are not hopeless even though worst-case complexity remains exponential.
The Gaussian-threshold workflow (CoMBO) is exact only when the pairwise tetrachoric calibration yields a positive semidefinite matrix; otherwise any PSD repair changes at least one induced Bernoulli correlation.
The order-k truncated moment completion provides a reduced model class with polynomial-time sampling for fixed k, but it is not a universal feasibility test: infeasibility may reflect an overly small truncation order.
The pairwise-only (k=2) construction has explicit closed-form atoms, and its failure to satisfy nonnegativity does not imply global Bernoulli infeasibility.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

A natural extension is to use the separating quadratic certificate as a diagnostic: the coefficients (c, a_i, a_ij) indicate which moment constraints conflict most, potentially guiding users to relax specific pairwise correlations rather than abandoning the target.
The convex-hull support bound suggests a column-generation approach where the pricing subproblem is a binary quadratic problem; if that QUBO can be solved or approximated well, exact laws for much larger N may be reachable than full enumeration allows.
The paper's exactness claim is in real arithmetic; in floating-point practice, boundary targets may be mis-certified, so a rational-arithmetic or tolerance-aware implementation would be a direct testable refinement.
The same PMF-LP machinery could be extended to higher-order moments (e.g., third-order cross-moments) as additional linear constraints, giving a flexible way to fit richer dependence structures beyond pairwise.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Desk Editor's Note

A useful and mathematically sound clarification of when Gaussian-threshold methods are exact, packaged with an exact LP feasibility test for moderate N; the issues are reproducibility details, not conceptual flaws.

read the letter

Candid take: the paper's central mathematical content is right, and it gives practitioners something they genuinely need—a clear statement of when Gaussian-threshold (CoMBO) simulation is exact and an exact LP check for Bernoulli moment feasibility. Send it to a knowledgeable referee.

The genuinely new pieces are the PSD condition for calibrated CoMBO exactness, the PMF-LP as the universal feasibility test, the support bound of 1+N+C(N,2) with a separating quadratic certificate, and the truncated-moment completion. The LP formulation itself is natural—it is the obvious LP over the 2^N atoms—but the paper states it precisely and draws the consequences. The derivations I checked are correct: tetrachoric calibration is monotone, the Möbius identities are standard, and the convex-geometric certificate is clean. The paper is also honest about exponential complexity and about the truncated method being a model class, not a universal fix. That deserves credit.

The soft spots are real but minor. Appendix B reports a 'calibrated Bernoulli' simulation without stating how the non-PSD calibrated matrix was repaired, so the headline illustration is not reproducible from the text. That is a concrete omission for a methods paper, and the fix is simple: state the repair rule or drop the table. The code repository has no commit hash, which makes the practical claims harder to pin down. As to the stress-test complaint: I don't buy that the numerical certification gap is load-bearing. The PMF-LP's exactness is a real-arithmetic statement, like every LP-based construction; the paper already acknowledges tolerances and clipping. A note on checking moment equalities after solving and on extracting Farkas certificates from standard solvers would strengthen the practical claim, but the absence does not undercut the theory.

The central argument holds up. The paper is a methods clarification rather than a breakthrough; its audience is risk modelers and anyone simulating dependent binaries. I would cite it for the PSD condition and the feasibility-LP results. It deserves peer review, with the reproducibility items addressed. I would accept after minor revision.

Referee Report

2 major / 5 minor

Summary. The paper addresses simulation of N dependent Bernoulli outcomes with prescribed means and pairwise Pearson correlations. It critiques the Gaussian-threshold construction (CoMBO), showing that naive CoMBO is generally biased and that pairwise tetrachoric calibration is exact only when the assembled latent correlation matrix is positive semidefinite. The main contribution is a direct linear program over the 2^N atomic probabilities of the joint Bernoulli PMF (PMF-LP) that matches first and second moments, with the claimed properties of exact feasibility characterization, infeasibility certificates, a support bound of 1+N+C(N,2) states for feasible targets, and a separating quadratic certificate for infeasible targets. The paper also develops a truncated-moment completion scheme with sequential conditioning and a sparse-support working-set refinement. The mathematical core is largely correct, but the advertised exactness and certification guarantees are stronger than what the implemented floating-point pipeline actually delivers, and the main empirical illustration omits a key implementation detail.

Significance. The PMF-LP formulation is a clean and useful contribution: it separates Bernoulli feasibility from representability within a latent Gaussian family, gives a rigorous convex-hull characterization, and provides concrete certificates. The proof of the support bound and the separating certificate is sound, and the clarification of when calibrated CoMBO is exact is valuable for practitioners. The paper also provides code, which aids reproducibility. The two extensions (truncated completion and sparse support) are reasonable but are model-class tools rather than exact feasibility tests, and the paper mostly acknowledges this. If the numerical-exactness caveats are resolved, this would be a solid practical and conceptual contribution to the simulation/risk-modeling literature.

major comments (2)

[Abstract and §4.2, Eq. (13)] The advertised deliverables—'returns an exact law ... or certifies infeasibility'—are not guaranteed by the implementation described. The feasibility LP is solved with a floating-point solver; §4.3 and App. D.8 explicitly allow clipping 'tiny negative values caused by solver tolerances to zero', which can violate the equality constraints (13d)–(13e). No a posteriori verification of the returned α against the moment equalities is described, and no procedure is given for extracting a Farkas/separating certificate from the solver when the LP is declared infeasible. For targets on the boundary of the feasible polytope or ill-conditioned instances, a floating-point solver can return α with nonzero residuals or mis-certify infeasibility. Because the exactness claim is the paper's central novelty, this needs to be fixed: either qualify the claims as 'up to solver tolerance' throughout the abstr
[Appendix B, Table 2] The 'calibrated Bernoulli' experiment is not reproducible as reported. The computed R*_lat has determinant ≈ −0.419 < 0 and is not PSD, so one cannot sample directly from N(0, R*_lat). The table reports calibrated Bernoulli correlations that differ from the targets, but the repair/regularization used to create a PSD matrix is not stated. This matters because the experiment is meant to illustrate that 'pairwise calibration alone does not yield an exact realizable multivariate Bernoulli generator.' Without specifying the repair (projection, eigenvalue clipping, etc.) and its effect, the reader cannot tell whether the residual error is due to the mathematical obstruction or to the arbitrary repair. Please state the repair procedure and its impact on the reported numbers, or replace the table with a deterministic demonstration.

minor comments (5)

[§3.2, Eq. (9)] The calibration interval should be r∈[−1,1] or endpoints should be handled separately; for a target correlation at the boundary of the admissible range, the unique latent correlation is ±1, not in (−1,1).
[§1 and §4.4] 'The accompanying code can be assessed at ...' should be 'accessed'; also consider including a versioned snapshot of the repository for reproducibility.
[References] Reference [17] lists the author as 'PK Triverdi'; the correct spelling is 'Trivedi' (P.K. Trivedi and D.M. Zimmer).
[Appendix B] The sentence 'Thus pairwise calibration thus provides' contains a duplicated 'thus'.
[Algorithms 1 and 2] Algorithms 1 and 2 are nearly identical; consider consolidating them to avoid redundancy.

Circularity Check

0 steps flagged

No significant circularity: the PMF-LP and its feasibility certificates are derived directly from probability axioms and external convex-geometric results, with no fitted input relabeled as prediction.

full rationale

The paper's central construction, the PMF-LP (Section 4.2, Eqs. 13a-13e), is obtained by writing normalization, nonnegativity, mean constraints, and pairwise cross-moment constraints as linear equalities in the 2^N atomic probabilities. The target means and correlations enter only as the right-hand sides of these constraints; feasibility of the LP is exactly membership of the target vector b in conv{phi(x)}, as proved in Appendix D, Proposition D.7 using Carathéodory and the strict separating-hyperplane theorem. No parameter is fitted to a subset of data and then used to predict a closely related quantity: the LP either produces one compatible law or an infeasibility witness, and the paper explicitly disclaims the truncated completion as a model class rather than a universal test ('infeasibility may mean either that the target moments are globally impossible or that the chosen truncation order is too small', Sec 4.3). The only self-reference is the authors' code repository (Sec 4.4), which is not load-bearing for the mathematical derivation. The numerical-tolerance caveat in Sec 4.3/D.8 ('solver tolerances can produce tiny negative values that are clipped to zero') is a floating-point reliability concern, not a circularity; the advertised exactness is conditional on solver accuracy, but that is a correctness risk outside the derivation-chain definition. External citations (Sklar, Chaganty-Joe, Emrich-Piedmonte, Qaqish, Hoeffding, etc.) are standard and do not import the paper's conclusions. Hence score 0.

Axiom & Free-Parameter Ledger

2 free parameters · 6 axioms · 0 invented entities

The paper introduces no new physical or mathematical entities. Its free choices are the truncation order k (a model-class knob) and an unstated PSD repair in the illustrative simulation. All derivations rely on standard probability, convex analysis, and Boolean-lattice identities.

free parameters (2)

truncation order k = user-chosen (>=2)
The reduced moment-completion family keeps moments up to order k and zeroes the rest; the choice of k determines whether a feasible law exists within the truncated family, and D.12 shows k must be Omega(N) when the mean sum is Theta(N).
unstated PSD repair in Appendix B = not specified
The 'calibrated Bernoulli' simulation uses a latent matrix R*_lat that is not PSD; to generate draws the authors must have replaced it with a nearby PSD matrix, but the repair rule is not described, so the reported numbers embed an implicit free choice.

axioms (6)

standard math Standard probability axioms and the Bernoulli state space {0,1}^N with a PMF that is a convex combination of atoms.
Used in Section 4.1 to write the LP (13) over atomic probabilities; any joint law can be represented this way.
standard math Sklar's theorem and Gaussian copula machinery for the CoMBO description.
Invoked in Section 3.1 to define the copula and the Gaussian-threshold construction; standard.
standard math Carathéodory's theorem for the support bound m = 1+N+C(N,2).
Used in Proposition D.7 to bound the support of a feasible law.
standard math Möbius inversion on the Boolean lattice (zeta transform).
Used throughout Appendix D and Section 4.3 to convert between atoms and cross-moments.
domain assumption The target means and correlations satisfy the necessary pairwise bounds of Eq. (3).
The PMF-LP uses these only as pre-checks; the LP is the global test. The user is assumed to supply a symmetric correlation matrix with unit diagonal.
ad hoc to paper For the truncated completion, the model class assumption gamma_A=0 for |A|>k.
This is an explicit modeling restriction introduced in Section 4.3; it is not a property of the target. Infeasibility in the truncated LP may reflect this restriction rather than global infeasibility.

pith-pipeline@v1.3.0-alltime-deepseek · 20994 in / 20751 out tokens · 187771 ms · 2026-08-01T19:53:46.336486+00:00 · methodology

0 comments

read the original abstract

Simulating dependent Bernoulli outcomes with prescribed means and pairwise Pearson correlations is a common task in risk modeling. A familiar approach is the Gaussian-threshold workflow for binary outcomes, often viewed as a Bernoulli analogue of the Gaussian copula construction. We show that setting latent Gaussian correlations equal to target Bernoulli correlations is generally incorrect after thresholding, and that pairwise tetrachoric calibration is exact only when the calibrated latent matrix is positive semidefinite. We therefore formulate the problem directly over the joint Bernoulli probability mass function. Given target means and pairwise correlations, we impose normalization, nonnegativity, mean constraints, and pairwise cross-moment constraints as a linear program over the $2^N$ atomic probabilities. The resulting PMF formulation either returns an exact law matching the requested first and second moments or certifies infeasibility. A convex-hull characterization further shows that every feasible target admits a law supported on at most $1+N+\binom{N}{2}$ states, while every infeasible target admits a separating quadratic certificate. We then develop a truncated-moment completion scheme that fits a reduced cross-moment table and generates samples by sequential conditioning, together with a sparse-support working-set refinement that can reduce memory usage on structured instances, although the worst-case complexity remains exponential. Together, these constructions provide an exact PMF-based framework for feasibility and simulation at moderate dimension and structured alternatives when the full atomic representation is impractical, while clarifying the limits of Gaussian-threshold constructions.

discussion (0)

Reference graph

Works this paper leans on

29 extracted references

[1]

A representation of the joint distribution of responses to n dichotomous items

Raghu Raj Bahadur. A representation of the joint distribution of responses to n dichotomous items. The Annals of Mathematical Statistics, 1961

1961
[2]

Cario and Barry L

Celso M. Cario and Barry L. Nelson. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Technical report, Department of Industrial Engineering and Management Sciences, Northwestern University, 1997. Technical Report

1997
[3]

Chaganty and Harry Joe

Narasimhan R. Chaganty and Harry Joe. Range of correlation matrices for dependent Bernoulli random variables.Biometrika, 2006

2006
[4]

Solving linear programs in the current matrix multiplication time.Journal of the ACM (JACM), 68(1):1–39, 2021

Michael B Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time.Journal of the ACM (JACM), 68(1):1–39, 2021

2021
[5]

Emrich and Marion R

Laurie J. Emrich and Marion R. Piedmonte. A method for generating high-dimensional correlated binary variates with specified marginal means and correlations.The American Statistician, 45(4):302–304, 1991

1991
[6]

Representation of multivariate Bernoulli distributions with a given set of specified moments.Journal of Multivariate Analysis, 168:290–303, 2018

Roberto Fontana and Patrizia Semeraro. Representation of multivariate Bernoulli distributions with a given set of specified moments.Journal of Multivariate Analysis, 168:290–303, 2018

2018
[7]

Probability inequalities for sums of bounded random variables.Journal of the American Statistical Association, 58(301):13–30, 1963

Wassily Hoeffding. Probability inequalities for sums of bounded random variables.Journal of the American Statistical Association, 58(301):13–30, 1963

1963
[8]

Iman and William J

Ronald L. Iman and William J. Conover. A distribution-free approach to inducing rank correlation among input variables.Communications in Statistics - Simulation and Computation, 11(3):311–334, 1982

1982
[9]

CRC Press, 2014

Harry Joe.Dependence Modeling with Copulas. CRC Press, 2014

2014
[10]

Joseph B. Kadane. Sums of possibly associated Bernoulli variables: The Conway–Maxwell–Binomial distribution.Bayesian Analysis, 11(2):403–420, 2016

2016
[11]

Nelsen.An Introduction to Copulas

Roger B. Nelsen.An Introduction to Copulas. Springer, 2006

2006
[12]

Consistent families of measures and their extensions.Theory of Probability and its Applications, 7(2):147–163, 1962

Vorob’ev Nikolai N. Consistent families of measures and their extensions.Theory of Probability and its Applications, 7(2):147–163, 1962

1962
[13]

Tree bounds for sums of Bernoulli random variables: A linear optimization approach.INFORMS Journal on Optimization, 3(1):23–45, 2021

Divya Padmanabhan and Karthik Natarajan. Tree bounds for sums of Bernoulli random variables: A linear optimization approach.INFORMS Journal on Optimization, 3(1):23–45, 2021

2021
[14]

Bahjat F. Qaqish. A family of multivariate binary distributions for simulating correlated binary variables with specified marginal means and correlations.Biometrika, 90(2):455–463, 2003

2003
[15]

Fonctions de r´ epartition ` a n dimensions et leurs marges

M Sklar. Fonctions de r´ epartition ` a n dimensions et leurs marges. InAnnales de l’ISUP, volume 8, pages 229–231, 1959

1959
[16]

Jozef L. Teugels. Some representations of the multivariate Bernoulli and binomial distributions.Journal of Multivariate Analysis, 32(2):256–268, 1990

1990
[17]

Copula modeling: An introduction for practitioners.Foundations and Trends in Econometrics, 1(1):1–111, 2005

PK Triverdi and David M Zimmer. Copula modeling: An introduction for practitioners.Foundations and Trends in Econometrics, 1(1):1–111, 2005. 12 A. Comparison of the Main Constructions Table 1 summarizes the main constructions discussed in the paper and highlights the trade-offs among exactness, representation, and computational cost. The point is not that...

2005
[18]

Setγ ∅ = 1,γ{i} =µ i, and γ{i,j} =µ iµj +ρij √ µi(1−µ i)µj(1−µ j) (1≤i<j≤N)
[19]

Introduce decision variablesγ A for allA⊆[N] with 3≤|A|≤k, and setγ A = 0 for all|A|>k
[20]

Solve the linear feasibility problem pA = ∑ B⊇A |B|≤k (−1)|B|−|A|γB ≥0 for allA⊆[N] with|A|≤k, together with any chosen bounds on the unknown higher-order moments
[21]

Otherwise store the fitted truncated moment table Γk ={γ A :|A|≤k}

If the LP is infeasible, either stop or increase k. Otherwise store the fitted truncated moment table Γk ={γ A :|A|≤k}. Online sampling step
[22]

For each sampleℓ= 1,...,L, initializeO 0 =Z 0 =∅
[23]

Form= 0,1,...,N−1, compute w(Om,Zm) = ∑ B⊆Zm |Om|+|B|≤k (−1)|B|γOm∪B and qm = ∑ B⊆Zm |Om|+1+|B|≤k (−1)|B|γOm∪{π(m+1)}∪B ∑ B⊆Zm |Om|+|B|≤k (−1)|B|γOm∪B
[24]

DrawU m∼Unif(0,1) and setX π(m+1) = 1 ifU m≤q m, otherwise setX π(m+1) = 0
[25]

Output:drawsX (1),...,X (L)∈{0,1} N from the fitted order-ktruncated model

Update Om+1 =O m∪{π(m+ 1) :X π(m+1) = 1}, Z m+1 =Z m∪{π(m+ 1) :X π(m+1) = 0}, and continue until all coordinates are assigned. Output:drawsX (1),...,X (L)∈{0,1} N from the fitted order-ktruncated model. Any valid law requires botht=p {1,2,3}≥0 andp {i,j}≥0, so necessarilyt= 0. But then p∅ = 1− ∑ i γ{i} + ∑ i<j γ{i,j}−t=− 1 2 <0. Thus no choice of third-or...
[26]

InitializeS 0 with random and/or structured states
[27]

Solve the restricted master LP
[28]

Generate one or more new states via the pricing subproblem
[29]

This refinement does not remove worst-case exponential behavior, but it can reduce memory use and runtime substantially on structured instances

Add them toS t and repeat until slack is zero or progress stalls. This refinement does not remove worst-case exponential behavior, but it can reduce memory use and runtime substantially on structured instances. A related idea is to sample a subset of variables first (for example X1,...,X k) and then solve an LP for the remainder. In general, however, that...

[1] [1]

A representation of the joint distribution of responses to n dichotomous items

Raghu Raj Bahadur. A representation of the joint distribution of responses to n dichotomous items. The Annals of Mathematical Statistics, 1961

1961

[2] [2]

Cario and Barry L

Celso M. Cario and Barry L. Nelson. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Technical report, Department of Industrial Engineering and Management Sciences, Northwestern University, 1997. Technical Report

1997

[3] [3]

Chaganty and Harry Joe

Narasimhan R. Chaganty and Harry Joe. Range of correlation matrices for dependent Bernoulli random variables.Biometrika, 2006

2006

[4] [4]

Solving linear programs in the current matrix multiplication time.Journal of the ACM (JACM), 68(1):1–39, 2021

Michael B Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time.Journal of the ACM (JACM), 68(1):1–39, 2021

2021

[5] [5]

Emrich and Marion R

Laurie J. Emrich and Marion R. Piedmonte. A method for generating high-dimensional correlated binary variates with specified marginal means and correlations.The American Statistician, 45(4):302–304, 1991

1991

[6] [6]

Representation of multivariate Bernoulli distributions with a given set of specified moments.Journal of Multivariate Analysis, 168:290–303, 2018

Roberto Fontana and Patrizia Semeraro. Representation of multivariate Bernoulli distributions with a given set of specified moments.Journal of Multivariate Analysis, 168:290–303, 2018

2018

[7] [7]

Probability inequalities for sums of bounded random variables.Journal of the American Statistical Association, 58(301):13–30, 1963

Wassily Hoeffding. Probability inequalities for sums of bounded random variables.Journal of the American Statistical Association, 58(301):13–30, 1963

1963

[8] [8]

Iman and William J

Ronald L. Iman and William J. Conover. A distribution-free approach to inducing rank correlation among input variables.Communications in Statistics - Simulation and Computation, 11(3):311–334, 1982

1982

[9] [9]

CRC Press, 2014

Harry Joe.Dependence Modeling with Copulas. CRC Press, 2014

2014

[10] [10]

Joseph B. Kadane. Sums of possibly associated Bernoulli variables: The Conway–Maxwell–Binomial distribution.Bayesian Analysis, 11(2):403–420, 2016

2016

[11] [11]

Nelsen.An Introduction to Copulas

Roger B. Nelsen.An Introduction to Copulas. Springer, 2006

2006

[12] [12]

Consistent families of measures and their extensions.Theory of Probability and its Applications, 7(2):147–163, 1962

Vorob’ev Nikolai N. Consistent families of measures and their extensions.Theory of Probability and its Applications, 7(2):147–163, 1962

1962

[13] [13]

Tree bounds for sums of Bernoulli random variables: A linear optimization approach.INFORMS Journal on Optimization, 3(1):23–45, 2021

Divya Padmanabhan and Karthik Natarajan. Tree bounds for sums of Bernoulli random variables: A linear optimization approach.INFORMS Journal on Optimization, 3(1):23–45, 2021

2021

[14] [14]

Bahjat F. Qaqish. A family of multivariate binary distributions for simulating correlated binary variables with specified marginal means and correlations.Biometrika, 90(2):455–463, 2003

2003

[15] [15]

Fonctions de r´ epartition ` a n dimensions et leurs marges

M Sklar. Fonctions de r´ epartition ` a n dimensions et leurs marges. InAnnales de l’ISUP, volume 8, pages 229–231, 1959

1959

[16] [16]

Jozef L. Teugels. Some representations of the multivariate Bernoulli and binomial distributions.Journal of Multivariate Analysis, 32(2):256–268, 1990

1990

[17] [17]

Copula modeling: An introduction for practitioners.Foundations and Trends in Econometrics, 1(1):1–111, 2005

PK Triverdi and David M Zimmer. Copula modeling: An introduction for practitioners.Foundations and Trends in Econometrics, 1(1):1–111, 2005. 12 A. Comparison of the Main Constructions Table 1 summarizes the main constructions discussed in the paper and highlights the trade-offs among exactness, representation, and computational cost. The point is not that...

2005

[18] [18]

Setγ ∅ = 1,γ{i} =µ i, and γ{i,j} =µ iµj +ρij √ µi(1−µ i)µj(1−µ j) (1≤i<j≤N)

[19] [19]

Introduce decision variablesγ A for allA⊆[N] with 3≤|A|≤k, and setγ A = 0 for all|A|>k

[20] [20]

Solve the linear feasibility problem pA = ∑ B⊇A |B|≤k (−1)|B|−|A|γB ≥0 for allA⊆[N] with|A|≤k, together with any chosen bounds on the unknown higher-order moments

[21] [21]

Otherwise store the fitted truncated moment table Γk ={γ A :|A|≤k}

If the LP is infeasible, either stop or increase k. Otherwise store the fitted truncated moment table Γk ={γ A :|A|≤k}. Online sampling step

[22] [22]

For each sampleℓ= 1,...,L, initializeO 0 =Z 0 =∅

[23] [23]

Form= 0,1,...,N−1, compute w(Om,Zm) = ∑ B⊆Zm |Om|+|B|≤k (−1)|B|γOm∪B and qm = ∑ B⊆Zm |Om|+1+|B|≤k (−1)|B|γOm∪{π(m+1)}∪B ∑ B⊆Zm |Om|+|B|≤k (−1)|B|γOm∪B

[24] [24]

DrawU m∼Unif(0,1) and setX π(m+1) = 1 ifU m≤q m, otherwise setX π(m+1) = 0

[25] [25]

Output:drawsX (1),...,X (L)∈{0,1} N from the fitted order-ktruncated model

Update Om+1 =O m∪{π(m+ 1) :X π(m+1) = 1}, Z m+1 =Z m∪{π(m+ 1) :X π(m+1) = 0}, and continue until all coordinates are assigned. Output:drawsX (1),...,X (L)∈{0,1} N from the fitted order-ktruncated model. Any valid law requires botht=p {1,2,3}≥0 andp {i,j}≥0, so necessarilyt= 0. But then p∅ = 1− ∑ i γ{i} + ∑ i<j γ{i,j}−t=− 1 2 <0. Thus no choice of third-or...

[26] [26]

InitializeS 0 with random and/or structured states

[27] [27]

Solve the restricted master LP

[28] [28]

Generate one or more new states via the pricing subproblem

[29] [29]

This refinement does not remove worst-case exponential behavior, but it can reduce memory use and runtime substantially on structured instances

Add them toS t and repeat until slack is zero or progress stalls. This refinement does not remove worst-case exponential behavior, but it can reduce memory use and runtime substantially on structured instances. A related idea is to sample a subset of variables first (for example X1,...,X k) and then solve an LP for the remainder. In general, however, that...