Repository: pgmpy/pgmpy

Issue:
Missing `get_random` for `LinearGaussianCPD`
We have `get_random` for `TabularCPD`
but not for `LinearGaussianCPD`.

We do have `get_random_cpds` in BaysianNetwork and LinearGaussianBaysianNetwork.

With a small refactoring, we could expose `get_random` for `LinearGaussianCPD`.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) Report
Attention: Patch coverage is `25.00000%` with `9 lines` in your changes missing coverage. Please review.
> Project coverage is 34.56%. Comparing base [(`720b0ae`)](https://app.codecov.io/gh/pgmpy/pgmpy/commit/720b0ae4fd57b4f8a9ebeb076ea860e799f4d39d?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) to head [(`dd2c625`)](https://app.codecov.io/gh/pgmpy/pgmpy/commit/dd2c6253154f29bbac369696b24ceb554ef87556?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy).
> Report is 15 commits behind head on dev.

| [Files with missing lines](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) | Patch % | Lines |
|---|---|---|
| [pgmpy/factors/continuous/LinearGaussianCPD.py](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&filepath=pgmpy%2Ffactors%2Fcontinuous%2FLinearGaussianCPD.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy#diff-cGdtcHkvZmFjdG9ycy9jb250aW51b3VzL0xpbmVhckdhdXNzaWFuQ1BELnB5) | 28.57% | [5 Missing :warning: ](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) |
| [pgmpy/models/LinearGaussianBayesianNetwork.py](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&filepath=pgmpy%2Fmodels%2FLinearGaussianBayesianNetwork.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy#diff-cGdtcHkvbW9kZWxzL0xpbmVhckdhdXNzaWFuQmF5ZXNpYW5OZXR3b3JrLnB5) | 0.00% | [2 Missing :warning: ](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) |
| [...actors/test\_continuous/test\_Linear\_Gaussian\_CPD.py](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&filepath=pgmpy%2Ftests%2Ftest_factors%2Ftest_continuous%2Ftest_Linear_Gaussian_CPD.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy#diff-cGdtcHkvdGVzdHMvdGVzdF9mYWN0b3JzL3Rlc3RfY29udGludW91cy90ZXN0X0xpbmVhcl9HYXVzc2lhbl9DUEQucHk=) | 33.33% | [2 Missing :warning: ](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy) |

> :exclamation:  There is a different number of reports uploaded between BASE (720b0ae) and HEAD (dd2c625). Click for more details.
> 
> <details><summary>HEAD has 17 uploads less than BASE</summary>
>
>| Flag | BASE (720b0ae) | HEAD (dd2c625) |
>|------|------|------|
>||18|1|
></details>

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##              dev    #1931       +/-   ##
===========================================
- Coverage   93.61%   34.56%   -59.06%     
===========================================
  Files         143      144        +1     
  Lines       18455    18692      +237     
===========================================
- Hits        17277     6461    -10816     
- Misses       1178    12231    +11053     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pgmpy/pgmpy/pull/1931?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pgmpy).


Patch:
diff --git a/examples/How to define CPDs.ipynb b/examples/How to define CPDs.ipynb
new file mode 100644
index 00000000..023627ee
--- /dev/null
+++ b/examples/How to define CPDs.ipynb
@@ -0,0 +1,318 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# How to define TabularCPD and LinearGaussianCPD\n",
+    "\n",
+    "One can implement `TabularCPD` for discrete variables inside `BayesianNetwork`.\n",
+    "and `LinearGaussianCPD` for continuous variables inside `LinearGaussianBayesianNetwork`. \n",
+    "\n",
+    "In this tutorial, we will demonstrate how to define each CPD.\n",
+    "\n",
+    "## TabularCPD for discrete variables\n",
+    "\n",
+    "In tabular CPD, the probability for discrete variable is given as a table. Let us start with examples for independent variables."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pgmpy.factors.discrete import TabularCPD\n",
+    "\n",
+    "cpd_coin_fair = TabularCPD(variable=\"coin\", variable_card=2, values=[[0.5], [0.5]], state_names={'Coin': ['Head', 'Tail']})\n",
+    "cpd_coin_biased = TabularCPD(variable=\"coin\", variable_card=2, values=[[0.9], [0.1]], state_names={'Coin': ['Head', 'Tail']})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For a coin flip, we have discrete number (2) possible outcomes, therefore `variable_card=2` is passed. `values` pass probabilities for each outcome, which sum up to 1. `state_names` is optional, it gives names for each outcome."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cpd_smoke = TabularCPD(variable=\"Smoker\", variable_card=2, values=[[0.3], [0.7]], state_names={'Smoker': ['Non-smoker', 'Smoker']})\n",
+    "cpd_pollution = TabularCPD(variable=\"Pollution\", variable_card=3, values=[[0.7], [0.29], [0.01]], state_names={'Pollution': ['Clean', 'Bad', 'Fatal']})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `Pollution` and `Smoker` variables do not depend on other variables. They take categorical values. `variable_card` denotes how many categories the variable can take. For example `Smoker` can take a binary value since `variable_card=2`. `values` is an array of probability values for each category. Note that the probability `values` sum up to 1.\n",
+    "\n",
+    "### Tabular CPD with multiple evidence.\n",
+    "\n",
+    "We next consider another discrete variable `Cancer`, with a model assumption that it depends on `Pollution` and `Smoker`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cpd_cancer = TabularCPD(\n",
+    "    variable=\"Cancer\",\n",
+    "    variable_card=2,\n",
+    "    values=[[0.20, 0.15, 0.03, 0.05, 0.001, 0.02],\n",
+    "            [0.80, 0.85, 0.97, 0.95, 0.999, 0.98]],\n",
+    "    evidence=[\"Smoker\", \"Pollution\"],\n",
+    "    evidence_card=[2, 3],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For `Cancer` variable, we pass another `TabularCPD`. `evidence_card` denotes cardinality of each evidence variable. There are total `2*3=6` different combinations of evidence. Note the values of `evidence` affects probability of `Cancer`. We have 6 columns in `values` denoting conditional probabilities. `values` sums up to 1 columnwise. The columns are partitioned by the first evidence first. The first three columns correspond to non-smokers.\n",
+    "\n",
+    "We next consider another discrete variable `D`, with a model assumption that it depends on `A`, `B` and `C`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cpd = TabularCPD(\n",
+    "    variable=\"D\",\n",
+    "    variable_card=2,\n",
+    "    values=[[0.20, 0.15, 0.93, 0.05, 0.001, 0.02, 0.10, 0.25 ],\n",
+    "            [0.80, 0.85, 0.07, 0.95, 0.999, 0.98, 0.90, 0.75 ]],\n",
+    "    evidence=[\"A\", \"B\", \"C\"],\n",
+    "    evidence_card=[2, 2, 2],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+------+------+------+------+------+-------+------+------+------+\n",
+      "| A    | A(0) | A(0) | A(0) | A(0) | A(1)  | A(1) | A(1) | A(1) |\n",
+      "+------+------+------+------+------+-------+------+------+------+\n",
+      "| B    | B(0) | B(0) | B(1) | B(1) | B(0)  | B(0) | B(1) | B(1) |\n",
+      "+------+------+------+------+------+-------+------+------+------+\n",
+      "| C    | C(0) | C(1) | C(0) | C(1) | C(0)  | C(1) | C(0) | C(1) |\n",
+      "+------+------+------+------+------+-------+------+------+------+\n",
+      "| D(0) | 0.2  | 0.15 | 0.93 | 0.05 | 0.001 | 0.02 | 0.1  | 0.25 |\n",
+      "+------+------+------+------+------+-------+------+------+------+\n",
+      "| D(1) | 0.8  | 0.85 | 0.07 | 0.95 | 0.999 | 0.98 | 0.9  | 0.75 |\n",
+      "+------+------+------+------+------+-------+------+------+------+\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(cpd)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each column of `values` correspond to a different combinations of `evidence` variables. If the evidence variables were lexicographic as above `(A, B, C)`, then the columns will be also lexicographic. \n",
+    "\n",
+    "```(A(0)B(0)C(0), A(0)B(0)C(1), A(0)B(1)C(0), A(0)B(1)C(1), ...)```\n",
+    "\n",
+    "It is as though we loop over A, and B, and C, in a nested order. \n",
+    "\n",
+    "### CPD with random values\n",
+    "It is also possible to create a random tabular CPD."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.55203474],\n",
+       "       [0.44796526]])"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "cpd_coin_random = TabularCPD.get_random(variable=\"coin\", cardinality = {'coin':2}, state_names={'Coin': ['Head', 'Tail']})\n",
+    "cpd_coin_random.get_values()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note the probabilities can be retrieved by `get_values` and they add up to 1 again. If `cardinality` is missing, the discrete variable is assumed to be binary. \n",
+    "\n",
+    "When there are `evidence`, the cardinality of evidence needs to be also specified."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.30404579, 0.99344749, 0.47933983, 0.58564998, 0.54269003,\n",
+       "        0.28164646],\n",
+       "       [0.69595421, 0.00655251, 0.52066017, 0.41435002, 0.45730997,\n",
+       "        0.71835354]])"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "cpd_cancer_random = TabularCPD.get_random(\n",
+    "    variable=\"Cancer\",\n",
+    "    cardinality = {\"Cancer\": 2, \"Smoker\": 2, \"Pollution\": 3},\n",
+    "    evidence=[\"Smoker\", \"Pollution\"],\n",
+    ")\n",
+    "cpd_cancer_random.get_values()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([1., 1., 1., 1., 1., 1.])"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "cpd_cancer_random.get_values().sum(axis = 0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note again the probability table sums up to 1 for each column.\n",
+    "\n",
+    "## LinearGaussianCPD for continuous variables\n",
+    "\n",
+    "We define CPDs for each variable. `LinearGaussianCPD` assumes that each variable takes normal distribution. \n",
+    "Assume `Healthy` and `Wealthy` variables have no parents (empty evidence). "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<LinearGaussianCPD: P(Healthy) = N(4; 2) at 0x15735ed20\n",
+      "<LinearGaussianCPD: P(Wealthy) = N(2; 3) at 0x1566c2930\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Step 2: Define the CPDs.\n",
+    "from pgmpy.factors.continuous import LinearGaussianCPD\n",
+    "cpd_healthy = LinearGaussianCPD(variable=\"Healthy\", beta=[4], std=2, evidence=[])\n",
+    "cpd_wealthy = LinearGaussianCPD(variable=\"Wealthy\", beta=[2], std=3, evidence=[])\n",
+    "\n",
+    "import pprint\n",
+    "pprint.pp(cpd_healthy)\n",
+    "pprint.pp(cpd_wealthy)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Above we defined, Healthy as normal distribution with mean 4 and standard deviation 2, and wealthy with mean 2, and standard deviation 3. \n",
+    "(Assume bigger variation in people's wealth than health, bigger mean for health than wealth)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<LinearGaussianCPD: P(Happy | Healthy, Wealthy) = N(3*Healthy + 2*Wealthy + 1; 5) at 0x154160050\n"
+     ]
+    }
+   ],
+   "source": [
+    "cpd_happy = LinearGaussianCPD(\n",
+    "    variable=\"Happy\",\n",
+    "    beta=[1, 3, 2],\n",
+    "    std=5,\n",
+    "    evidence=[\"Healthy\", \"Wealthy\"],\n",
+    ")\n",
+    "pprint.pp(cpd_happy)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `Happy` variable has mean `3*Healthy + 2*Wealthy + 1`, this formula is determined by passing `beta` and `evidence` variables. Note that the first element of `beta` is the intercept (constant term). The rest of the elements in `beta` each match the evidence.\n",
+    "The standard deviation of normal distribution of `Happy` is set by `std`."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pgmpy",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pgmpy/estimators/ExpertKnowledge.py b/pgmpy/estimators/ExpertKnowledge.py
index 9f51999f..190834e0 100644
--- a/pgmpy/estimators/ExpertKnowledge.py
+++ b/pgmpy/estimators/ExpertKnowledge.py
@@ -1,3 +1,5 @@
+from itertools import chain
+
 from pgmpy.global_vars import logger
 
 
@@ -22,18 +24,22 @@ class ExpertKnowledge:
             graph structure. Refer to the algorithm documentation for details
             on how the argument is handled.
 
-    temporal order: list of lists (default: None)
+    temporal order: iterator (default: None)
             The temporal ordering of variables according to prior knowledge.
-            Each list in the list of lists contains variables with the same
-            temporal significance; the more prior (parental) variables (list) are at
-            the start while the priority decreases as we go down the list.
+            Each list/structure in the (2 dimensional) iterator contains
+            variables with the same temporal significance; the more prior
+            (parental) variables are at the start while the priority decreases
+            as we go move towards the end of the structure (iterator).
 
     Examples
     --------
     Import an example model from pgmpy.utils
 
     >>> from pgmpy.utils import get_example_model
+    >>> from pgmpy.estimators import ExpertKnowledge, PC
+    >>> from pgmpy.sampling import BayesianModelSampling
     >>> asia_model = get_example_model("asia")
+    >>> cancer_model = get_example_model("cancer")
 
     **Required and forbidden edges**
 
@@ -43,7 +49,6 @@ class ExpertKnowledge:
 
     **Use during structure learning**
 
-    >>> from pgmpy.estimators import PC
     >>> data = BayesianModelSampling(asia_model).forward_sample(size=int(1e4))
     >>> est = PC(data)
     >>> est.estimate(
@@ -51,6 +56,20 @@ class ExpertKnowledge:
     ...         expert_knowledge=expert_knowledge,
     ...         show_progress=False,
     ...     )
+
+    **Temporal order**
+
+    >>> expert_knowledge = ExpertKnowledge(temporal_order=[["Pollution", "Smoker"], ["Cancer"], ["Dyspnoea", "Xray"]])
+
+    **Use during structure learning**
+
+    >>> data = BayesianModelSampling(cancer_model).forward_sample(size=int(1e4))
+    >>> est = PC(data)
+    >>> est.estimate(
+    ...         variant="stable",
+    ...         expert_knowledge=expert_knowledge,
+    ...         show_progress=False,
+    ...     )
     """
 
     def _validate_edges(self, edge_list):
@@ -63,6 +82,86 @@ class ExpertKnowledge:
         else:
             return edge_list
 
+    def _validate_temporal_order(self, nodes):
+        """
+        Method to check consistency of temporal order with nodes of a graph.
+
+        The temporal order, if specified by the user, is currently used by the PC
+        algorithm. The temporal order of all nodes present in the graph/dataset
+        need to be present in the ExpertKnowledge instance.
+
+        Parameters
+        ----------
+        nodes: iterable
+            A collection of nodes present in a dataset/graph object.
+        """
+        if self.temporal_order == [[]]:
+            return
+
+        # Check if no node is present in multiple tiers
+        if len(set.intersection(*map(set, self.temporal_order))) != 0:
+            raise ValueError("Node found in multiple tiers of temporal order.")
+
+        # Check if all nodes are present in the temporal order
+        if set(chain(*self.temporal_order)) != set(nodes):
+            raise ValueError(
+                f"Missing nodes in temporal order - {set(nodes) - set(chain(*self.temporal_order))}"
+            )
+
+    def _get_temporal_ordering(self, temporal_order):
+        """
+        Method to check consistency of temporal order with nodes of a graph.
+
+        The temporal order, if specified by the user, is currently used by the PC
+        algorithm. The temporal order of all nodes present in the graph/dataset
+        need to be present in the ExpertKnowledge instance.
+
+        Parameters
+        ----------
+        temporal_order: iterator
+            The temporal ordering of variables according to prior knowledge.
+
+        Returns
+        --------
+        temporal_ordering: dict
+            Dictionary with the tier (0, 1, 2, 3 etc.) for each node.
+        """
+        if not hasattr(temporal_order, "__iter__"):
+            raise TypeError(
+                f"Expected iterator type for temporal order. Got {type(temporal_order)} instead."
+            )
+
+        temporal_ordering = dict()
+        for order, tier in enumerate(self.temporal_order):
+            for node in tier:
+                if node in temporal_ordering:
+                    raise ValueError(
+                        f"Variable {node} present in multiple tiers. Aborting"
+                    )
+                temporal_ordering[node] = order
+
+        return temporal_ordering
+
+    def _orient_temporal_forbidden_edges(self, graph):
+        """
+        Add edge directions forbidden by the temporal order to forbidden_edges.
+
+        Parameters
+        ----------
+        graph: variable
+            The graph for which temporal order is specified.
+        """
+        if self.temporal_ordering == dict():
+            return
+
+        forbidden_edges = []
+        for node in graph.nodes:
+            for neighbor in graph.neighbors(node):
+                if self.temporal_ordering[neighbor] < self.temporal_ordering[node]:
+                    forbidden_edges.append((node, neighbor))
+
+        self.forbidden_edges = self.forbidden_edges.union(forbidden_edges)
+
     def __init__(
         self,
         forbidden_edges=None,
@@ -81,22 +180,23 @@ class ExpertKnowledge:
             else set()
         )
 
-        if temporal_order is not None:
-            raise ValueError(f"Specification of temporal order isn't supported yet.")
+        self.temporal_order = temporal_order if temporal_order is not None else [[]]
+        self.temporal_ordering = self._get_temporal_ordering(self.temporal_order)
 
     def apply_expert_knowledge(self, pdag):
         """
         Method to check consistency and orient edges in a graph based on expert knowledge.
 
-        The required and forbidden edges, if specified by the user, are correctly
-        oriented in the graph object passed. In case of any conflict between the
-        graph structure and a required/forbidden edge, the edge is ignored and
-        a warning is raised.
+        The required and forbidden edges, if specified by the user,
+        are correctly oriented in the graph object passed. Temporal order,
+        as specified, is also taken into account. In case of any conflict
+        between the graph structure and a required/forbidden edge, the edge is
+        ignored and a warning is raised.
 
         Parameters
         ----------
         pdag: pgmpy.base.PDAG
-            A  partial DAG with directed and undirected edges.
+            A partial DAG with directed and undirected edges.
 
         Returns
         --------
@@ -108,6 +208,8 @@ class ExpertKnowledge:
         ----------
         [1] https://doi.org/10.48550/arXiv.2306.01638
         """
+        self._validate_temporal_order(pdag.nodes())
+        self._orient_temporal_forbidden_edges(pdag)
 
         for edge in self.forbidden_edges:
             u, v = edge
diff --git a/pgmpy/estimators/HillClimbSearch.py b/pgmpy/estimators/HillClimbSearch.py
index aab52391..e7acbb4e 100644
--- a/pgmpy/estimators/HillClimbSearch.py
+++ b/pgmpy/estimators/HillClimbSearch.py
@@ -215,6 +215,7 @@ class HillClimbSearch(StructureEstimator):
 
         score, score_c = get_scoring_method(scoring_method, self.data, self.use_cache)
         score_fn = score_c.local_score
+
         # Step 1.2: Check the start_dag
         if start_dag is None:
             start_dag = DAG()
@@ -231,16 +232,12 @@ class HillClimbSearch(StructureEstimator):
             expert_knowledge = ExpertKnowledge()
 
         # Step 1.4: Check if required edges cause a cycle
-        start_dag.add_edges_from(
-            expert_knowledge.required_edges
-        )  # check for adding redundant edges?
+        start_dag.add_edges_from(expert_knowledge.required_edges)
         if not nx.is_directed_acyclic_graph(start_dag):
             raise ValueError(
                 "required_edges create a cycle in start_dag. Please modify either required_edges or start_dag."
             )
-        start_dag.remove_edges_from(
-            expert_knowledge.forbidden_edges
-        )  # check cycle here?
+        start_dag.remove_edges_from(expert_knowledge.forbidden_edges)
 
         # Step 1.5: Initialize max_indegree, tabu_list, and progress bar
         if max_indegree is None:
diff --git a/pgmpy/estimators/PC.py b/pgmpy/estimators/PC.py
index 3cde742b..a177ccce 100644
--- a/pgmpy/estimators/PC.py
+++ b/pgmpy/estimators/PC.py
@@ -195,15 +195,24 @@ class PC(StructureEstimator):
         if return_type.lower() == "skeleton":
             return skel, separating_sets
 
-        # Step 2: Orient the edges based on Meek's rules to build the PDAG/CPDAG.
-        pdag = self.orient_colliders(skel, separating_sets)
-        pdag = self.apply_orientation_rules(pdag)
+        # Step 2: Orient the edges based on collider structures.
+        pdag = self.orient_colliders(
+            skel, separating_sets, expert_knowledge.temporal_ordering
+        )
 
         # Step 3: Either return the CPDAG, integrate expert knowledge or fully orient the edges to build a DAG.
-        if not enforce_expert_knowledge:
+        if expert_knowledge.temporal_order != [[]]:
+            pdag = expert_knowledge.apply_expert_knowledge(pdag)
+            pdag = self.apply_orientation_rules(pdag, apply_r4=True)
+
+        elif not enforce_expert_knowledge:
+            pdag = self.apply_orientation_rules(pdag)
             pdag = expert_knowledge.apply_expert_knowledge(pdag)
             pdag = self.apply_orientation_rules(pdag, apply_r4=True)
 
+        else:
+            pdag = self.apply_orientation_rules(pdag)
+
         if self.data is not None:
             pdag.add_nodes_from(set(self.data.columns) - set(pdag.nodes()))
 
@@ -275,6 +284,7 @@ class PC(StructureEstimator):
 
         # Step 1: Initialize a fully connected undirected graph
         graph = nx.complete_graph(n=self.variables, create_using=nx.Graph)
+        temporal_ordering = expert_knowledge.temporal_ordering
         if enforce_expert_knowledge:
             graph.remove_edges_from(expert_knowledge.forbidden_edges)
 
@@ -290,13 +300,8 @@ class PC(StructureEstimator):
                     if (enforce_expert_knowledge is False) or (
                         (u, v) not in expert_knowledge.required_edges
                     ):
-                        for separating_set in chain(
-                            combinations(
-                                set(graph.neighbors(u)) - set([v]), lim_neighbors
-                            ),
-                            combinations(
-                                set(graph.neighbors(v)) - set([u]), lim_neighbors
-                            ),
+                        for separating_set in PC._get_potential_sepsets(
+                            u, v, temporal_ordering, graph, lim_neighbors
                         ):
                             # If a conditioning set exists remove the edge, store the separating set
                             # and move on to finding conditioning set for next edge.
@@ -320,9 +325,8 @@ class PC(StructureEstimator):
                     if (enforce_expert_knowledge is False) or (
                         (u, v) not in expert_knowledge.required_edges
                     ):
-                        for separating_set in chain(
-                            combinations(set(neighbors[u]) - set([v]), lim_neighbors),
-                            combinations(set(neighbors[v]) - set([u]), lim_neighbors),
+                        for separating_set in PC._get_potential_sepsets(
+                            u, v, temporal_ordering, graph, lim_neighbors
                         ):
                             # If a conditioning set exists remove the edge, store the
                             # separating set and move on to finding conditioning set for next edge.
@@ -343,9 +347,8 @@ class PC(StructureEstimator):
                 neighbors = {node: set(graph[node]) for node in graph.nodes()}
 
                 def _parallel_fun(u, v):
-                    for separating_set in chain(
-                        combinations(set(graph.neighbors(u)) - set([v]), lim_neighbors),
-                        combinations(set(graph.neighbors(v)) - set([u]), lim_neighbors),
+                    for separating_set in PC._get_potential_sepsets(
+                        u, v, temporal_ordering, graph, lim_neighbors
                     ):
                         if ci_test(
                             u,
@@ -394,6 +397,57 @@ class PC(StructureEstimator):
             pbar.close()
         return graph, separating_sets
 
+    @staticmethod
+    def _get_potential_sepsets(u, v, temporal_ordering, graph, lim_neighbors):
+        """
+        Return the temporally consistent superset of separating set of u, v.
+
+        The temporal order (if specified) of the superset can only be smaller
+        ("earlier") than the particular node. The neighbors of 'u' satisfying
+        this condition are returned.
+
+        Parameters
+        ----------
+        u: variable
+            The node whose neighbors are being considered for separating set.
+
+        v: variable
+            The node along with u whose separating set is being calculated.
+
+        temporal_ordering: dict
+            The temporal ordering of variables according to prior knowledgee.
+
+        graph: UndirectedGraph
+            The graph where separating sets are being calculated for the edges.
+
+        lim_neighbors: int
+            The maximum number of neighbours (conditioning variables) for u, v.
+
+        Returns
+        --------
+        separating_set: set
+            Set containing the superset of separating set of u, v.
+        """
+        separating_set_u = set(graph.neighbors(u))
+        separating_set_v = set(graph.neighbors(v))
+        separating_set_u.discard(v)
+        separating_set_v.discard(u)
+
+        if temporal_ordering != dict():
+            max_order = min(temporal_ordering[u], temporal_ordering[u])
+            for neigh in list(separating_set_u):
+                if temporal_ordering[neigh] > max_order:
+                    separating_set_u.discard(neigh)
+
+            for neigh in list(separating_set_v):
+                if temporal_ordering[neigh] > max_order:
+                    separating_set_v.discard(neigh)
+
+        return chain(
+            combinations(separating_set_u, lim_neighbors),
+            combinations(separating_set_v, lim_neighbors),
+        )
+
     @staticmethod
     def _check_incoming_edges(pdag, u, v):
         "Used for checking whether a new v-structure is getting formed"
@@ -413,7 +467,7 @@ class PC(StructureEstimator):
         return False
 
     @staticmethod
-    def orient_colliders(skeleton, separating_sets):
+    def orient_colliders(skeleton, separating_sets, temporal_ordering=dict()):
         """
         Orients the edges that form v-structures in a graph skeleton
         based on information from `separating_sets` to form a DAG pattern (PDAG).
@@ -467,7 +521,11 @@ class PC(StructureEstimator):
             if not skeleton.has_edge(X, Y):
                 for Z in set(skeleton.neighbors(X)) & set(skeleton.neighbors(Y)):
                     if Z not in separating_sets[frozenset((X, Y))]:
-                        pdag.remove_edges_from([(Z, X), (Z, Y)])
+                        if (temporal_ordering == dict()) or (
+                            (temporal_ordering[Z] >= temporal_ordering[X])
+                            and (temporal_ordering[Z] >= temporal_ordering[Y])
+                        ):
+                            pdag.remove_edges_from([(Z, X), (Z, Y)])
 
         edges = set(pdag.edges())
         undirected_edges = []
diff --git a/pgmpy/factors/continuous/LinearGaussianCPD.py b/pgmpy/factors/continuous/LinearGaussianCPD.py
index c9ff785c..4dc2a0ff 100644
--- a/pgmpy/factors/continuous/LinearGaussianCPD.py
+++ b/pgmpy/factors/continuous/LinearGaussianCPD.py
@@ -125,3 +125,54 @@ class LinearGaussianCPD(BaseFactor):
     def __repr__(self):
         str_repr = self.__str__()
         return f"<LinearGaussianCPD: {str_repr} at {hex(id(self))}"
+
+    @staticmethod
+    def get_random(variable, evidence, loc=0.0, scale=1.0, seed=None):
+        """
+        Generates a LinearGaussianCPD instance with random values on `variable` with
+        parents/evidence `evidence` with beta and std sampled from loc and scale
+
+        Parameters
+        ----------
+        variable: str, int or any hashable python object.
+            The variable on which to define the TabularCPD.
+
+        evidence: list, array-like
+            A list of variable names which are the parents/evidence of `variable`.
+
+        loc: float
+            The mean of the normal distribution from which the coefficients are
+            sampled.
+
+        scale: float
+            The standard deviation of the normal distribution from which the
+            coefficients are sampled.
+
+        seed: int (default: None)
+            The seed for the random number generator.
+
+        Returns
+        -------
+        Random CPD: pgmpy.factors.continuous.LinearGaussianCPD
+            A LinearGaussianCPD object on `variable` with `evidence` as evidence with random values.
+
+        Examples
+        --------
+        >>> from pgmpy.factors.continuous import LinearGaussianCPD
+        >>> LinearGaussianCPD.get_random(variable='Income', evidence=['Age', 'Experience'],
+        ...            loc=2.0, scale=0.5, seed=5)
+        <LinearGaussianCPD: P(Income | Age, Experience) = N(1.338*Age + 1.876*Experience + 1.599; 2.21) at 0x1795561e0
+        """
+        rng = np.random.default_rng(seed=seed)
+
+        beta = rng.normal(loc=loc, scale=scale, size=(len(evidence) + 1))
+        std = abs(rng.normal(loc=loc, scale=scale))
+
+        node_cpd = LinearGaussianCPD(
+            variable=variable,
+            beta=beta,
+            std=std,
+            evidence=evidence,
+        )
+
+        return node_cpd
diff --git a/pgmpy/models/BayesianNetwork.py b/pgmpy/models/BayesianNetwork.py
index d1ae7b42..fee4e2b5 100644
--- a/pgmpy/models/BayesianNetwork.py
+++ b/pgmpy/models/BayesianNetwork.py
@@ -91,12 +91,17 @@ class BayesianNetwork(DAG):
     3
     """
 
-    def __init__(self, ebunch=None, latents=set()):
-        super(BayesianNetwork, self).__init__(ebunch=ebunch, latents=latents)
+    def __init__(self, ebunch=None, latents=set(), lavaan_str=None, dagitty_str=None):
+        super(BayesianNetwork, self).__init__(
+            ebunch=ebunch,
+            latents=latents,
+            lavaan_str=lavaan_str,
+            dagitty_str=dagitty_str,
+        )
         self.cpds = []
         self.cardinalities = defaultdict(int)
 
-    def add_edge(self, u, v, **kwargs):
+    def add_edge(self, u, v, w=None, **kwargs):
         """
         Add an edge between u and v.
 
@@ -123,7 +128,10 @@ class BayesianNetwork(DAG):
                 % (u, v)
             )
         else:
-            super(BayesianNetwork, self).add_edge(u, v, **kwargs)
+            if w:
+                super(BayesianNetwork, self).add_edge(u, v, w, **kwargs)
+            else:
+                super(BayesianNetwork, self).add_edge(u, v, **kwargs)
 
     def remove_node(self, node):
         """
diff --git a/pgmpy/models/LinearGaussianBayesianNetwork.py b/pgmpy/models/LinearGaussianBayesianNetwork.py
index 483297cd..2b4a5bd1 100644
--- a/pgmpy/models/LinearGaussianBayesianNetwork.py
+++ b/pgmpy/models/LinearGaussianBayesianNetwork.py
@@ -21,9 +21,12 @@ class LinearGaussianBayesianNetwork(DAG):
 
     """
 
-    def __init__(self, ebunch=None, latents=set()):
+    def __init__(self, ebunch=None, latents=set(), lavaan_str=None, dagitty_str=None):
         super(LinearGaussianBayesianNetwork, self).__init__(
-            ebunch=ebunch, latents=latents
+            ebunch=ebunch,
+            latents=latents,
+            lavaan_str=lavaan_str,
+            dagitty_str=dagitty_str,
         )
         self.cpds = []
 
@@ -162,17 +165,19 @@ class LinearGaussianBayesianNetwork(DAG):
         seed: int
             The seed for the random number generator.
         """
-        rng = np.random.default_rng(seed=seed)
+        # We want to provide a different seed for each cpd, therefore we force it to be integer and increment in a loop.
+        seed = seed if seed else 42
 
         cpds = []
-        for var in self.nodes():
+        for i, var in enumerate(self.nodes()):
             parents = self.get_parents(var)
             cpds.append(
-                LinearGaussianCPD(
+                LinearGaussianCPD.get_random(
                     variable=var,
-                    beta=rng.normal(loc=loc, scale=scale, size=(len(parents) + 1)),
-                    std=abs(rng.normal(loc=loc, scale=scale)),
                     evidence=parents,
+                    loc=loc,
+                    scale=scale,
+                    seed=(seed + i),
                 )
             )
         if inplace:


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_PC.py b/pgmpy/tests/test_estimators/test_PC.py
index ba7dd411..d117de72 100644
--- a/pgmpy/tests/test_estimators/test_PC.py
+++ b/pgmpy/tests/test_estimators/test_PC.py
@@ -636,3 +636,83 @@ class TestPCRealModels(unittest.TestCase):
                 "WARNING:pgmpy:Specified expert knowledge conflicts with learned structure. Ignoring edge xray->either from required edges"
             ],
         )
+
+    def test_temporal_pc_cancer(self):
+        cancer_model = get_example_model("cancer")
+        data = BayesianModelSampling(cancer_model).forward_sample(
+            size=int(5e4), seed=42
+        )
+        est = PC(data)
+        background = ExpertKnowledge(  # e.g. we only know "Pollution", "Smoker", "Cancer" can be the causes of others
+            temporal_order=[["Pollution", "Smoker", "Cancer"], ["Dyspnoea", "Xray"]],
+            max_cond_vars=4,
+        )
+        pdag = est.estimate(
+            variant="stable",
+            expert_knowledge=background,
+            n_jobs=2,
+            show_progress=False,
+        )
+        self.assertSetEqual(
+            set(pdag.edges()),
+            set(
+                [
+                    ("Cancer", "Xray"),
+                    ("Cancer", "Dyspnoea"),
+                    ("Smoker", "Cancer"),
+                    ("Pollution", "Cancer"),
+                ]
+            ),
+        )
+
+    def test_temporal_pc_sachs(self):
+        temporal_order = [
+            ["PKC", "Plcg"],
+            [
+                "PKA",
+                "Raf",
+                "Jnk",
+                "P38",
+                "PIP3",
+                "PIP2",
+                "Mek",
+                "Erk",
+            ],
+            ["Akt"],
+        ]
+        temporal_forbidden_edges = set(
+            [
+                ("PKA", "PKC"),
+                ("PKA", "Plcg"),
+                ("Raf", "PKC"),
+                ("Raf", "Plcg"),
+                ("Jnk", "PKC"),
+                ("Jnk", "Plcg"),
+                ("P38", "PKC"),
+                ("P38", "Plcg"),
+                ("PIP3", "PKC"),
+                ("PIP3", "Plcg"),
+                ("PIP2", "PKC"),
+                ("PIP2", "Plcg"),
+                ("Mek", "PKC"),
+                ("Mek", "Plcg"),
+                ("Erk", "PKC"),
+                ("Erk", "Plcg"),
+                ("Akt", "PKC"),
+                ("Akt", "Plcg"),
+                ("Akt", "PKA"),
+                ("Akt", "Raf"),
+                ("Akt", "Jnk"),
+                ("Akt", "P38"),
+                ("Akt", "PIP3"),
+                ("Akt", "PIP2"),
+                ("Akt", "Mek"),
+                ("Akt", "Erk"),
+            ]
+        )
+
+        model = get_example_model("sachs")
+        df = model.simulate(int(1e3))
+        expert = ExpertKnowledge(temporal_order=temporal_order)
+        pdag = PC(df).estimate(ci_test="chi_square", expert_knowledge=expert)
+        self.assertTrue(temporal_forbidden_edges.isdisjoint(set(pdag.edges())))
diff --git a/pgmpy/tests/test_factors/test_continuous/test_Linear_Gaussian_CPD.py b/pgmpy/tests/test_factors/test_continuous/test_Linear_Gaussian_CPD.py
index 53339a49..c97391cc 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_Linear_Gaussian_CPD.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_Linear_Gaussian_CPD.py
@@ -24,3 +24,7 @@ class TestLGCPD(unittest.TestCase):
             cpd2.__str__(),
             "P(y | x1, x2, x3) = N(1.0*x1 + 4.56*x2 + 8.0*x3 + 0.67; 2)",
         )
+
+    def test_get_random(self):
+        cpd_random = LinearGaussianCPD.get_random("x", ["x1", "x2", "x3"], 0.23, 0.56)
+        self.assertIn("P(x | x1, x2, x3) = N(", cpd_random.__str__())
diff --git a/pgmpy/tests/test_models/test_BayesianNetwork.py b/pgmpy/tests/test_models/test_BayesianNetwork.py
index 4a5b96d8..8fc9dd69 100644
--- a/pgmpy/tests/test_models/test_BayesianNetwork.py
+++ b/pgmpy/tests/test_models/test_BayesianNetwork.py
@@ -42,6 +42,26 @@ class TestBaseModelCreation(unittest.TestCase):
     def test_class_init_with_data_nonstring(self):
         BayesianNetwork([(1, 2), (2, 3)])
 
+    def test_class_init_with_adj_matrix_dict_of_dict(self):
+        adj = {"a": {"b": 4, "c": 3}, "b": {"c": 2}}
+        self.graph = BayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+        self.assertEqual(self.graph.adj["a"]["c"]["weight"], 3)
+
+    def test_class_init_with_adj_matrix_dict_of_list(self):
+        adj = {"a": ["b", "c"], "b": ["c"]}
+        self.graph = BayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+
+    def test_class_init_with_pd_adj_df(self):
+        df = pd.DataFrame([[0, 3], [0, 0]])
+        self.graph = BayesianNetwork(df, latents=set([0]))
+        self.assertEqual(self.graph.latents, set([0]))
+        self.assertListEqual(sorted(self.graph.nodes()), [0, 1])
+        self.assertEqual(self.graph.adj[0][1]["weight"], {"weight": 3})
+
     def test_add_node_string(self):
         self.G.add_node("a")
         self.assertListEqual(list(self.G.nodes()), ["a"])
@@ -118,6 +138,23 @@ class TestBaseModelCreation(unittest.TestCase):
         del self.G
 
 
+class TestBayesianNetworkParser(unittest.TestCase):
+    def test_from_lavaan(self):
+        model_str = "i =~ x1 + x2 + x3"
+        model_from_str = BayesianNetwork.from_lavaan(string=model_str)
+        expected_edges = set([("i", "x1"), ("i", "x2"), ("i", "x3")])
+        expected_latents = set(["i"])
+        self.assertEqual(set(model_from_str.edges()), expected_edges)
+        self.assertEqual(set(model_from_str.latents), expected_latents)
+
+    def test_from_daggitty(self):
+        dag = BayesianNetwork.from_dagitty(
+            'dag{ bb="0,0,1,1" X [l, pos="-1.228,-1.145"] X-> {Y Z}  Z->A}'
+        )
+        self.assertEqual(set(dag.edges()), set([("X", "Z"), ("X", "Y"), ("Z", "A")]))
+        self.assertEqual(set(dag.latents), set(["X"]))
+
+
 class TestBayesianNetworkMethods(unittest.TestCase):
     def setUp(self):
         self.G = BayesianNetwork([("a", "d"), ("b", "d"), ("d", "e"), ("b", "c")])
diff --git a/pgmpy/tests/test_models/test_FunctionalBayesianNetwork.py b/pgmpy/tests/test_models/test_FunctionalBayesianNetwork.py
index 020cc57a..a60ba5f6 100644
--- a/pgmpy/tests/test_models/test_FunctionalBayesianNetwork.py
+++ b/pgmpy/tests/test_models/test_FunctionalBayesianNetwork.py
@@ -161,3 +161,39 @@ class TestFBNMethods(unittest.TestCase):
 
         self.assertTrue(np.all(samples["lognormal"] > 0))
         self.assertTrue(np.all(samples["gamma"] > 0))
+
+
+class TestFBNCreation(unittest.TestCase):
+    def test_class_init_with_adj_matrix_dict_of_dict(self):
+        adj = {"a": {"b": 4, "c": 3}, "b": {"c": 2}}
+        self.graph = FunctionalBayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+        self.assertEqual(self.graph.adj["a"]["c"]["weight"], 3)
+
+    def test_class_init_with_adj_matrix_dict_of_list(self):
+        adj = {"a": ["b", "c"], "b": ["c"]}
+        self.graph = FunctionalBayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+
+    def test_class_init_with_pd_adj_df(self):
+        df = pd.DataFrame([[0, 3], [0, 0]])
+        self.graph = FunctionalBayesianNetwork(df, latents=set([0]))
+        self.assertEqual(self.graph.latents, set([0]))
+        self.assertListEqual(sorted(self.graph.nodes()), [0, 1])
+        self.assertEqual(self.graph.adj[0][1]["weight"], {"weight": 3})
+
+
+class TestDAGParser(unittest.TestCase):
+    def test_from_lavaan(self):
+        model_str = "ind60 =~ x1"
+        model_from_str = FunctionalBayesianNetwork.from_lavaan(string=model_str)
+        expected_edges = set([("ind60", "x1")])
+        self.assertEqual(set(model_from_str.edges()), expected_edges)
+
+    def test_from_dagitty(self):
+        model_str = """dag{X -> Y}"""
+        model_from_str = FunctionalBayesianNetwork.from_dagitty(string=model_str)
+        expected_edges = set([("X", "Y")])
+        self.assertEqual(set(model_from_str.edges()), expected_edges)
diff --git a/pgmpy/tests/test_models/test_LinearGaussianBayesianNetwork.py b/pgmpy/tests/test_models/test_LinearGaussianBayesianNetwork.py
index ea2aacf1..a1a62ad5 100644
--- a/pgmpy/tests/test_models/test_LinearGaussianBayesianNetwork.py
+++ b/pgmpy/tests/test_models/test_LinearGaussianBayesianNetwork.py
@@ -223,3 +223,39 @@ class TestLGBNMethods(unittest.TestCase):
 
     def tearDown(self):
         del self.model, self.cpd1, self.cpd2, self.cpd3
+
+
+class TestLGBNCreation(unittest.TestCase):
+    def test_class_init_with_adj_matrix_dict_of_dict(self):
+        adj = {"a": {"b": 4, "c": 3}, "b": {"c": 2}}
+        self.graph = LinearGaussianBayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+        self.assertEqual(self.graph.adj["a"]["c"]["weight"], 3)
+
+    def test_class_init_with_adj_matrix_dict_of_list(self):
+        adj = {"a": ["b", "c"], "b": ["c"]}
+        self.graph = LinearGaussianBayesianNetwork(adj, latents=set(["a"]))
+        self.assertEqual(self.graph.latents, set("a"))
+        self.assertListEqual(sorted(self.graph.nodes()), ["a", "b", "c"])
+
+    def test_class_init_with_pd_adj_df(self):
+        df = pd.DataFrame([[0, 3], [0, 0]])
+        self.graph = LinearGaussianBayesianNetwork(df, latents=set([0]))
+        self.assertEqual(self.graph.latents, set([0]))
+        self.assertListEqual(sorted(self.graph.nodes()), [0, 1])
+        self.assertEqual(self.graph.adj[0][1]["weight"], {"weight": 3})
+
+
+class TestDAGParser(unittest.TestCase):
+    def test_from_lavaan(self):
+        model_str = "d ~ i"
+        model_from_str = LinearGaussianBayesianNetwork.from_lavaan(string=model_str)
+        expected_edges = set([("i", "d")])
+        self.assertEqual(set(model_from_str.edges()), expected_edges)
+
+    def test_from_dagitty(self):
+        model_str = """dag{ smoking -> "carry matches" }"""
+        model_from_str = LinearGaussianBayesianNetwork.from_dagitty(string=model_str)
+        expected_edges = set([("smoking", "carry matches")])
+        self.assertEqual(set(model_from_str.edges()), expected_edges)