How to sample from joint Bernoulli distribution?

I have done searches on google.com and here but haven't found a related one. I apologize if this is a re-post.

Here is my question:

I have a set $n$ of Bernoulli random numbers, and I know the marginal distributions $p_i$ and pairwise joint distributions $p_{ij}$, I would like to generate $m$ i.i.d. samples from this random numbers. How should I do this? Any input is appreciated.

Maybe I should have asked this question in the first place. In practice, when talking about correlation, is the pair-wise joint distribution or pair-wise correlation more often seen?

Thanks!

-

First, knowing the pairwise distributions $p_{ij}$already gives you $p_i$. On the other hand, that data is not enough to determine the full join distribution.

If you only are interested in generating a Bernoulli that fit those $p_{ij}$, withouth regard for higher distributions, then simply assume a Markov process, and generate $x_1$ according to $p_1$, then $x_2$ according to $p_{2|1}=p_{12}/p_{1}$

Regarding you last question: no special preference. You can either give $p_1$ $p_2$ and $\rho_{12}$ or either $p_{12}$. In both cases there are three degrees of freedom, and it's trivial to convert from one representaion to another.

Update: The above only is useful to fit the pair distributions of consecutive elements, ($p_{1,2},p_{2,3},\dots$) but it does not fit other pair correlations. A complete treatment is more complicated.

A multivariate Bernolli is fully (and univocally) specified by the $2^n$ joint probabilities $p_{\bf x}$ restricted to $0 \le p_{\bf x} \le 1$ and $\sum p_{\bf x}=1$. An alternate (and perhaps more convenient) representation is given (refer to this answer) by the $2^n-1$ coefficients $a_{i...k}=E[x_i \cdots x_k]$, where $x_i=\pm1$ and $a_{\emptyset}=1$. These both representations are linearly related through a Hadamard matrix.

In our case, we are given the $n$ first order coefficients $a_i$ and the $n(n-1)/2$ second order coefficients $a_{i,j}$. If we are at liberty to choose the other coefficients, we can try by setting them to zero, and checking that the resulting ${\bf p}={\bf M a}$ falls into the allowed hypercube. Once you have find this (or some other acceptable solution), we have the full joint probability, and we can proceed to generate samples.

For samples generation, given the full joint probability, we can either generate the components in order using $P(x_i|x_{i-1} x_{i-2} \cdots x_1)$ or we can use some Gibss sampling or Metropolis algorithm

-
This does not guarantee the desired $p_{13}$. – Did Jan 19 '13 at 16:32
Hi Leonbloy, Thanks for your reply. Your comment "that data is not enough to determine the full join distribution" makes a good point. That's one motivation for me to study a distributionally robust model. I thought about the idea: first sample $i$ using $p_i$. But for $j$, should I use $p_j$ or $p_{j|i}$? and suppose I have sampled $i$ and $j$, then for $k$, I can use $p_{k|i}$ and $p_{k|j}$, which one should I use? this confused me. – Feng Jan 19 '13 at 16:39
@did: Your right, I'm only fitting the $n-1$ consecutive pairs $p_{i,i+1}$ instead of the full $n(n-1)/2$ pairs. Gotta think of it. – leonbloy Jan 19 '13 at 18:57
Hi Leonbloy, I apologize if I asked my questions before giving enough thinking to them. 1, why knowing all $p_{ij}$ would give you the information about $p_i$? 2, I'm kind of confused by the definition of $a_{i,...,k}$. Isn't it $p_{i,...,k}$. Why $x_i=\pm 1$? instead of 0 or 1? I might have seen this somewhere before, but couldn't remember where. Thanks! – Feng Jan 20 '13 at 19:07
@Feng 1) a marginal can be obtained by a joint denisty: $p_i= \sum_j p_{i,j}$. 2) the explanation is in the link, it might be useful to you or not. – leonbloy Jan 20 '13 at 19:11