# I want to simulate a random sample of length n from DAG of correlated Bernoulli's

Suppose I have a DAG of 4 vertices. Each vertex consists of a Bernoulli of parameter $$p$$.

It is the following: (Z) ---> (Y) (Z) ---> (W) (X) ---> (Y) ---> (W)

I hope it is clear. Anyway, I managed to find the joint density function which is the following: $$f(x,y,z,w) = f_x(x) * f_z(z) * f_y|(x,z)(y | x,z) * f_w|(y,z)(w|y,z)$$ Now my problem is that I want to find a random sample of length n from this graph. I am struggling to see what it means to have a sample of length $$n$$: do I have to find $$n$$ values of the density function in random points?

Anyway, do you know how I can do this in R? Especially the conditional part.

You seem to be confused, so maybe you should look first at some simpler situations. I will interpret random sample of size $$n$$ to mean to simulate (or "draw") $$n$$ (some positive integer) realization of some random variable. So, if that random variable (rv) is Bernoulli with some probability $$p$$, a realization must be one of the possible values 0 or 1. A random sample of size 10 in R:

 rbinom(10,1,0.75)
[1] 1 1 1 0 1 1 1 0 1 1


Here the probability of the value 1 was choosen as 0.75.

A sample from your DAG will be a vector of length 4, which in your case will be all 0 or 1's. A random sample of size 10 will be a collection of 10 such 4-vectors, which could conveniently be given as a $$10\times 4$$-matrix.

To draw one such 4-vector, you must use the structure of the DAG, as given in your factored probability mass function (pmf). In your formula $$f(x,y,z,w) = f_x(x) * f_z(z) * f_y|(x,z)(y | x,z) * f_w|(y,z)(w|y,z)$$
there is no conditioning in the first two factors, so you start by drawing $$X,Z$$ (you must specify two probabilities first). The third factor conditions on $$X,Z$$, so can be drawn now (but you need to specify a probability parameter that depends on $$x,z$$). Finally you can draw $$W$$ from the last factor in the same way. An example where I use as full specification $$\DeclareMathOperator{\P}{\mathbb{P}} \P(X=1)=p=0.75,\quad \P(Z=1)=q=0.5,\quad \P(Y=1\mid X=x,Z=z)=p_y(x,z),\quad \P(W=1\mid Y=y,Z=z)=p_w(y,z)$$ and then I will use for example: $$p_y(x,z)= 0.1+\frac{x+z}{3},\\ p_w(y,z)= 0.2+\frac{y+z}{3}$$ R code for one draw:

draw_1_from_DAG <- function() {
x  <-  rbinom(1,1,0.75)
z  <-  rbinom(1,1,0.5)
y  <-  rbinom(1,1,0.1+(x+z)/3)
w  <-  rbinom(1,1,0.2+(y+z)/3)
c(x=x,z=z,y=y, w=w)
}


One example call:

  draw_1_from_DAG()
x z y w
1 0 1 1


And for a sample of size 10:

  replicate(10,draw_1_from_DAG())
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
x    1    0    1    0    1    1    1    1    1     1
z    0    1    0    1    1    0    1    0    1     1
y    0    1    1    1    1    0    1    1    0     1
w    0    1    1    1    1    0    0    0    1     1

• Wow, this is the clearest answer I have ever seen. Thank you very much. I just have one last doubt: Is there a "standard" way to define $p_y(x,z)$ and $p_w(y,z)$? Like maybe using Bayes theorem to find the probability of the conditional event? – qcc101 Jan 9 at 19:39