# Reconstruct multivariate binary distribution from marginals

I'm have a random vector $\bf a$ with binary entries, $a_i \in \{0,1\}$. The probability distribution $P({\bf a})$ is not fully specified, but I have the marginals $p_i$, which are the probabilities that the $i$-th entry is 1. Additionally, I also have the second-order marginals $p_{ij}$, which give the probabilities that both $a_i = 1$ and $a_j = 1$ at the same time. I'm looking for a way to derive a consistent probability mass function for $P({\bf a})$.

Generally, not all pairs of $p_i$ and $p_{ij}$ lead to a consistent distribution $P(\bf a)$, but in my case I derive them from e.g. numerical simulations of an underlying $P({\bf a})$ and can thus be sure that a solution must exist. Of course, there could be multiple $P({\bf a})$ consistent with the requirements. In this case I'd like to have the distribution with maximal entropy, although an approximation would be good enough.

If only $p_i$ was given, or if $p_{ij} = p_ip_j$, I'm pretty sure the solution would read $P({\bf a}) = \prod_i p_i^{a_i} (1 - p_i)^{1 - a_i}$, where I defined $0^0=1$. However, I'm not sure what the solution looks like with a more general $p_{ij}$.

Does anybody know the solution, a numerical method, or search terms I could use to find out more about the problem?