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. 
 A: 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

