How to simulate a directed acyclic graph?

My goal is to have a code to identify and estimate the regression models associated with a directed acyclic graph (DAG), which is obtained through the reduction of a given weighted, directed graph G. The only input for this operation is the graph G, and the only output expected is a set of estimated regression models that simulate the edge weights of G. Here is the 4-step procedure (working fine) to achieve this goal. I like to have an efficient code using the 4-step operations described below to automate the simulation method.

STEP 1: Construct a directed graph:

el = {"stel" $DirectedEdge] "sfin", "scst" \[DirectedEdge] "stel", "scst" \[DirectedEdge] "sfin", "scst" \[DirectedEdge] "sbus", "scst" \[DirectedEdge] "swhl", "scst" \[DirectedEdge] "ma6", "sfin" \[DirectedEdge] "stel", "sfin" \[DirectedEdge] "sbus", "sfin" \[DirectedEdge] "swhl", "sbus" \[DirectedEdge] "stel", "sbus" \[DirectedEdge] "sfin", "sbus" \[DirectedEdge] "ma6", "swhl" \[DirectedEdge] "scst", "swhl" \[DirectedEdge] "sbus", "swhl" \[DirectedEdge] "ma6", "ma6" \[DirectedEdge] "swhl"}; ew = {0.0299, 0.0372, 0.0223, 0.0424, 0.0445, 0.0221, 0.0215, 0.0590, 0.0533, 0.0441, 0.0438, 0.0221, 0.0987, 0.0394, 0.0208, 0.0317}; wUpstream = Graph[VertexList[el], el, EdgeWeight -> ew, VertexLabels -> Placed["Name", Center], VertexSize -> Medium]  STEP 2: Construct a directed acyclic graph (DAG) using the above graph: ClearAll[wagNew, EL, dag]; EL[dir_ : Up] := g |-> EdgeList[ g, _?((dir /. {Up -> Less, Down -> Greater}) @@ PropertyValue[{g, VertexList[{#}]}, VertexCoordinates][[All, 2]] &)]; (* by @kglr *) wagNew = Graph[ EdgeList[wUpstream], EdgeWeight -> PropertyValue[wUpstream, EdgeWeight], EdgeShapeFunction -> ({Arrowheads[{{.035, .75}}], Arrow[#]} &), DirectedEdges -> True, options, GraphLayout -> "LayeredDigraphEmbedding", ImagePadding -> 10, ImageSize -> 450 ]; dag = HighlightGraph[wagNew, Graph@EL[Down]@wagNew, GraphHighlightStyle -> "DehighlightHide", EdgeLabels -> "EdgeWeight", ImagePadding -> 20]  STEP 3: Design a data-generating mechanism based on the DAG above. Here is the procedure to generate random data sets for simulations: (1) start with layer 1 in the graph dag to produce simulation data for scst first as it is placed at the top of the layered graph; (2) then move to layer 2, which includes sfin receiving input from scst only. Hence, the data generated for scst should be used in the generation of data for sfin; (3) then move to layer 3, including swhl only, receiving inputs from sfin and scst. The data generation for swhl should use the data already generated for sfin and scst. This process goes on layer by layer to fully cover the data generation for all vertices in the entire network. Note: in-degree edges for vertex i determines the set of data to be generated for vertex I. The following code simply applies the procedure described. ClearAll[vl,nDegree,Scst, Stel, Sfin, Sbus, Swhl, Ma6]; vl = VertexList[Graph@EL[Down]@wagNew]; nDegree = Thread[vl -> Normalize[VertexDegree[Graph@EL[Down]@wagNew], Total]] // N SeedRandom[203]; n = 10^6; brD[a_] := BernoulliDistribution[a]; pScst = "scst" /. nDegree; sampScst = RandomVariate[brD[pScst], n]; pSfin = ("sfin" + 0.0223 sampScst) /. nDegree; sampSfin = Table[RandomVariate[brD[pSfin[[i]]], 1], {i, n}] // Flatten; pSwhl = ("swhl" + 0.0445 sampScst + 0.0533 sampSfin) /. nDegree; sampSwhl = Table[RandomVariate[brD[pSwhl[[i]]], 1], {i, n}] // Flatten; pSbus = ("sbus" + 0.0424 sampScst + 0.0590 sampSfin + 0.0396 sampSwhl) /. nDegree ; sampSbus = Table[RandomVariate[brD[pSbus[[i]]], 1], {i, n}] // Flatten; pStel = ("stel" + 0.0372 sampScst + 0.0215 sampSfin + 0.0441 sampSbus) /. nDegree; sampStel = Table[RandomVariate[brD[pStel[[i]]], 1], {i, n}] // Flatten; pMa6 = ("ma6" + 0.0221 sampScst + 0.0208 sampSwhl + 0.0221 sampSbus) /. nDegree; sampMa6 = Table[RandomVariate[brD[pMa6[[i]]], 1], {i, n}] // Flatten; means = Thread[ vl -> {Mean[sampScst], Mean[sampStel], Mean[sampSfin], Mean[sampSbus], Mean[sampSwhl], Mean[sampMa6]}] // N  STEP 4: Construct linear regression models for the simulation of the DAG. The construction of the regression models below is implemented by following the exact same procedure described in Step 3. Regression coefficients of the following models simulate the edge weights of the DAG. TableForm\[ Transpose@{ {"Sfin=f(Scst)", "Swhl=f(Scst, Sfin)", "Sbus=f(Scst, Sfin, Swhl)", "Stel=f(Scst, Sfin, Sbus)", "Ma6=f(Scst, Swhl, Sbus)"}, {LinearModelFit\[Transpose@{sampScst, sampSfin}, {1, Scst}, Scst$,
LinearModelFit$Transpose@{sampScst, sampSfin, sampSwhl}, {1, Scst, Sfin}, {Scst, Sfin}$,
LinearModelFit$Transpose@{sampScst, sampSfin, sampSwhl, sampSbus}, {1, Scst, Sfin, Swhl}, {Scst, Sfin, Swhl}$,
LinearModelFit$Transpose@{sampScst, sampSfin, sampSbus, sampStel}, {1, Scst, Sfin, Sbus}, {Scst, Sfin, Sbus}$,
LinearModelFit$Transpose@{sampScst, sampSwhl, sampSbus, sampMa6}, {1, Scst, Swhl, Sbus}, {Scst, Swhl, Sbus}$}
}]


The code given works as expected and all the coefficients of the regressions simulate the edgewights of the DAG.

My interest is to automate this 4-step simulation method for large graphs. A good start would be to construct a function, such as simDAG[lg_, df_, rf_]:=... to perform the simulation at one-go. For this, the inputs constructed at each step will be used as arguments of the simDAG[...]:=...:

1. lg denotes a layered, weighted, directed acyclic graph;
2. df denotes a set of distribution functions, each one of which is associated with a vertex in the DAG. In the above example, for sake of simplicity, I used Bernoulli Distribution for all vertices.
3. rf denotes regression functions. In the example, I construct the regression functions based on the layered structure of the DAG. For example, vertex sfin = f(scst) or vertex swhl = f(scst,sfin) or vertex ma6 = f(sbus,swhl, scst), etc. Namely, the in-degree edges of a vertex "swhl" are used as arguments of the function swhl.

EDIT 1

The question in this edit was not in the original request but it appears that the simulation method will be completed with an additional code aimed to automate single variable regression models associated with each edge in the DAG.

Given a directed edge from vertex i to vertex j (i.e.g i-->j), construct a regression model of the form:

vertex j = f(vertex i),

where f(.) is a linear function similar to the ones in the original question. These regression models should be constructed for all of the edges in the DAG given. The estimations of these models should use the already simulated data in Step 3. For example,

If i denotes scst and j denotes sfin, then we can construct a regression model of the form: sampSfin = f(sampScst). Likewise, construct the regression models for all the edges:

• sampSwhl = f(sampSfin)
• sampSwhl = f(sampScst)
• sampStel = f(sampScst), and so on.

Thanks.

• There is a lot of code here, and very little explanation. No wonder that no one commented. Can you explain what you are looking for in clear and concise plain language? Commented Nov 13, 2023 at 18:28
• @Szabolcs: Thanks for your interest in the question. I will explain my objective in a plain language by editing my question. Commented Nov 13, 2023 at 18:59

An idea in order to generalize the steps 3 and 4, considering the definition of the weighted graph given by edges and weights:

edges = DirectedEdge @@@ {
"stel"->"sfin", "scst"->"stel", "scst"->"sfin", "scst"->"sbus", "scst"->"swhl",
"scst"->"ma6", "sfin"->"stel", "sfin"->"sbus", "sfin"->"swhl", "sbus"->"stel",
"sbus"->"sfin", "sbus"->"ma6", "swhl"->"scst", "swhl"->"sbus", "swhl"->"ma6", "ma6"->"swhl"};
weights = {0.0299, 0.0372, 0.0223, 0.0424, 0.0445, 0.0221, 0.0215,
0.0590, 0.0533, 0.0441, 0.0438, 0.0221, 0.0987, 0.0394, 0.0208, 0.0317};


Borrowing the definitions from step 1 & 2 in order to generate the directed acyclic graph (dag).

wag = Graph[edges, EdgeWeight -> weights, VertexLabels -> Placed["Name", Center]
, VertexSize -> Large, GraphLayout -> "LayeredDigraphEmbedding"];

edgeFilter[dir_ : Up] := g |-> EdgeList[g, _?((dir /. {Up -> Less, Down -> Greater}) @@
PropertyValue[{g, VertexList[{#}]}, VertexCoordinates][[All,2]] &)];(* kglr *)

ew = KeyTake[AssociationThread[EdgeList@wag -> AnnotationValue[wag, EdgeWeight]], edgeFilter[Down]@wag];

dag = Graph[Keys@ew, EdgeWeight -> Values@ew
, EdgeLabels -> "EdgeWeight"
, VertexLabels -> Placed["Name", Center]
, VertexLabelStyle -> {Directive[14, White, Bold]}
, VertexStyle -> Red
, VertexSize -> .5
, ImageSize -> Large]


Automating data generating mechanism:

ClearAll[ sample, regression, distribution, rv, aggregation, degree, var, f];
n = 10^6;
degree = AssociationThread[VertexList@dag -> Normalize[VertexDegree@dag, Total]] // N;
distribution = <|# -> BernoulliDistribution & /@ VertexList@dag|>;

aggregation[vertex_, func_ : sample] := Module[{},
degree[vertex] +
, Rule[e : DirectedEdge[u_, d_], w_] /; d == vertex :>
Times @@ {func@u, w }]];

# -> aggregation[#, sample] & /@ TopologicalSort@dag
(* Output *)

(*
{"scst" -> 0.208333, "sfin" -> 0.166667 + 0.0223 sample["scst"],
"swhl" -> 0.166667 + 0.0445 sample["scst"] + 0.0533 sample["sfin"],
"sbus" -> 0.208333 + 0.0424 sample["scst"] + 0.059 sample["sfin"] + 0.0394 sample["swhl"],
"ma6" -> 0.125 + 0.0221 sample["sbus"] + 0.0221 sample["scst"] + 0.0208 sample["swhl"],
"stel" -> 0.125 + 0.0441 sample["sbus"] + 0.0372 sample["scst"] + 0.0215 sample["sfin"]}
*)

(* helper functions for RandomVariate calculation *)

rv[vertex_][num_?NumericQ] := RandomVariate[distribution[vertex][num], n];
rv[vertex_][expr_] := Hold[Flatten[RandomVariate[distribution[vertex]@#, 1] & /@ expr]];

BlockRandom[
ts = Association[# -> aggregation[ #, sample] & /@ TopologicalSort@dag];
KeyValueMap[Set[sample@#, ReleaseHold@rv[#][#2]] &, ts];
, RandomSeeding -> 203
];

means = Thread[# -> Mean@*sample /@ #] &[VertexList@dag] // N
(* Output *)
(*
{"scst" -> 0.208025, "stel" -> 0.146883, "sfin" -> 0.17183,
"sbus" -> 0.234246, "swhl" -> 0.185425, "ma6" -> 0.13858}
*)


Automating the linear regression functions generation:

regression[vertex_, func_] := Module[{},
LinearModelFit[
Transpose@Join[ #@func, {func@vertex}]
, {1, Splice@#@var}, #@var, WorkingPrecision -> 6] &
[Cases[EdgeList@dag, DirectedEdge[u_, d_] /; d == vertex :> #@u] &]
];

Grid[
{# -> f @@ Cases[VertexInComponent[dag, #,  1], v_ /; v != #], regression[#, sample]} & /@
Cases[TopologicalSort@dag, v_ /; Length@VertexInComponent[dag, v,  1] > 1]
, Alignment -> Left]


UPDATE (related to EDIT 1 from OP)

Sort@Cases[EdgeList@dag, u_ \[DirectedEdge] d_ :>
{d -> f[u]
, LinearModelFit[
Transpose@{sample[u], sample[d]}
, {1, var[u]}, var[u], WorkingPrecision -> 6]}]  // TableForm


• Thanks very much for the code. It produces the expected output. The code runs very slow, though. Furthermore, can you tell me how to run simple regressions between any two vertices (i.e., single edge) separately in order to re-establish the DAG edge-by-edge using the data generated. With your code, I am having difficulty to design a regression model with a single regressor of interest. Commented Nov 19, 2023 at 23:06
• @TugrulTemel Thank you for your feedback! Regarding performance, the timings should be comparable with the ones from the initial version and they are strongly related to the way the random data sets are generated in step 3. Regarding your other question, if you could add a short example (similar with those in your original post) I might be able to try to generalize it as well. Commented Nov 20, 2023 at 7:56
• @TugrulTemel Please see the updated answer related to EDIT 1. Commented Nov 21, 2023 at 20:48
• Thank you very much for your time and effort in putting the code together. Excellent... Commented Nov 21, 2023 at 21:03