In an application scenario, I have two variables that I have to assume are mutually influencing each other. Let’s say the model could look like this, with all variables binary and all variables except for B observed:

```
A ---> (B) <---> D
| |
v |
C <----- +
```

How would I go about modeling the bidirectional edge between B and D in a (num)pyro model?

Ideally, I would like to do something like

```
...
B = sample("B", dist.Bernoulli(logits=k1 + k2*A + k3*D))
D = sample("D", dist.Bernoulli(logits=k4 + k5*B), obs=...)
...
```

but this obviously does not work because D is not yet defined in the first line. (The graph is loopy.)

I guess one approach would be to merge B and D into a single 2D binary variable and model their joint distribution instead? That’s a) a little “ugly” to me personally but ok, but more importantly b) how would I then handle the fact that D is observed but B is not?

I could also imagine that I could do something using a `factor`

, but I have never worked with that and am a bit unsure about how I would have to apply that here.

(I do realize that this is at least as much a general modeling question as it is a pyro-specific question…)