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.


1 Answer 1


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:

 [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:

x z y w 
1 0 1 1 

And for a sample of size 10:

  [,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
  • 1
    $\begingroup$ 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? $\endgroup$
    – qcc101
    Jan 9, 2019 at 19:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.