I'm looking for a "nice" way to parametrize the joint distribution of multiple, possibly correlated discrete random variables on {0,1}. Even for N=2, there doesn't seem to be an obvious way to do it.
The kind of thing I'm looking for, is how bi-normal distributions can be simply parametrized by the means and variances of the individual distributions and by a single "correlation" parameter. Things get more complicated in higher dimensions of course, but is there such a nice analog in my case?
I have no deep insight to share for general $N$, but for the simplest case of $N=2$ (arbitrary coupling / joint distribution of two fixed binary discrete uniform random variables), there's a 1-to-1 correspondence between possible values of the correlation coefficient and possible joint distributions/couplings.
Cf. this related question on CrossValidated (statistics stackexchange). In this very special case nice parameterizations in terms of the correlation do exist.
Specifically, given discrete $N=2$ binary RVs $X, Y$ with $P(X=0, Y=0):=p_{00}$, $P(X=0, Y=1):=p_{01}$, $P(X=1, Y=0):=p_{10}$, and $P(X=1, Y=1):=p_{11}$, it's a perhaps tedious but not difficult computation to show that the covariance is $$p_{11}p_{00} - p_{01}p_{10} \,. $$ Interestingly this formula for the correlation in terms of the joint distribution does not depend on the particular marginal distributions of $X$ and $Y$.
Given $N=2$ discrete $M$-ary random variables $X, Y$ (i.e. both having "state space" of cardinality $M$), there are $2M - 1$ linearly independent constraints on the $M^2$ parameters defining an arbitrary joint distribution/coupling, so the space of all possible joint distributions/couplings corresponds to $M^2 - 2M + 1 = (M-1)(M-1)$ free parameters.
Note that each of the two $M$-ary random variables can be completely charactetrized by $(M-1)$ binary random variables (the indicator RVs for the first $M-1$ possible values, with the indicator RV for the $M$th possible value obviously being a linear combination of the others).
Note also that in the case that $M=2$, $(M-1)^2 = 1$.
So I conjecture (based on these combinatorial/dimensional arguments and because it seems intuitively plausible) that the joint distribution of an $M$-ary discrete RV $X$ and an $M$-ary discrete RV $Y$ can be completely characterized by the $(M-1)^2$ joint distributions of the indicator RV pairs $$(I\{X=1\}, I\{Y=1\}), (I\{X=2\}, I\{Y=1\}), \dots, (I\{X=M-1\}, I\{Y=1\}), (I\{X=1\}, I\{Y=2\}) , \dots, (I\{X=M-1\}, I\{Y=2\}), \dots, (I\{X=M-1\}, I\{Y=M-1\}). $$ Given how each of these $(M-1)^2$ joint distributions can each be completely characterized by a single parameter (e.g. the corresponding covariance or correlation), I strongly suspect that the $(M-1)^2$ parameters required to characterize the original joint distribution of $X$ and $Y$ can be chosen to be such covariances or correlations of indicator RV pairs.
Note that this is basically the same idea as that behind the "one-vs-rest" approach to multiclass classification in statistics/machine learning.
Also the arguments leading to the above don't actually seem to require the RVs $X$ and $Y$ to have the same number of states, i.e. I think it's appropriate to also conjecture that if $X$ has $M_X$ possible states and $Y$ has $M_Y$ possible states, then exactly $(M_X - 1)(M_Y - 1)$ parameters should be needed to characterize the space of all possible joint distributions / couplings of those two random variables, and that the parameters again could probably be chosen to be "one-vs-rest" covariances or correlations as discussed above.
So TL;DR I agree with Douglas Zare that first and second moments should be sufficient in the case that $N=2$. I have done some of the tedious linear algebra calculations to substantiate some of these claims over the weekend, but have a non-math-related day job, am not strong at calculations, and so probably won't be able to write out on paper (much less LaTeX up) full descriptions of calculations/proofs of all of these results any time soon. But hopefully I stated the conjectures cleanly enough that a motivated reader could straightforwardly verify or falsify these claims themselves.
For $N \ge 3$, the niceness of being able to use 2nd moments to characterize joint distributions for $N=2$ would naturally lead one to conjecture that $N$th order tensors of $N$th order moments (possibly again of "one-vs-rest" indicators) are always sufficient. But that is conjecture mostly based on wishful thinking, rather than evidence. Anyway the multilinear algebra calculations involved are surely more tedious. The obvious starting example to investigate first would be the space of all possible couplings of three uniform binary RVs -- I hope to do that myself one day and report back.
In any case, given that all of the computations are (multi)linear, they seem to be special cases of the kinds of computations that people do in algebraic statistics where restricting to joint distributions with certain (marginal or conditional) independences leads to multiplicative constraints and thus solution spaces that are algebraic varieties corresponding to higher degree polynomials, rather than the linear spaces discussed above. (Cf. e.g. Example 3.1.6 in Drton's Lectures on Algebraic Statistics, PDF link.) So I'm not sure whether the questions being asked are too trivial to be discussed in most of the introductory books on algebraic statistics that are out there. That being said, I would encourage the motivated reader to look there for possible relevant references.
P.S. This kind of seems similar to the "mean parameterization" of exponential families, for definitions cf. section 1.10 of these notes, section 4 of these notes, or section 8.4 of these notes. I am not sure whether exponential families are actually related to this problem, however, so don't quote me on this (so to speak).
The probabilities on $\{0,1\}^N$ have to add up to $1$, so they give the barycentric coordinates of a point in a simplex of dimension $2^N-1$. For $N=2$, you can use other coordinates, including triples
$(E(X),E(Y), E(XY))$
or
$(E(X),E(Y), \text{cov}(X,Y))$
where $\text{cov}(X,Y)=E(XY)-E(X)E(Y)) = \rho(X,Y)\sigma(X)\sigma(Y)$.
