Consider a set of random variables $\{X_i\}$ with joint pdf $f(x_1 ... x_n)$. Given the marginal pdfs $f_i(x_i)$, we can construct a joint distribution $$g(x_1 ... x_n) = \prod_i f_i(x_i)$$ which has the same marginals as $f$. (In particular, $g$ is the maximum entropy distribution satisfying this constraint.) $g$ is in some sense an approximation to $f$.

If we know the pairwise joint marginals $f_{ij}(x_i,x_j)$, is it possible to do something similar? Specifically:

  1. What conditions do the $f_{ij}$ need to satisfy in order for there to exist a $g$ giving pairwise joint marginals $f_{ij}$?
  2. Is there an explicit formula or simple algorithm that gives any such $g$?
  3. Is there an explicit formula that gives the maximum-entropy such $g$?

I'm also interested in generalizations to $k$-wise joint marginals, and the conditions under which the corresponding sequence of approximations converges to $f$, but that seems like a separate question. An answer for the special case where all the $f_{ij}$ are the same would also be interesting.

See also:
Constructing a joint distribution from pairwise bivariate marginal distributions? (No answers, but a comment mentions "pair-copula model")

  • $\begingroup$ The pairwise marginal densities need be compatible with at least one joint$$f_{ij}(x_i,x_j)=\int f(\mathbf x)\text d\mathbf x_{-(i,j)}$$ $\endgroup$
    – Xi'an
    Commented Jan 8, 2021 at 8:26
  • $\begingroup$ @Xi'an right, of course, I guess really I'm asking for a more explicit condition (efficiently computable, maybe). $\endgroup$
    – Daniel
    Commented Jan 9, 2021 at 5:43
  • 1
    $\begingroup$ I don't have a direct answer to your question, but the "pair-copula" construction mentioned in the other question is probably not what you want (e.g. for 3 variables, it constructs the joint from 3 pairs, but it's not $(X_1,X_2)$, $(X_2,X_3)$ and $(X_1,X_3)$, it's $(X_1,X_2)$, $(X_2,X_3)$, and $(X_1,X_3|X_2)$). $\endgroup$
    – Chris Haug
    Commented Jun 12, 2021 at 14:49
  • $\begingroup$ In practice, the copula approach mentioned above makes a lot of sense. See e.g. arxiv.org/abs/2102.02852 $\endgroup$
    – Björn
    Commented Aug 18, 2021 at 21:01

4 Answers 4


For discrete random variables the raking/iterative proportional fitting algorithm constructs a joint distribution if one exists (under some additional assumption about zero cells). It works for marginal distributions of any order, and not necessarily the same order for each margin.

IPF is (or was) used to fit loglinear models -- its guaranteed convergence rate is only linear, but it use less memory than Newton-type methods. Raking is used in survey statistics to match sample margins to known population margins.

I don't know if this generalises to continuous distributions.

  • $\begingroup$ The algorithm described in the Wikipedia article you link looks like it only enforces that the single-variable marginals match a target distribution. Is there some way of generalizing it to enforce a constraint on the two-variable joint marginals? $\endgroup$
    – Daniel
    Commented Jan 8, 2021 at 4:01
  • $\begingroup$ Yes. It's the same algorithm (classic IPF): you reweight the observations so one pairwise marginal distribution is correct and then iterate over all the pairwise margins. It wouldn't be much use if it only worked for univariate margins. There's an example here: stats.stackexchange.com/questions/37771/… $\endgroup$ Commented Jan 8, 2021 at 4:21
  • $\begingroup$ Ok, this makes sense, and probably one can do something vaguely similar with continuous distributions using a sequence of finitely-parameterized spaces. Unfortunately, I don't think this algorithm provides much insight into the structure of the resulting distribution, since it's nearly a black-box optimization algorithm. $\endgroup$
    – Daniel
    Commented Jan 14, 2021 at 18:19
  • 1
    $\begingroup$ In the discrete case, the resulting distribution is the log-linear model whose sufficient statistics are the supplied margins. $\endgroup$ Commented Jan 27, 2021 at 2:19
  • $\begingroup$ Ok, that's interesting, thanks! $\endgroup$
    – Daniel
    Commented Jan 30, 2021 at 0:10

I had the same question, and was surprised to learn the problem of combining multivariate marginal distributions to a full, joint distribution doesn't have a straightforward analytical solution, as it does with univariate marginals (assuming independence).

I don't know if iterative proportional fitting is the only or best solution, but it well worked for me. I'm posting some R code using the (unusually excellent) mipfp package with an example of implementation for posterity.

# simulate some survey data
# ethnicity (hispanic / not), sex (male / female), & race (white / other)
simdat <- array(rpois(n = 8, 10), dim= c(2, 2, 2), 
                dimnames = list(ethn = c("H", "NH"),
                                sex = c("M", "F"),
                                race = c("W", "O")))

# proportions table; compare to p.hat below
simpro <- simdat / sum(simdat) 

# all bivariate marginal distributions
sex_race <- apply(simdat, MARGIN = c(2, 3), sum)
ethn_race <- apply(simdat, MARGIN = c(1, 3), sum)
ethn_sex <- apply(simdat, MARGIN = c(1, 2), sum)
target_marginals <- list(sex_race, 

# set a seed (initial values) to all 1s 
# or use known / prior cell counts when available
seed <- (simdat + 1) / (simdat + 1)

# indicate to which indices bivariate marginals refer
# (see apply functions above. order matters!)
target_dims <- list(c(2, 3), c(1, 3), c(1, 2))

# run IPF algorithm
outp<- Ipfp(seed, target_dims, target_marginals,
               iter=50, print=TRUE, tol=1e-5)

outp$x.hat # estimated contingency table (array)
outp$p.hat # estimated multinomial probability distribution

I believe that by setting the initial values to 1, we are essentially assuming whatever independence is left over after specifying the bivariate marginals. Please correct me if I'm wrong!

Note that the algorithm converges in one iteration if you use simdat as the initial values, as expected!


$$ \newcommand{\avg}{\mathbb{E}} $$ The answers above are helpful for discrete distributions, but don't have much to say about continuous distributions. They also don't address the question of existence, which turns out to be pretty interesting. Here are a few things I have figured out:

  1. The pairwise marginals need to result in consistent univariate distributions, e.g. we must have $$ \int_y f_{XY}(x,y) = \int_z f_{XZ}(x,z) $$
    In addition, the covariance matrix of any set of functions of the variables, which has elements given by $$ M_{ij}^{(h_1 ... h_n)} = \avg[h_i(X_i)h_j(X_j)] $$ must be positive definite. Notice that $ M_{ij} $ can be computed entirely from the pairwise marginals. I suspect that these conditions are sufficient, but haven't found a way to prove it yet. Also, I think for the second condition it might suffice to check only the functions in any complete basis set, e.g. $h_i(x) = x^{n_i}$. This case looks a little like the Hamburger Moment Problem, and it's possible that a similar proof would be useful here.
  2. Under some regularity conditions, specifying all pairwise marginals is equivalent to specifying all two-variable moments $$ \mathbb{E}[x_i^n x_j^m] $$ The problem is then to make some choice for the remaining moments. This choice needs to satisfy positive-definiteness (and possibly other constraints I've missed?). This paper gives some conditions for the existence of a solution. Isserlis' theorem itself does this for the special case where the $ f_{ij} $ are Gaussian. It seems plausible that there's a generalization of Isserlis' theorem that will work for the general case, where we make some choice like like $$ \avg[X^2 Y^2 Z] = \frac{1}{3}\bigg(\avg[X^2 Y^2]\avg[Z] + \avg[X^2 Z]\avg[Y^2] + \avg[Y^2 Z]\avg[X^2]\bigg) $$
  3. The maximum-entropy solution is of the form $$ g(x_1 ... x_n) = \prod_{i,j} g_{ij}(x_i,x_j) $$ and any solution of this form is maximum-entropy. This can be shown using Lagrange multipliers, and is equivalent to the log-linear model observation above. Note that Isserlis' theorem does give the maximum-entropy solution in the Gaussian case.

The general framework is to select from the set of all joint distributions matching a set of statistics the unique distribution with maximum entropy.

For pair-wise covariance statistics, see:
Dempster, Arthur P. Covariance selection. Biometrics (1972): 157-175.

For arbitrary "filters" (functions of many variables), see:
Zhu, S. C., Y. Wu, and D. Mumford. FRAME: Filters, Random Fields and Maximum Entropy - Towards a Unified Theory for Texture Modeling. IEEE CVPR 1993.

(Note similarities in approach. Dempster, Mumford, Wu and Zhu were all at Harvard in the early 1990s.)


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.