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