How to approximately sample an unknown non-parametric joint distribution given a complete set of partial conditional distributions? This question is related somewhat to Bayesian networks. In a BN, you have a DAG (directed acyclic graph). By supplying the root nodes with a sample, you can then follow the directed arcs to sample each of their children, and children's children, etc. until you have sampled the whole graph.
But what about a scenario where I have a complete, undirected (and thus cyclic) graph?  For example, if I have 4 discrete variables $a, b, c, d$, and I have the complete set of empirical non-parametric partial conditional distributions $C = \{ p(a|b), p(a|c), p(a|d), p(b|a), p(b|c), ... \}$, how could I best sample a random vector $V_n = \{a_n, b_n, c_n, d_n\}$ satisfying all partial conditional distributions $C$? Since there is no notion of a set of "starting nodes", the cycles become an issue.
Here are two ideas I have so far:


*

*Generate a vast number of random $V_n =  \{ p(a), p(b), p(c), p(d) \}$ and compute some approximation of likelihood for each like $\sim MLE(V_n) = \prod_{p(x|y)\in C, (x,y)\in V_n}{p(y)p(x|y)}$, and then stochastically select a winner by sampling over $\sim MLE(V_n)$. This seems very inefficient and unlikely to generalize well if the variable count and/or number of factors begins to increase substantially.

*A more complex stochastic method where you start with sample $V_{n_0} =  \{ p(a), p(b), p(c), p(d) \}$, and iteratively attempt to converge toward the actual joint distribution by randomly selecting and replacing one of the variables, lets say for example $c$, using only the conditionals $p(c|\ast)$ in $C$ to sample a new $c_1$, and thus creating a new sample vector $V_{n_1} =  \{ a_0, b_0, c_1, d_0 \}$. Then we could compare the likelihood estimates $\sim MLE(V_{n_0})$ and $\sim MLE(V_{n_1})$ and choose the more likely (deterministic or stochastic choice?) as the starting point for our next iteration. Hopefully this would converge after some reasonable number of iterations (even if the variable count and/or number of factors was high). However, how many iterations to take (or when to conditionally stop) is unclear to me.


Would either of my two ideas work at all? What better ways are there to efficiently generate quality samples approximating a non-parametric joint distribution given the complete set of partial conditional distributions? Maybe some other types of acceptance/rejection methods?
 A: The approximate joint distribution
If I understand the question correctly, the set of conditional distributions you have isn't enough to uniquely specify the full joint distribution. Instead, you want a way to sample from an approximate joint distribution that's consistent with the given conditional distributions. So, the first question is: what should this approximate joint distribution be? There are uncountably many valid possibilities, but the principle of maximum entropy (MaxEnt) argues for a particular choice: the distribution with maximum entropy that satisfies the constraints. This choice is consistent with all known information while otherwise imposing the minimum possible structure (i.e. it doesn't assume any information we don't actually have).
The MaxEnt distribution
Here's how to find the MaxEnt distribution consistent with the given conditional and marginal distributions. First, define the pairwise distributions $p(a,b), p(a,c), p(a,d), p(b,c), p(b,d), p(c,d)$. For example, $p(a,b) = p(a \mid b) p(b)$. These pairwise distributions completely capture the given marginal and conditional distributions.
The goal is then to find the joint distribution with maximum entropy, among those that are consistent with the pairwise distributions above. This is an optimization problem with linear constraints, and can be expressed as follows. There are 4 discrete variables, so the joint distribution can be expressed as a 4-dimensional table/array. Vectorize this table to represent the joint distribution as a vector $\hat{p}$. The optimization problem is then:
$$\max_\hat{p} \ -\sum_i \hat{p}_i \log \hat{p}_i
\quad \text{subject to:}$$
$$\hat{p}_i \ge 0
\quad \vec{1}^T \hat{p} = 1
\quad W \hat{p} = v$$
The objective function is the Shannon entropy (which assumes the definition $0 \log 0 = 0$). The first two constraints require that $\hat{p}$ is a valid probability distribution (i.e. it's non-negative and sums to one; $\vec{1}$ denotes a vector of ones). The last set of constraints enforce consistency with the pairwise distributions--marginalizing out each pair of variables from the joint distribution must yield the corresponding pairwise distribution. Each row of matrix $W$ contains a pattern of zeros and ones that specify how to sum over $\hat{p}$, and the corresponding entry of vector $v$ specifies the value this sum must take. For example, if $A$ and $B$ can take particular values $a'$ and $b'$, then $v$ must contain the element $p(a',b')$ and the corresponding row of $W$ must be set to sum $\hat{p}$ over all values where $A=a', B=b'$. This enforces the constraint $\sum_{c,d} \hat{p}(a',b',c,d) = p(a',b')$. There will be one constraint for every possible pair of values that each pair of variables can take.
This problem can be solved numerically using a standard constrained optimization solver. Then, the vector $\hat{p}$ can be reshaped back to a 4-dimensional probability table. Note that this approach won't scale well with the number of variables, as the size of $\hat{p}$ will increase multiplicatively.
Sampling from the MaxEnt distribution
Given a probability table representing the full joint distribution, sampling from it is trivial (e.g. using the inverse transform method).
