Joris,

So we assume we know only the marginals $p(x_i)$ and the probabilities that $p(x_i=x_j)$. In terms of the physics' "spin" notation, $s_i=\pm 1$, this means that we know $\left<s_i \right> \equiv \sum_{s_i} s_i p(s_i)$ and $\left<s_i s_j \right> \equiv \sum_{s_i s_j} s_i s_j p(s_i,s_j)$.

There are many different joint probability distributions $p(s) \equiv p(s_1, \ldots, s_N)$ that have these single node and pairwise statistics. Indeed Egosphere only asks for *some* distribution.

One common approach to model distributions with partial knowledge in the form of constraints is to consider a special distribution, namely the maximum entropy (maxent) distribution.

Here, the maxent distribution $p^*(s)$ maximizes the entropy
$$
H = -\sum_{s} p^*(s) \log p^*(s)
$$
under the constraints $$\sum_{s_i} s_i p^*(s_i) = \left<s_i \right> \forall i~\mbox{and}~ \sum_{s_i s_j} s_i s_j p^*(s_i,s_j) = \left<s_i s_j \right> \forall i<j.$$
Using Lagrange multipliers $w = \{w_i, w_{ij} \}$, it can be shown that this distribution has the following parametrized form
$$p(s|w) = \frac{1}{Z(w)} \exp(\sum_i w_i s_i + \sum_{i<j} w_{ij} s_{i} s_{j})$$
Such a distribution is known as a "Boltzmann Machine" (BM). The parameters $w$ are called weights.

The maxent distribution $p^*(s)=p(s|w^*)$ has optimal weights $w^*$, which are such that the constraints are satisfied:
$$\left<s_i \right>_{w^*} = \left<s_i \right> ~\mbox{and}~\left<s_i s_j \right>_{w^*} = \left<s_i s_j \right>,$$
with notation $\left< s_i \right>_w \equiv \sum_{s_i} s_i p(s_i|w)$ etc.

These constraint equations can also be interpreted as the stationary equations $\nabla L(w) = 0$ of the loglikelihood $$L(w) = \sum_s p(s) \log p(s|w),$$ which can be written equivalently as $$L(w) = \sum_i w_i \left<s_i \right> + \sum_{i<j} w_{ij} \left<s_i s_j \right> - \log(Z(w)).$$ As a consequence, weights can be found by maximizing $L(w)$. E.g. gradient ascent leads to the well known Boltzmann Machine learning rule.
$$
\Delta w_i = \eta \frac{\partial}{\partial w_i} L(w) = \eta (\left<s_i \right>-\left<s_i \right>_w)
$$
and
$$
\Delta w_{ij} = \eta \frac{\partial}{\partial w_{ij}} L(w) = \eta( \left<s_i s_j \right>-\left<s_i s_j \right>_w)
$$
The optimal weights are unique because $L(w)$ is concave. To see this, consider the matrix of second derivatives of $\log Z(w)$. This matrix has the form of a covariance matrix, which is positive definite. Therefore $\log Z(w)$ is convex and $L(w)$ is concave. However, finding the optimal $w^*$ can still be demanding since the computation of e.g. $\left<s_i \right>_w$, $\left<s_i s_j \right>_w$ or $Z(w)$ require the summation over exponentially many states.

Sampling from a BM can be done by e.g. Gibbs sampling. The nice thing with Gibbs sampling in BM is that computing a probability of a single spin $s_i$ given the others $s_{\{j \neq i\}}$ is computationally cheap
$$
p(s_i | s_{\{j \neq i\}}) \propto \exp( (w_i + \sum_{j < i} w_{ji} s_j + \sum_{j>i} w_{ij} s_j ) s_i )
$$

An introductory page about BMs is e.g. http://www.scholarpedia.org/article/Boltzmann_machine

An interesting reference about the duality between statistics $\{\left<s_i \right>_w, \left<s_i s_j \right>_w\}$ and weights $\{w_i, w_{ij}\}$ is chapter 3 in
M. J. Wainwright and M. I. Jordan (2008). Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, Vol. 1, Numbers 1--2, pp. 1--305, December 2008

I hope this is of some use,

Wim