Drawing random variates from a partially described probability distribution

Question

I have a probability distribution over $\{0,1\}^n$ but instead of knowing the full joint distribution $p(x_1,\dots,x_n)$, I only know $p(x_i=x_j)$ for each $i,j$. How could I draw a random binary vector $x$ from some distribution that has these marginals? Since $n$ is large I would rather not search over all distributions (in the $(2^n-1)$-dimensional simplex) to find one that matches (this can be hard even though by Caratheodory's theorem there is such a distribution with support only $\binom n 2$, because what's in the support is not known). Hopefully there's a more efficient method if all I want to do is draw a random variate.

If it helps, assume that $p(x_i=0)=p(x_i=1)=0.5$ for all $i$.

Thanks!

Answer 1

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

Answer 2

It would help if you know the marginal of at least one variable $x_i$. Suppose this is the case, i.e. without loss of generality suppose that you know $p(x_1=0)$.

Then $p(x_1=i, x_2=j) = \left\{ \begin{array}{ll} p(x_1 = i) p(x_1=x_2) \quad &\mbox{if} \ i = j \\ p(x_1 = i) (1 - p(x_1 = x_2)) \quad & \mbox{if} \ i \neq j. \end{array} \right.$

Inductively, suppose we know $p(x_1 = i_1, ... , x_n = i_n)$. Then \begin{equation*} p(x_1 = i_1, ... x_n = i_n, x_{n+1} = i_{n+1}) = \left\{ \begin{array}{ll} p(x_1 = i_1,...,x_n=i_n) p(x_n=x_{n+1}) \quad &\mbox{if} \ i_n = i_{n+1} \\ p(x_1 = i_1,...,x_n=i_n) (1 - p(x_n=x_{n+1})) \quad & \mbox{if} \ i_n \neq i_{n+1}. \end{array} \right.\end{equation*} So this procedure gives the full distribution of $x_1,...,x_n$.

From this construction, you immediately notice that the family of distributions is fixed if you only have the marginal distribution of $x_1$ and all probabilities $p(x_n = x_{n+1})$. In particular, this raises the question what condition you require on your $p_{ij} = p(x_i = x_j)$, so that these are consistent.

This idea may be a start for a general solution, i.e. without the mentioned assumption.

Remark: there is an assumption of independence in my construction, which may not hold in your case. ($\{x_i=x_{i+1}\}$ is assumed to be independent of $x_1,...x_i$).