Repository: pgmpy/pgmpy

Issue:
Move BayesianNetwork.fit method to DAG.fit
As structure learning algorithms  output either a `pgmpy.base.PDAG` or `pgmpy.base.DAG` object, to fit parameters to these models users currently need to convert them to a `pgmpy.models.BayesianNetwork` object. By moving the `BayesianNetwork.fit` method to `DAG.fit`, it would make the workflow simpler and more straightforward. As BayesianNetwork inherits DAG, there would be no change in the Bayesian Network API.

Discussion:

Solution patch:
diff --git a/CITATION.cff b/CITATION.cff
index 9ad6484a6..c8bac8974 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -3,22 +3,21 @@ message: "If you use this software, please cite it as below."
 authors:
 - family-names: "Ankan"
   given-names: "Ankur"
-- family-names: "Panda"
-  given-names: "Abinash"
-title: "pgmpy: Probabilistic Graphical Models using Python"
-version: 0.1.25
+- family-names: "Textor"
+  given-names: "Johannes"
+title: "pgmpy: A Python Toolkit for Bayesian Networks"
+version: 0.1.26
 url: "https://github.com/pgmpy/pgmpy"
 preferred-citation:
-  type: conference-paper
+  type: article
   authors:
   - family-names: "Ankan"
     given-names: "Ankur"
-  - family-names: "Panda"
-    given-names: "Abinash"
-  doi: "10.25080/majora-7b98e3ed-001"
-  journal: "Proceedings of the 14th Python in Science Conference"
-  start: 6
-  end: 11
-  title: "pgmpy: Probabilistic Graphical Models using Python"
-  issue: 14
-  year: 2015
+  - family-names: "Textor"
+    given-names: "Johannes"
+  journal: "Journal Of Machine Learning Research"
+  start: 1
+  end: 8
+  title: "pgmpy: A Python Toolkit for Bayesian Networks"
+  issue: 25
+  year: 2024
diff --git a/README.md b/README.md
index 3a759bde8..c0ebe65dc 100644
--- a/README.md
+++ b/README.md
@@ -50,19 +50,20 @@ Citing
 If you use `pgmpy` in your scientific work, please consider citing us:
 
 ```
-Ankur Ankan, & Abinash Panda ( 2015 ). pgmpy: Probabilistic Graphical Models using Python. In Proceedings of the 14th Python in Science Conference (pp. 6 - 11 ).
+Ankur Ankan, & Johannes Textor (2024). pgmpy: A Python Toolkit for Bayesian Networks. Journal of Machine Learning Research, 25(265), 1–8.
 ```
 
 Bibtex:
 ```
-@InProceedings{ Ankan2015,
-  author    = { {A}nkur {A}nkan and {A}binash {P}anda },
-  title     = { pgmpy: {P}robabilistic {G}raphical {M}odels using {P}ython },
-  booktitle = { {P}roceedings of the 14th {P}ython in {S}cience {C}onference },
-  pages     = { 6 - 11 },
-  year      = { 2015 },
-  editor    = { {K}athryn {H}uff and {J}ames {B}ergstra },
-  doi       = { 10.25080/Majora-7b98e3ed-001 }
+@article{Ankan2024,
+  author  = {Ankur Ankan and Johannes Textor},
+  title   = {pgmpy: A Python Toolkit for Bayesian Networks},
+  journal = {Journal of Machine Learning Research},
+  year    = {2024},
+  volume  = {25},
+  number  = {265},
+  pages   = {1--8},
+  url     = {http://jmlr.org/papers/v25/23-0487.html}
 }
 ```
 
diff --git a/pgmpy/base/DAG.py b/pgmpy/base/DAG.py
index 841f55077..9e6a243ff 100644
--- a/pgmpy/base/DAG.py
+++ b/pgmpy/base/DAG.py
@@ -1090,6 +1090,77 @@ def to_graphviz(self):
         """
         return nx.nx_agraph.to_agraph(self)
 
+    def fit(self, data, estimator=None, state_names=[], n_jobs=1, **kwargs):
+        """
+        Estimates the CPD for each variable based on a given data set.
+
+        Parameters
+        ----------
+        data: pandas DataFrame object
+            DataFrame object with column names identical to the variable names of the network.
+            (If some values in the data are missing the data cells should be set to `numpy.nan`.
+            Note that pandas converts each column containing `numpy.nan`s to dtype `float`.)
+
+        estimator: Estimator class
+            One of:
+            - MaximumLikelihoodEstimator (default)
+            - BayesianEstimator: In this case, pass 'prior_type' and either 'pseudo_counts'
+            or 'equivalent_sample_size' as additional keyword arguments.
+            See `BayesianEstimator.get_parameters()` for usage.
+            - ExpectationMaximization
+
+        state_names: dict (optional)
+            A dict indicating, for each variable, the discrete set of states
+            that the variable can take. If unspecified, the observed values
+            in the data set are taken to be the only possible states.
+
+        n_jobs: int (default: 1)
+            Number of threads/processes to use for estimation. Using n_jobs > 1
+            for small models or datasets might be slower.
+
+        Returns
+        -------
+        Fitted Model: BayesianNetwork
+            Returns a BayesianNetwork object with learned CPDs.
+            The DAG structure is preserved, and parameters (CPDs) are added.
+            This allows the DAG to represent both the structure and the parameters of a Bayesian Network.
+
+        Examples
+        --------
+        >>> import pandas as pd
+        >>> from pgmpy.models import BayesianNetwork
+        >>> from pgmpy.base import DAG
+        >>> data = pd.DataFrame(data={'A': [0, 0, 1], 'B': [0, 1, 0], 'C': [1, 1, 0]})
+        >>> model = DAG([('A', 'C'), ('B', 'C')])
+        >>> fitted_model = model.fit(data)
+        >>> fitted_model.get_cpds()
+        [<TabularCPD representing P(A:2) at 0x17945372c30>,
+        <TabularCPD representing P(B:2) at 0x17945a19760>,
+        <TabularCPD representing P(C:2 | A:2, B:2) at 0x17944f42690>]
+        """
+        from pgmpy.estimators import BaseEstimator, MaximumLikelihoodEstimator
+        from pgmpy.models import BayesianNetwork
+
+        if isinstance(self, BayesianNetwork):
+            bn = self
+        else:
+            bn = BayesianNetwork(self.edges())
+
+        if estimator is None:
+            estimator = MaximumLikelihoodEstimator
+        else:
+            if not issubclass(estimator, BaseEstimator):
+                raise TypeError("Estimator object should be a valid pgmpy estimator.")
+
+        _estimator = estimator(
+            bn,
+            data,
+            state_names=state_names,
+        )
+        cpds_list = _estimator.get_parameters(n_jobs=n_jobs, **kwargs)
+        bn.add_cpds(*cpds_list)
+        return bn
+
 
 class PDAG(nx.DiGraph):
     """
diff --git a/pgmpy/estimators/BayesianEstimator.py b/pgmpy/estimators/BayesianEstimator.py
index 0af556371..97ad12642 100644
--- a/pgmpy/estimators/BayesianEstimator.py
+++ b/pgmpy/estimators/BayesianEstimator.py
@@ -6,6 +6,7 @@
 import numpy as np
 from joblib import Parallel, delayed
 
+from pgmpy.base import DAG
 from pgmpy.estimators import ParameterEstimator
 from pgmpy.factors.discrete import TabularCPD
 from pgmpy.global_vars import logger
@@ -19,15 +20,20 @@ class BayesianEstimator(ParameterEstimator):
     """
 
     def __init__(self, model, data, **kwargs):
-        if not isinstance(model, BayesianNetwork):
+        if not isinstance(model, (DAG, BayesianNetwork)):
             raise NotImplementedError(
-                "Bayesian Parameter Estimation is only implemented for BayesianNetwork"
-            )
-        elif len(model.latents) != 0:
-            raise ValueError(
-                f"Bayesian Parameter Estimation works only on models with all observed variables. Found latent variables: {model.latents}"
+                "Bayesian Parameter Estimation is only implemented for DAG or BayesianNetwork"
             )
 
+        if isinstance(model, (DAG, BayesianNetwork)):
+            if len(model.latents) != 0:
+                raise ValueError(
+                    f"Bayesian Parameter Estimation works only on models with all observed variables. Found latent variables: {model.latents}"
+                )
+
+            if isinstance(model, DAG):
+                model = BayesianNetwork(model.edges())
+
         super(BayesianEstimator, self).__init__(model, data, **kwargs)
 
     def get_parameters(
diff --git a/pgmpy/estimators/EM.py b/pgmpy/estimators/EM.py
index 25539da94..83850a5a7 100644
--- a/pgmpy/estimators/EM.py
+++ b/pgmpy/estimators/EM.py
@@ -10,6 +10,7 @@
 from pgmpy.estimators import MaximumLikelihoodEstimator, ParameterEstimator
 from pgmpy.factors.discrete import TabularCPD
 from pgmpy.models import BayesianNetwork
+from pgmpy.base import DAG
 
 
 class ExpectationMaximization(ParameterEstimator):
@@ -50,11 +51,17 @@ class ExpectationMaximization(ParameterEstimator):
     """
 
     def __init__(self, model, data, **kwargs):
-        if not isinstance(model, BayesianNetwork):
+        if not isinstance(model, (DAG, BayesianNetwork)):
             raise NotImplementedError(
-                "Expectation Maximization is only implemented for BayesianNetwork"
+                "Expectation Maximization is only implemented for DAG or BayesianNetwork"
             )
 
+        if isinstance(model, DAG):
+            model_bn = BayesianNetwork(model.edges())
+            model_bn.add_nodes_from(model.nodes())
+            model_bn.latents = model.latents
+            model = model_bn
+
         super(ExpectationMaximization, self).__init__(model, data, **kwargs)
         self.model_copy = self.model.copy()
 
diff --git a/pgmpy/estimators/ExhaustiveSearch.py b/pgmpy/estimators/ExhaustiveSearch.py
index 077d007df..086f40b90 100644
--- a/pgmpy/estimators/ExhaustiveSearch.py
+++ b/pgmpy/estimators/ExhaustiveSearch.py
@@ -19,8 +19,8 @@ class ExhaustiveSearch(StructureEstimator):
     ----------
     data: pandas DataFrame object
         dataframe object where each column represents one variable.
-        (If some values in the data are missing the data cells should be set to `numpy.nan`.
-        Note that pandas converts each column containing `numpy.nan`s to dtype `float`.)
+        (If some values in the data are missing the data cells should be set to `numpy.NaN`.
+        Note that pandas converts each column containing `numpy.NaN`s to dtype `float`.)
 
     scoring_method: Instance of a `StructureScore`-subclass (`K2Score` is used as default)
         An instance of `K2Score`, `BDeuScore`, `BicScore` or 'AICScore'.
diff --git a/pgmpy/estimators/MLE.py b/pgmpy/estimators/MLE.py
index bceb4b1e1..3076c1f7f 100644
--- a/pgmpy/estimators/MLE.py
+++ b/pgmpy/estimators/MLE.py
@@ -9,6 +9,7 @@
 from pgmpy.factors import FactorDict
 from pgmpy.factors.discrete import TabularCPD
 from pgmpy.models import BayesianNetwork, JunctionTree
+from pgmpy.base import DAG
 
 
 class MaximumLikelihoodEstimator(ParameterEstimator):
@@ -42,24 +43,30 @@ class MaximumLikelihoodEstimator(ParameterEstimator):
     """
 
     def __init__(self, model, data, **kwargs):
-        if not isinstance(model, BayesianNetwork) and not isinstance(
-            model, JunctionTree
-        ):
+        if not isinstance(model, (BayesianNetwork, JunctionTree, DAG)):
             raise NotImplementedError(
-                "Maximum Likelihood Estimate is only implemented for BayesianNetwork and JunctionTree"
+                "Maximum Likelihood Estimate is only implemented for BayesianNetwork, JunctionTree, and DAG"
             )
 
-        elif set(model.nodes()) > set(data.columns):
-            if isinstance(model, BayesianNetwork):
+        if isinstance(model, (DAG, BayesianNetwork)):
+            if len(model.latents) > 0:
                 raise ValueError(
-                    f"Found latent variables: {model.latents}. Maximum Likelihood doesn't support latent variables, please use ExpectationMaximization"
+                    f"Found latent variables: {model.latents}. Maximum Likelihood doesn't support latent variables, please use ExpectationMaximization."
                 )
-            else:
+
+            elif set(model.nodes()) > set(data.columns):
                 raise ValueError(
-                    "Nodes detected in the model that are not present in the dataset. "
-                    + "Refine the model so that all parameters can be estimated from the data."
+                    f"Nodes detected in the model that are not present in the dataset: {set(model.nodes) - set(data.columns)}. "
+                    "Refine the model so that all parameters can be estimated from the data."
                 )
 
+            if isinstance(model, JunctionTree):
+                if len(set(data.columns) - set(chain(*model.nodes()))) != 0:
+                    raise ValueError(
+                        f"Nodes detected in the JunctionTree that are not present in the dataset: {set(data.columns) - set(chain(*model.nodes()))} "
+                        "Refine the model to ensure all parameters can be estimated."
+                    )
+
         super(MaximumLikelihoodEstimator, self).__init__(model, data, **kwargs)
 
     def get_parameters(self, n_jobs=1, weighted=False):
diff --git a/pgmpy/factors/continuous/LinearGaussianCPD.py b/pgmpy/factors/continuous/LinearGaussianCPD.py
index c7abaf6f2..dc9d7c9c8 100644
--- a/pgmpy/factors/continuous/LinearGaussianCPD.py
+++ b/pgmpy/factors/continuous/LinearGaussianCPD.py
@@ -46,7 +46,8 @@ def __init__(
         variable: any hashable python object
             The variable whose CPD is defined.
 
-        evidence_mean: Mean vector (numpy array) of the joint distribution, X
+        evidence_mean: list (array-like)
+            Mean vector (numpy array) of the joint distribution, X
 
         evidence_variance: int, float
             The variance of the multivariate gaussian, X = ['x1', 'x2', ..., 'xn']
@@ -58,7 +59,6 @@ def __init__(
             An iterable representing the coefficient vector of the linear equation.
             The first term represents the constant term in the linear equation.
 
-
         Examples
         --------
 
diff --git a/pgmpy/inference/ExactInference.py b/pgmpy/inference/ExactInference.py
index 21c979a3c..3067b9a71 100644
--- a/pgmpy/inference/ExactInference.py
+++ b/pgmpy/inference/ExactInference.py
@@ -732,8 +732,10 @@ def _update_beliefs(self, sending_clique, receiving_clique, operation):
         ----------
         sending_clique: node (as the operation is on junction tree, node should be a tuple)
             Node sending the message
+
         receiving_clique: node (as the operation is on junction tree, node should be a tuple)
             Node receiving the message
+
         operation: str ('marginalize' | 'maximize')
             The operation to do for passing messages between nodes.
 
diff --git a/pgmpy/metrics/__init__.py b/pgmpy/metrics/__init__.py
index a33d19cce..e50e519bb 100644
--- a/pgmpy/metrics/__init__.py
+++ b/pgmpy/metrics/__init__.py
@@ -1,5 +1,6 @@
 from .bn_inference import BayesianModelProbability
 from .metrics import (
+    SHD,
     correlation_score,
     fisher_c,
     implied_cis,
@@ -13,4 +14,5 @@
     "structure_score",
     "implied_cis",
     "fisher_c",
+    "SHD",
 ]
diff --git a/pgmpy/metrics/metrics.py b/pgmpy/metrics/metrics.py
index 0e9702465..8a190bffc 100644
--- a/pgmpy/metrics/metrics.py
+++ b/pgmpy/metrics/metrics.py
@@ -1,6 +1,7 @@
 import math
 from itertools import combinations
 
+import networkx as nx
 import numpy as np
 import pandas as pd
 from scipy import stats
@@ -394,3 +395,71 @@ def fisher_c(model, data, ci_test, show_progress=True):
     C = -2 * np.log(cis.loc[:, "p_value"]).sum()
     p_value = 1 - stats.chi2.cdf(C, df=2 * cis.shape[0])
     return p_value
+
+
+def SHD(true_model, est_model):
+    """
+    Computes the Structural Hamming Distance between `true_model` and `est_model`.
+
+    SHD is defined as total number of basic operations: adding edges, removing
+    edges, and reversing edges required to transform one graph to the other. It
+    is a symmetrical measure.
+
+    The code first accounts for edges that need to be deleted (from true_model),
+    added (to true_model) and finally edges that need to be reversed. All operations
+    count as 1.
+
+    Parameters
+    ----------
+    true_model: pgmpy.base.DAG or pgmpy.base.CPDAG or pgmpy.models.BayesianNetwork
+        The first model to compare.
+    est_model: pgmpy.base.DAG or pgmpy.base.CPDAG or pgmpy.models.BayesianNetwork
+        The second model to compare.
+
+    Returns
+    -------
+    int:
+        If both true_model and est_model are DAGs or Bayesian Networks returns
+        an integer.
+
+    Examples
+    --------
+    >>> from pgmpy.metrics import SHD
+    >>> from pgmpy.models import BayesianNetwork
+    >>> dag1 = BayesianNetwork([(1, 2), (2, 3)])
+    >>> dag2 = BayesianNetwork([(2, 1), (2, 3)])
+    >>> SHD(dag1, dag2)
+    1
+    """
+    if set(true_model.nodes()) != set(est_model.nodes()):
+        raise ValueError("The graphs must have the same nodes.")
+
+    nodes_list = true_model.nodes()
+
+    dag_true = nx.DiGraph(true_model.edges())
+    m1 = nx.adjacency_matrix(dag_true, nodelist=nodes_list).todense()
+
+    dag_est = nx.DiGraph(est_model.edges())
+    m2 = nx.adjacency_matrix(dag_est, nodelist=nodes_list).todense()
+
+    shd = 0
+
+    s1 = m1 + m1.T
+    s2 = m2 + m2.T
+
+    # Edges that are in m1 but not in m2 (deletions from m1)
+    ds = s1 - s2
+    ind = np.where(ds > 0)
+    m1[ind] = 0
+    shd = shd + (len(ind[0]) / 2)
+
+    # Edges that are in m2 but not in m1 (additions to m1)
+    ind = np.where(ds < 0)
+    m1[ind] = m2[ind]
+    shd = shd + (len(ind[0]) / 2)
+
+    # Edges that need to be simply reversed
+    d = np.abs(m1 - m2)
+    shd = shd + (np.sum((d + d.T) > 0) / 2)
+
+    return int(shd)
diff --git a/pgmpy/models/BayesianNetwork.py b/pgmpy/models/BayesianNetwork.py
index f167b4400..486a828f4 100644
--- a/pgmpy/models/BayesianNetwork.py
+++ b/pgmpy/models/BayesianNetwork.py
@@ -279,18 +279,32 @@ def get_cpds(self, node=None):
 
         Returns
         -------
-        A list of TabularCPDs: list
+        cpd : TabularCPD object or list of TabularCPD objects
+            If 'node' is specified, returns the 'TabularCPD' object corresponding to the node.
+            If 'node' is not specified, returns a list of all 'TabularCPD' objects added to the model.
+
+        Raises
+        ------
+        ValueError
+            If the specified node is not present in the model.
 
         Examples
         --------
-        >>> from pgmpy.models import BayesianNetwork
-        >>> from pgmpy.factors.discrete import TabularCPD
-        >>> student = BayesianNetwork([('diff', 'grade'), ('intel', 'grade')])
-        >>> cpd = TabularCPD('grade', 2, [[0.1, 0.9, 0.2, 0.7],
-        ...                               [0.9, 0.1, 0.8, 0.3]],
-        ...                  ['intel', 'diff'], [2, 2])
-        >>> student.add_cpds(cpd)
-        >>> student.get_cpds()
+        >>> from pgmpy.utils import get_example_model
+        >>> model = get_example_model('asia')
+        >>> cpds = model.get_cpds()
+        >>> cpds
+        [<TabularCPD representing P(asia:2) at 0x7dbbd9bdbb80>,
+        <TabularCPD representing P(bronc:2 | smoke:2) at 0x7dbbd9bda3e0>,
+        <TabularCPD representing P(dysp:2 | bronc:2, either:2) at 0x7dbbd9bd8550>,
+        <TabularCPD representing P(either:2 | lung:2, tub:2) at 0x7dbbd9bda800>,
+        <TabularCPD representing P(lung:2 | smoke:2) at 0x7dbbd9bd89d0>,
+        <TabularCPD representing P(smoke:2) at 0x7dbbd9bd8f70>,
+        <TabularCPD representing P(tub:2 | asia:2) at 0x7dbbd9bda860>,
+        <TabularCPD representing P(xray:2 | either:2) at 0x7dbbd9bd9a80>]
+        >>> cpd = model.get_cpds('bronc')
+        >>> cpd
+        <TabularCPD representing P(bronc:2 | smoke:2) at 0x7dbbd9bda3e0>
         """
         if node is not None:
             if node not in self.nodes():
@@ -510,75 +524,6 @@ def to_junction_tree(self):
         mm = self.to_markov_model()
         return mm.to_junction_tree()
 
-    def fit(
-        self,
-        data,
-        estimator=None,
-        state_names=[],
-        n_jobs=1,
-        **kwargs,
-    ):
-        """
-        Estimates the CPD for each variable based on a given data set.
-
-        Parameters
-        ----------
-        data: pandas DataFrame object
-            DataFrame object with column names identical to the variable names of the network.
-            (If some values in the data are missing the data cells should be set to `numpy.nan`.
-            Note that pandas converts each column containing `numpy.nan`s to dtype `float`.)
-
-        estimator: Estimator class
-            One of:
-            - MaximumLikelihoodEstimator (default)
-            - BayesianEstimator: In this case, pass 'prior_type' and either 'pseudo_counts'
-            or 'equivalent_sample_size' as additional keyword arguments.
-            See `BayesianEstimator.get_parameters()` for usage.
-            - ExpectationMaximization
-
-        state_names: dict (optional)
-            A dict indicating, for each variable, the discrete set of states
-            that the variable can take. If unspecified, the observed values
-            in the data set are taken to be the only possible states.
-
-        n_jobs: int (default: 1)
-            Number of threads/processes to use for estimation. Using n_jobs > 1
-            for small models or datasets might be slower.
-
-        Returns
-        -------
-        Fitted Model: None
-            Modifies the network inplace and adds the `cpds` property.
-
-        Examples
-        --------
-        >>> import pandas as pd
-        >>> from pgmpy.models import BayesianNetwork
-        >>> from pgmpy.estimators import MaximumLikelihoodEstimator
-        >>> data = pd.DataFrame(data={'A': [0, 0, 1], 'B': [0, 1, 0], 'C': [1, 1, 0]})
-        >>> model = BayesianNetwork([('A', 'C'), ('B', 'C')])
-        >>> model.fit(data)
-        >>> model.get_cpds()
-        [<TabularCPD representing P(A:2) at 0x7fb98a7d50f0>,
-        <TabularCPD representing P(B:2) at 0x7fb98a7d5588>,
-        <TabularCPD representing P(C:2 | A:2, B:2) at 0x7fb98a7b1f98>]
-        """
-        from pgmpy.estimators import BaseEstimator, MaximumLikelihoodEstimator
-
-        if estimator is None:
-            estimator = MaximumLikelihoodEstimator
-        else:
-            if not issubclass(estimator, BaseEstimator):
-                raise TypeError("Estimator object should be a valid pgmpy estimator.")
-
-        _estimator = estimator(
-            self,
-            data,
-            state_names=state_names,
-        )
-        cpds_list = _estimator.get_parameters(n_jobs=n_jobs, **kwargs)
-        self.add_cpds(*cpds_list)
-
     def fit_update(self, data, n_prev_samples=None, n_jobs=1):
         """
         Method to update the parameters of the BayesianNetwork with more data.


Test patch:
diff --git a/pgmpy/tests/test_base/test_DAG.py b/pgmpy/tests/test_base/test_DAG.py
index 002139128..38d26a0d8 100644
--- a/pgmpy/tests/test_base/test_DAG.py
+++ b/pgmpy/tests/test_base/test_DAG.py
@@ -4,9 +4,16 @@
 
 import networkx as nx
 import numpy as np
+import pandas as pd
 
 import pgmpy.tests.help_functions as hf
 from pgmpy.base import DAG, PDAG
+from pgmpy.estimators import (
+    BayesianEstimator,
+    MaximumLikelihoodEstimator,
+    ExpectationMaximization,
+)
+from pgmpy.factors.discrete import TabularCPD
 
 
 class TestDAGCreation(unittest.TestCase):
@@ -292,6 +299,43 @@ def test_random_dag(self):
 
         dag_latents = DAG.get_random(n_nodes=n_nodes, edge_prob=0.5, latents=True)
 
+    def test_dag_fit(self):
+        self.model = DAG([("A", "C"), ("B", "C")])
+        self.data = pd.DataFrame(data={"A": [0, 0, 1], "B": [0, 1, 0], "C": [1, 1, 0]})
+        self.pseudo_counts = {
+            "A": [[9], [3]],
+            "B": [[9], [3]],
+            "C": [[9, 9, 9, 9], [3, 3, 3, 3]],
+        }
+
+        self.fitted_model_bayesian = self.model.fit(
+            self.data,
+            estimator=BayesianEstimator,
+            prior_type="dirichlet",
+            pseudo_counts=self.pseudo_counts,
+        )
+
+        self.fitted_model_mle = self.model.fit(
+            self.data, estimator=MaximumLikelihoodEstimator
+        )
+
+        self.fitted_model_em = self.model.fit(
+            self.data, estimator=ExpectationMaximization
+        )
+
+        self.assertEqual(
+            self.fitted_model_bayesian.get_cpds("B"),
+            TabularCPD("B", 2, [[11.0 / 15], [4.0 / 15]]),
+        )
+        self.assertEqual(
+            self.fitted_model_mle.get_cpds("B"),
+            TabularCPD("B", 2, [[2.0 / 3], [1.0 / 3]]),
+        )
+        self.assertEqual(
+            self.fitted_model_em.get_cpds("B"),
+            TabularCPD("B", 2, [[2.0 / 3], [1.0 / 3]]),
+        )
+
     def tearDown(self):
         del self.graph
 
diff --git a/pgmpy/tests/test_estimators/test_BayesianEstimator.py b/pgmpy/tests/test_estimators/test_BayesianEstimator.py
index 96a78c907..4704d3b85 100644
--- a/pgmpy/tests/test_estimators/test_BayesianEstimator.py
+++ b/pgmpy/tests/test_estimators/test_BayesianEstimator.py
@@ -8,12 +8,14 @@
 from pgmpy.estimators import BayesianEstimator
 from pgmpy.factors.discrete import TabularCPD
 from pgmpy.models import BayesianNetwork
+from pgmpy.base import DAG
 
 
 class TestBayesianEstimator(unittest.TestCase):
     def setUp(self):
         self.m1 = BayesianNetwork([("A", "C"), ("B", "C")])
         self.model_latent = BayesianNetwork([("A", "C"), ("B", "C")], latents=["C"])
+        self.dag_with_latents = DAG([("A", "B")], latents=["C"])
         self.d1 = pd.DataFrame(data={"A": [0, 0, 1], "B": [0, 1, 0], "C": [1, 1, 0]})
         self.d2 = pd.DataFrame(
             data={
@@ -30,6 +32,7 @@ def setUp(self):
 
     def test_error_latent_model(self):
         self.assertRaises(ValueError, BayesianEstimator, self.model_latent, self.d1)
+        self.assertRaises(ValueError, BayesianEstimator, self.dag_with_latents, self.d1)
 
     def test_estimate_cpd_dirichlet(self):
         cpd_A = self.est1.estimate_cpd(
diff --git a/pgmpy/tests/test_metrics/test_metrics.py b/pgmpy/tests/test_metrics/test_metrics.py
index f4bdc0608..74519c9ee 100644
--- a/pgmpy/tests/test_metrics/test_metrics.py
+++ b/pgmpy/tests/test_metrics/test_metrics.py
@@ -7,7 +7,14 @@
 
 from pgmpy.base import DAG
 from pgmpy.estimators.CITests import chi_square
-from pgmpy.metrics import *
+from pgmpy.metrics import (
+    SHD,
+    correlation_score,
+    fisher_c,
+    implied_cis,
+    log_likelihood_score,
+    structure_score,
+)
 from pgmpy.models import BayesianNetwork
 from pgmpy.utils import get_example_model
 
@@ -174,3 +181,34 @@ def test_fisher_c(self):
 
         p_value = fisher_c(self.model_alarm_random, self.df_alarm, chi_square)
         self.assertEqual(p_value, 0)
+
+
+class TestStructuralHammingDistance(unittest.TestCase):
+    def setUp(self):
+        self.dag_1 = BayesianNetwork([(1, 2)])
+        self.dag_2 = BayesianNetwork([(2, 1)])
+
+        self.dag_3 = BayesianNetwork([(1, 2), (2, 4), (1, 3), (3, 4)])
+        self.dag_4 = BayesianNetwork([(1, 2), (1, 3), (3, 2), (3, 4)])
+
+        self.dag_5 = BayesianNetwork([(1, 2), (1, 3), (3, 2), (3, 5)])
+
+        self.large_dag_1 = BayesianNetwork(
+            [(1, 2), (1, 3), (2, 4), (3, 5), (4, 5), (5, 6)]
+        )
+        self.large_dag_2 = BayesianNetwork(
+            [(1, 2), (1, 3), (4, 2), (3, 5), (4, 6), (5, 6)]
+        )
+
+    def test_shd(self):
+        self.assertEqual(SHD(self.dag_1, self.dag_2), 1)
+
+    def test_shd(self):
+        self.assertEqual(SHD(self.dag_3, self.dag_4), 2)
+
+    def test_shd(self):
+        self.assertEqual(SHD(self.large_dag_1, self.large_dag_2), 3)
+
+    def test_shd_unequal_graphs(self):
+        with self.assertRaises(ValueError, msg="The graphs must have the same nodes."):
+            SHD(self.dag_4, self.dag_5)