wolffd@0: function dag = sample_dag(P)
wolffd@0: % SAMPLE_DAG Create a random directed acyclic graph with edge probabilities P(i,j)
wolffd@0: % dag = sample_dag(P)
wolffd@0: %
wolffd@0: % This uses rejection sampling to reject graphs with directed cycles.
wolffd@0: 
wolffd@0: done = 0;
wolffd@0: directed = 1;
wolffd@0: iter = 1;
wolffd@0: while ~done
wolffd@0:   dag = binornd(1, P); % each edge is an indep Bernoulli (0/1) random variable
wolffd@0:   dag = setdiag(dag, 0);
wolffd@0:   done = acyclic(dag, directed);
wolffd@0:   iter = iter + 1
wolffd@0: end