Repository: pgmpy/pgmpy

Issue:
Add MU+BDs for Bayesian network structure learning
Add in the ability to use the Marginal Uniform and Bayesian Dirichlet sparse prior from http://proceedings.mlr.press/v52/scutari16.pdf for structure learning in Bayesian networks. 

I think this will be fairly straightforward by adding in a class that inherits from the BDeuScore class in `pgmpy.estimators.StructureScore`.


Discussion:
ankurankan: @mjmt05 Sorry for so much delay on this. I went through the PR but if I am understanding things correctly, there seems to be a problem with it right now. So, if you look at the `HillClimbSearch` class for structure learning, it uses the `local_score` method from the scoring classes to compute the change in score for any modification that it does in the network. But with the current implementation of BDsScore, the structure learning will get local scores same as BDeUScore as the `local_score` method doesn't use the `structure_prior`. I think there are multiple ways to fix it. The easiest one would be to implement `local_score` for BDsScore as well which computes the change in structure_prior as well. What do you think? 
mjmt05: Oh my bad, for some reason I thought it was calling the score method from the base StructureScore class. That sounds good if somewhat repetitive. In looking at this I also just noticed that we'll have to add in bdsscore to the list of supported_methods in HillClimbSearch.py
ankurankan: @mjmt05 From my understanding the code shouldn't be repetitive. We could call the `BDeUScore.local_score` from `BDsScore.local_score` and then manipulate the returned score to also include the structure prior?
mjmt05: Oh I see what you're saying, yep makes sense!

Patch:
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index d65b6fe9..00000000
--- a/.travis.yml
+++ /dev/null
@@ -1,209 +0,0 @@
-env:
-    global:
-        - secure: "gESRFSRjxdrv8EebgIhIutOHKowghOlNn9UKnNEulyAv3fTT26++dQW3kIP/TlKmIA8ClSLMfegV69XCYqjbSX6fZtE31g1VxGY6ZKt3wreKpHzLru4qN0/WoOtesga9i0pocFjBSl4bxFECJSqExCsMB+Yl9Gm7H6TO3597OP4=" 
-        - secure: "lsV2jOGu759lCZso9/SetSnCR/IzKZxhLzfgTxbMFnAhwat0dGGND2/OmzZsCqD9NeP4h2kO06MkFPmqcHLXALVucNA1rzK8sBrqCd3ef8sc2eM6bDODwHtnfy3/Xf+W1nzv6gONoVjON1SnnPSwtGkeVQGUw6k//5sVTGxPNO0=" 
-        - secure: "HVM08/B0+F3o81Xq3ay/RojRwOpSDvBLdEUMZFuX+510e87gIFtFH+fiZgsWFSkpZzQFT9/QrQaBioDIZD7fHiTOxypVHJPxqiEq5I8ktBrGNREXYHdQeFxMg+SproOQAaA0tjG6lB54OFiPekvIB3ZtfKeYe4OXeDz48Bg5DyI=" 
-
-language: python
-sudo: false
-
-cache: 
-    apt: true
-    directories:
-        - $HOME/.pip-cache
-        - $HOME/gpg_install_prefix
-        - $HOME/.cache/pip
-        - $HOME/download
-
-python:
-  - "3.6"
-  - "3.7"
-  - "3.8"
-  - "3.9"
-
-before_install:
-  - pip install pip -U
-  # run a linter before doing anything else
-  - pip install black -U
-  - black --diff .
-  - black --check .
-  # install the rest of the external requirements
-  - pip install -r requirements.txt -U
-
-install:
-  # Install the library (note we could do this in non-develop mode instead)
-  - travis_retry pip install -e .
-    
-script: 
-  # Execute the tests
-  - travis_wait ./run_tests.py 
-
-after_success: 
-    # Upload coverage results
-    - codecov 
-    # Prep GPG and encryption utilities
-    - gpg --version
-    - gpg2 --version
-    - export GPG_EXECUTABLE=gpg2
-    - openssl version
-    - | 
-        __heredoc__='''
-        # THE SCRIPT IN THIS HEREDOC IS USED TO SETUP SECRET REPO VARIABLES
-        #
-        # Load or generate secrets
-        source $(secret_loader.sh)
-        echo $TWINE_USERNAME
-        echo $TWINE_PASSWORD
-        echo $CI_GITHUB_SECRET
-
-        # encrypt relevant travis variables
-        #  sudo apt install ruby ruby-dev -y
-        #  sudo gem install travis
-        travis encrypt TWINE_USERNAME=$TWINE_USERNAME 
-        travis encrypt TWINE_PASSWORD=$TWINE_PASSWORD  
-        travis encrypt CI_GITHUB_SECRET=$CI_GITHUB_SECRET
-
-        # HOW TO ENCRYPT YOUR SECRET GPG KEY
-        IDENTIFIER="<YOUR-GPG-IDENTIFIER> 
-        KEYID=$(gpg --list-keys --keyid-format LONG "$IDENTIFIER" | head -n 2 | tail -n 1 | awk '{print $1}' | tail -c 9)
-        echo "KEYID = $KEYID"
-
-        # Export plaintext gpg public keys, private keys, and trust info
-        mkdir -p dev
-        gpg --armor --export-secret-keys $KEYID > dev/travis_secret_gpg_key.pgp
-        gpg --armor --export $KEYID > dev/travis_public_gpg_key.pgp
-        gpg --export-ownertrust > dev/gpg_owner_trust
-
-        # Encrypt gpg keys and trust with CI secret
-        TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -e -a -in dev/travis_public_gpg_key.pgp > dev/travis_public_gpg_key.pgp.enc
-        TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -e -a -in dev/travis_secret_gpg_key.pgp > dev/travis_secret_gpg_key.pgp.enc
-        TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -e -a -in dev/gpg_owner_trust > dev/gpg_owner_trust.enc
-        echo $KEYID > dev/public_gpg_key
-
-        source $(secret_unloader.sh)
-
-        # Look at what we did, clean up, and add it to git
-        ls dev/*.enc
-        rm dev/gpg_owner_trust dev/*.pgp
-        git status
-        git add dev/*.enc
-        git add dev/public_gpg_key
-
-        '''  # <hack vim "regex" parser> '
-    - |
-        # Install a more recent version of GPG
-        # https://gnupg.org/download/
-        export GPG_INSTALL_PREFIX=$HOME/gpg_install_prefix
-        export LD_LIBRARY_PATH=$GPG_INSTALL_PREFIX/lib:$LD_LIBRARY_PATH
-        export PATH=$GPG_INSTALL_PREFIX/bin:$PATH
-        export CPATH=$GPG_INSTALL_PREFIX/include:$CPATH
-        export GPG_EXECUTABLE=$GPG_INSTALL_PREFIX/bin/gpg
-        ls $GPG_INSTALL_PREFIX
-        ls $GPG_INSTALL_PREFIX/bin || echo "no bin"
-        # try and have travis cache this
-        if [[ ! -f "$GPG_INSTALL_PREFIX/bin/gpg" ]]; then
-            # This part of the script installs a newer version of GPG on the CI
-            # machine so we can sign our releases.
-            mkdir -p $GPG_INSTALL_PREFIX
-            echo $GPG_INSTALL_PREFIX
-            OLD=$(pwd)
-            cd $GPG_INSTALL_PREFIX
-            pip install ubelt
-
-            ERROR_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/libgpg-error/libgpg-error-1.36.tar.bz2',
-                hash_prefix='6e5f853f77dc04f0091d94b224cab8e669042450f271b78d0ea0219',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            GCRYPT_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/libgcrypt/libgcrypt-1.8.5.tar.bz2',
-                hash_prefix='b55e16e838d1b1208e7673366971ae7c0f9c1c79e042f41c03d1',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            KSBA_CRYPT_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/libksba/libksba-1.3.5.tar.bz2',
-                hash_prefix='60179bfd109b7b4fd8d2b30a3216540f03f5a13620d9a5b63f1f95',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            ASSUAN_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/libassuan/libassuan-2.5.3.tar.bz2',
-                hash_prefix='e7ccb651ea75b07b2e687d48d86d0ab83cba8e2af7f30da2aec',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            NTBLTLS_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/ntbtls/ntbtls-0.1.2.tar.bz2',
-                hash_prefix='54468208359dc88155b14cba37773984d7d6f0f37c7a4ce13868d',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            NPTH_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/npth/npth-1.6.tar.bz2',
-                hash_prefix='2ed1012e14a9d10665420b9a23628be7e206fd9348111ec751349b',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-            GPG_FPATH=$(python -c "import ubelt as ub; print(ub.grabdata(
-                'https://gnupg.org/ftp/gcrypt/gnupg/gnupg-2.2.17.tar.bz2',
-                hash_prefix='a3cd094addac62b4b4ec1683005a2bec761ea2aacf6daf904316b',
-                dpath=ub.ensuredir('$HOME/.pip-cache'), verbose=0))")
-
-            tar xjf $ERROR_FPATH
-            tar xjf $GCRYPT_FPATH
-            tar xjf $KSBA_CRYPT_FPATH
-            tar xjf $ASSUAN_FPATH
-            tar xjf $NTBLTLS_FPATH
-            tar xjf $NPTH_FPATH
-            tar xjf $GPG_FPATH
-            (cd libgpg-error-1.36 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd libgcrypt-1.8.5 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd libksba-1.3.5 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd libassuan-2.5.3 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd ntbtls-0.1.2 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd npth-1.6 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            (cd gnupg-2.2.17 && ./configure --prefix=$GPG_INSTALL_PREFIX && make install)
-            echo "GPG_EXECUTABLE = '$GPG_EXECUTABLE'"
-            cd $OLD
-        fi
-    # Decrypt and import GPG Keys / trust
-    - $GPG_EXECUTABLE --version
-    - openssl version
-    - $GPG_EXECUTABLE --list-keys
-    - TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -d -a -in dev/travis_public_gpg_key.pgp.enc | $GPG_EXECUTABLE --import 
-    - TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -d -a -in dev/gpg_owner_trust.enc | $GPG_EXECUTABLE --import-ownertrust
-    - TSP=$CI_GITHUB_SECRET openssl enc -aes-256-cbc -md MD5 -pass env:TSP -d -a -in dev/travis_secret_gpg_key.pgp.enc | $GPG_EXECUTABLE --import 
-    - $GPG_EXECUTABLE --list-keys
-    - MB_PYTHON_TAG=$(python -c "import setup; print(setup.MB_PYTHON_TAG)")
-    - VERSION=$(python -c "import setup; print(setup.VERSION)") 
-    - |
-        pip install twine
-        if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then
-          pip install six pyopenssl ndg-httpsclient pyasn1 -U --user
-          pip install requests[security] twine --user
-        elfi
-        if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
-          pip install six twine
-          pip install --upgrade pyOpenSSL
-        fi
-    # Package and publish to pypi (if on release)
-    - |
-        echo "TRAVIS_BRANCH = $TRAVIS_BRANCH"
-
-        KEYID=$(cat dev/public_gpg_key)
-        echo "KEYID = '$KEYID'"
-        if [[ "$TRAVIS_BRANCH" == "release" ]]; then
-            export CURRENT_BRANCH=$TRAVIS_BRANCH
-            TAG_AND_UPLOAD=yes
-        else
-            TAG_AND_UPLOAD=no
-        fi
-        MB_PYTHON_TAG=$MB_PYTHON_TAG \
-            USE_GPG=True \
-            GPG_KEYID=$KEYID \
-            CURRENT_BRANCH=$TRAVIS_BRANCH \
-            TWINE_PASSWORD=$TWINE_PASSWORD \
-            TWINE_USERNAME=$TWINE_USERNAME \
-            GPG_EXECUTABLE=$GPG_EXECUTABLE \
-            DEPLOY_BRANCH=release \
-            TAG_AND_UPLOAD=$TAG_AND_UPLOAD \
-            ./publish.sh
-
-notifications:
-    email: false
-    webhooks:
-        urls:
-            - https://webhooks.gitter.im/e/e2ce04225e5e5d1ab64c
-        on_success: change # options: [always|never|change] default: always
-        on_failure: change # options: [always|never|change] default: always
-        on_start: false     # default: false
diff --git a/README.md b/README.md
index 863ce22a..db1eb75e 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,6 @@
 pgmpy
 =====
-[![Build Status](https://api.travis-ci.com/pgmpy/pgmpy.svg?branch=dev)](https://travis-ci.com/pgmpy/pgmpy)
-[![Appveyor](https://ci.appveyor.com/api/projects/status/github/pgmpy/pgmpy?branch=dev)](https://www.appveyor.com/)
+![Build](https://github.com/pgmpy/pgmpy/actions/workflows/ci.yml/badge.svg?branch=dev)
 [![codecov](https://codecov.io/gh/pgmpy/pgmpy/branch/dev/graph/badge.svg)](https://codecov.io/gh/pgmpy/pgmpy)
 [![Codacy Badge](https://api.codacy.com/project/badge/Grade/78a8256c90654c6892627f6d8bbcea14)](https://www.codacy.com/gh/pgmpy/pgmpy?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=pgmpy/pgmpy&amp;utm_campaign=Badge_Grade)
 [![Downloads](https://img.shields.io/pypi/dm/pgmpy.svg)](https://pypistats.org/packages/pgmpy)
diff --git a/appveyor.yml b/appveyor.yml
deleted file mode 100644
index 50bb79b9..00000000
--- a/appveyor.yml
+++ /dev/null
@@ -1,32 +0,0 @@
-environment:
-  matrix:
-    # pytorch only seems to work on 64bit Python versions on Windows
-    - PYTHON: "C:\\Python39-x64"
-      PYTHON_VERSION: "3.9.x"
-      APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
-      PYTHON_ARCH: "64"
-    - PYTHON: "C:\\Python38-x64" 
-    - PYTHON: "C:\\Python37-x64"
-    - PYTHON: "C:\\Python36-x64"
-
-init:
-  - "ECHO %PYTHON% %PYTHON_VERSION% %PYTHON_ARCH%"
-
-install:
-  - "%PYTHON%\\python.exe -m pip install pip -U"
-  - "%PYTHON%\\python.exe -m pip install torch==1.8.1+cpu -f https://download.pytorch.org/whl/torch_stable.html"
-  - "%PYTHON%\\python.exe -m pip install -r requirements.txt -U"
-  - "%PYTHON%\\python.exe -m pip install -e ."
-
-build: off
-
-test_script:
-  - "%PYTHON%\\python.exe -m pytest"
-
-notifications:
-  - provider: Webhook
-    url: https://webhooks.gitter.im/e/c084306bedf9aa508af2
-    method: GET
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: true
diff --git a/docs/estimators.rst b/docs/estimators.rst
index 8e14a68b..21de5ce2 100644
--- a/docs/estimators.rst
+++ b/docs/estimators.rst
@@ -19,6 +19,12 @@ Bayesian Estimator
 .. autoclass:: pgmpy.estimators.BayesianEstimator
    :members:
 
+Expectation Maximization (EM)
+-----------------------------
+
+.. autoclass:: pgmpy.estimators.ExpectationMaximization
+   :members:
+
 BDeu Score
 ----------
 
diff --git a/pgmpy/estimators/BayesianEstimator.py b/pgmpy/estimators/BayesianEstimator.py
index 8af808ca..f1b42029 100644
--- a/pgmpy/estimators/BayesianEstimator.py
+++ b/pgmpy/estimators/BayesianEstimator.py
@@ -163,7 +163,11 @@ class BayesianEstimator(ParameterEstimator):
         # Throw a warning if pseudo_count is specified without prior_type=dirichlet
         #     cast to np.array first to use the array.size attribute, which returns 0 also for [[],[]]
         #     (where len([[],[]]) evaluates to 2)
-        if np.array(pseudo_counts).size > 0 and (prior_type != "dirichlet"):
+        if (
+            pseudo_counts is not None
+            and np.array(pseudo_counts).size > 0
+            and (prior_type != "dirichlet")
+        ):
             warn(
                 f"pseudo count specified with {prior_type} prior. It will be ignored, use dirichlet prior for specifying pseudo_counts"
             )
diff --git a/pgmpy/estimators/EM.py b/pgmpy/estimators/EM.py
new file mode 100644
index 00000000..735ff082
--- /dev/null
+++ b/pgmpy/estimators/EM.py
@@ -0,0 +1,230 @@
+import warnings
+from itertools import product, chain
+
+import numpy as np
+import pandas as pd
+from tqdm import tqdm
+
+from pgmpy.estimators import ParameterEstimator, MaximumLikelihoodEstimator
+from pgmpy.models import BayesianModel
+from pgmpy.factors.discrete import TabularCPD
+from pgmpy.global_vars import SHOW_PROGRESS
+
+
+class ExpectationMaximization(ParameterEstimator):
+    def __init__(self, model, data, **kwargs):
+        """
+        Class used to compute parameters for a model using Expectation
+        Maximization (EM).
+
+        EM is an iterative algorithm commonly used for
+        estimation in the case when there are latent variables in the model.
+        The algorithm iteratively improves the parameter estimates maximizing
+        the likelihood of the given data.
+
+        Parameters
+        ----------
+        model: A pgmpy.models.BayesianModel instance
+
+        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`.)
+
+        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.
+
+        complete_samples_only: bool (optional, default `True`)
+            Specifies how to deal with missing data, if present. If set to
+            `True` all rows that contain `np.NaN` somewhere are ignored. If
+            `False` then, for each variable, every row where neither the
+            variable nor its parents are `np.NaN` is used.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> import pandas as pd
+        >>> from pgmpy.models import BayesianModel
+        >>> from pgmpy.estimators import ExpectationMaximization
+        >>> data = pd.DataFrame(np.random.randint(low=0, high=2, size=(1000, 5)),
+        ...                       columns=['A', 'B', 'C', 'D', 'E'])
+        >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D'), ('B', 'E')])
+        >>> estimator = ExpectationMaximization(model, data)
+        """
+        if not isinstance(model, BayesianModel):
+            raise NotImplementedError(
+                "Expectation Maximization is only implemented for BayesianModel"
+            )
+
+        super(ExpectationMaximization, self).__init__(model, data, **kwargs)
+        self.model_copy = self.model.copy()
+
+    def _get_likelihood(self, datapoint):
+        """
+        Computes the likelihood of a given datapoint. Goes through each
+        CPD matching the combination of states to get the value and multiplies
+        them together.
+        """
+        likelihood = 1
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            for cpd in self.model_copy.cpds:
+                scope = set(cpd.scope())
+                likelihood *= cpd.get_value(
+                    **{key: value for key, value in datapoint.items() if key in scope}
+                )
+        return likelihood
+
+    def _compute_weights(self, latent_card):
+        """
+        For each data point, creates extra data points for each possible combination
+        of states of latent variables and assigns weights to each of them.
+        """
+        cache = []
+
+        data_unique = self.data.drop_duplicates()
+        n_counts = self.data.groupby(list(self.data.columns)).size().to_dict()
+
+        for i in range(data_unique.shape[0]):
+            v = list(product(*[range(card) for card in latent_card.values()]))
+            latent_combinations = np.array(v, dtype=int)
+            df = data_unique.iloc[[i] * latent_combinations.shape[0]].reset_index(
+                drop=True
+            )
+            for index, latent_var in enumerate(latent_card.keys()):
+                df[latent_var] = latent_combinations[:, index]
+
+            weights = df.apply(lambda t: self._get_likelihood(dict(t)), axis=1)
+            df["_weight"] = (weights / weights.sum()) * n_counts[
+                tuple(data_unique.iloc[i])
+            ]
+            cache.append(df)
+
+        return pd.concat(cache, copy=False)
+
+    def _is_converged(self, new_cpds, atol=1e-08):
+        """
+        Checks if the values of `new_cpds` is within tolerance limits of current
+        model cpds.
+        """
+        for cpd in new_cpds:
+            if not cpd.__eq__(self.model_copy.get_cpds(node=cpd.scope()[0]), atol=atol):
+                return False
+        return True
+
+    def get_parameters(
+        self,
+        latent_card=None,
+        max_iter=100,
+        atol=1e-08,
+        n_jobs=-1,
+        seed=None,
+        show_progress=True,
+    ):
+        """
+        Method to estimate all model parameters (CPDs) using Expecation Maximization.
+
+        Parameters
+        ----------
+        latent_card: dict (default: None)
+            A dictionary of the form {latent_var: cardinality} specifying the
+            cardinality (number of states) of each latent variable. If None,
+            assumes `2` states for each latent variable.
+
+        max_iter: int (default: 100)
+            The maximum number of iterations the algorithm is allowed to run for.
+            If max_iter is reached, return the last value of parameters.
+
+        atol: int (default: 1e-08)
+            The absolute accepted tolerance for checking convergence. If the parameters
+            change is less than atol in an iteration, the algorithm will exit.
+
+        n_jobs: int (default: -1)
+            Number of jobs to run in parallel. Default: -1 uses all the processors.
+
+        seed: int
+            The random seed to use for generating the intial values.
+
+        show_progress: boolean (default: True)
+            Whether to show a progress bar for iterations.
+
+        Returns
+        -------
+        list: A list of estimated CPDs for the model.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> import pandas as pd
+        >>> from pgmpy.models import BayesianModel
+        >>> from pgmpy.estimators import ExpectationMaximization as EM
+        >>> data = pd.DataFrame(np.random.randint(low=0, high=2, size=(1000, 3)),
+        ...                       columns=['A', 'C', 'D'])
+        >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D')], latents={'B'})
+        >>> estimator = EM(model, data)
+        >>> estimator.get_parameters(latent_card={'B': 3})
+        [<TabularCPD representing P(C:2) at 0x7f7b534251d0>,
+        <TabularCPD representing P(B:3 | C:2, A:2) at 0x7f7b4dfd4da0>,
+        <TabularCPD representing P(A:2) at 0x7f7b4dfd4fd0>,
+        <TabularCPD representing P(D:2 | C:2) at 0x7f7b4df822b0>]
+        """
+        # Step 1: Parameter checks
+        if latent_card is None:
+            latent_card = {var: 2 for var in self.model_copy.latents}
+
+        # Step 2: Create structures/variables to be used later.
+        n_states_dict = {key: len(value) for key, value in self.state_names.items()}
+        n_states_dict.update(latent_card)
+        for var in self.model_copy.latents:
+            self.state_names[var] = list(range(n_states_dict[var]))
+
+        # Step 3: Initialize random CPDs if starting values aren't provided.
+        if seed is not None:
+            np.random.seed(seed)
+
+        cpds = []
+        for node in self.model_copy.nodes():
+            parents = list(self.model_copy.predecessors(node))
+            cpds.append(
+                TabularCPD.get_random(
+                    variable=node,
+                    evidence=parents,
+                    cardinality={
+                        var: n_states_dict[var] for var in chain([node], parents)
+                    },
+                    state_names={
+                        var: self.state_names[var] for var in chain([node], parents)
+                    },
+                )
+            )
+
+        self.model_copy.add_cpds(*cpds)
+
+        if show_progress and SHOW_PROGRESS:
+            pbar = tqdm(total=max_iter)
+
+        # Step 4: Run the EM algorithm.
+        for _ in range(max_iter):
+            # Step 4.1: E-step: Expands the dataset and computes the likelihood of each
+            #           possible state of latent variables.
+            weighted_data = self._compute_weights(latent_card)
+            # Step 4.2: M-step: Uses the weights of the dataset to do a weighted MLE.
+            new_cpds = MaximumLikelihoodEstimator(
+                self.model_copy, weighted_data
+            ).get_parameters(n_jobs=n_jobs, weighted=True)
+
+            # Step 4.3: Check of convergence and max_iter
+            if self._is_converged(new_cpds, atol=atol):
+                if show_progress and SHOW_PROGRESS:
+                    pbar.close()
+                return new_cpds
+
+            else:
+                self.model_copy.cpds = new_cpds
+                if show_progress and SHOW_PROGRESS:
+                    pbar.update(1)
+
+        return cpds
diff --git a/pgmpy/estimators/HillClimbSearch.py b/pgmpy/estimators/HillClimbSearch.py
index 07473ce6..f8141d04 100644
--- a/pgmpy/estimators/HillClimbSearch.py
+++ b/pgmpy/estimators/HillClimbSearch.py
@@ -11,6 +11,7 @@ from pgmpy.estimators import (
     K2Score,
     ScoreCache,
     BDeuScore,
+    BDsScore,
     BicScore,
 )
 from pgmpy.base import DAG
@@ -56,7 +57,15 @@ class HillClimbSearch(StructureEstimator):
         super(HillClimbSearch, self).__init__(data, **kwargs)
 
     def _legal_operations(
-        self, model, score, tabu_list, max_indegree, black_list, white_list, fixed_edges
+        self,
+        model,
+        score,
+        structure_score,
+        tabu_list,
+        max_indegree,
+        black_list,
+        white_list,
+        fixed_edges,
     ):
         """Generates a list of legal (= not in tabu_list) graph modifications
         for a given model, together with their score changes. Possible graph modifications:
@@ -90,6 +99,7 @@ class HillClimbSearch(StructureEstimator):
                     new_parents = old_parents + [X]
                     if len(new_parents) <= max_indegree:
                         score_delta = score(Y, new_parents) - score(Y, old_parents)
+                        score_delta += structure_score("+")
                         yield (operation, score_delta)
 
         # Step 2: Get all legal operations for removing edges
@@ -100,6 +110,7 @@ class HillClimbSearch(StructureEstimator):
                 new_parents = old_parents[:]
                 new_parents.remove(X)
                 score_delta = score(Y, new_parents) - score(Y, old_parents)
+                score_delta += structure_score("-")
                 yield (operation, score_delta)
 
         # Step 3: Get all legal operations for flipping edges
@@ -127,6 +138,7 @@ class HillClimbSearch(StructureEstimator):
                             - score(X, old_X_parents)
                             - score(Y, old_Y_parents)
                         )
+                        score_delta += structure_score("flip")
                         yield (operation, score_delta)
 
     def estimate(
@@ -153,7 +165,7 @@ class HillClimbSearch(StructureEstimator):
         ----------
         scoring_method: str or StructureScore instance
             The score to be optimized during structure estimation.  Supported
-            structure scores: k2score, bdeuscore, bicscore. Also accepts a
+            structure scores: k2score, bdeuscore, bdsscore, bicscore. Also accepts a
             custom score but it should be an instance of `StructureScore`.
 
         start_dag: DAG instance
@@ -219,6 +231,7 @@ class HillClimbSearch(StructureEstimator):
         supported_methods = {
             "k2score": K2Score,
             "bdeuscore": BDeuScore,
+            "bdsscore": BDsScore,
             "bicscore": BicScore,
         }
         if (
@@ -228,7 +241,7 @@ class HillClimbSearch(StructureEstimator):
             )
         ) and (not isinstance(scoring_method, StructureScore)):
             raise ValueError(
-                "scoring_method should either be one of k2score, bdeuscore, bicscore, or an instance of StructureScore"
+                "scoring_method should either be one of k2score, bdeuscore, bicscore, bdsscore, or an instance of StructureScore"
             )
 
         if isinstance(scoring_method, str):
@@ -291,6 +304,7 @@ class HillClimbSearch(StructureEstimator):
                 self._legal_operations(
                     current_model,
                     score_fn,
+                    score.structure_prior_ratio,
                     tabu_list,
                     max_indegree,
                     black_list,
diff --git a/pgmpy/estimators/MLE.py b/pgmpy/estimators/MLE.py
index 7fd5ef1d..8c3dbace 100644
--- a/pgmpy/estimators/MLE.py
+++ b/pgmpy/estimators/MLE.py
@@ -45,21 +45,31 @@ class MaximumLikelihoodEstimator(ParameterEstimator):
         >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D'), ('B', 'E')])
         >>> estimator = MaximumLikelihoodEstimator(model, data)
         """
-
         if not isinstance(model, BayesianModel):
             raise NotImplementedError(
                 "Maximum Likelihood Estimate is only implemented for BayesianModel"
             )
-        elif len(model.latents) != 0:
+        elif set(model.nodes()) > set(data.columns):
             raise ValueError(
-                f"Maximum Likelihood Estimator works only for models with all observed variables. Found latent variables: {model.latents}"
+                f"Maximum Likelihood Estimator works only for models with all observed variables. Found latent variables: {model.latents}."
             )
 
         super(MaximumLikelihoodEstimator, self).__init__(model, data, **kwargs)
 
-    def get_parameters(self, n_jobs=-1):
+    def get_parameters(self, n_jobs=-1, weighted=False):
         """
-        Method to estimate the model parameters (CPDs) using Maximum Likelihood Estimation.
+        Method to estimate the model parameters (CPDs) using Maximum Likelihood
+        Estimation.
+
+        Parameters
+        ----------
+        n_jobs: int (default: -1)
+            Number of jobs to run in parallel. Default: -1 uses all the processors.
+
+        weighted: bool
+            If weighted=True, the data must contain a `_weight` column specifying the
+            weight of each datapoint (row). If False, assigns an equal weight to each
+            datapoint.
 
         Returns
         -------
@@ -87,12 +97,12 @@ class MaximumLikelihoodEstimator(ParameterEstimator):
         """
 
         parameters = Parallel(n_jobs=n_jobs, prefer="threads")(
-            delayed(self.estimate_cpd)(node) for node in self.model.nodes()
+            delayed(self.estimate_cpd)(node, weighted) for node in self.model.nodes()
         )
 
         return parameters
 
-    def estimate_cpd(self, node):
+    def estimate_cpd(self, node, weighted=False):
         """
         Method to estimate the CPD for a given variable.
 
@@ -101,6 +111,11 @@ class MaximumLikelihoodEstimator(ParameterEstimator):
         node: int, string (any hashable python object)
             The name of the variable for which the CPD is to be estimated.
 
+        weighted: bool
+            If weighted=True, the data must contain a `_weight` column specifying the
+            weight of each datapoint (row). If False, assigns an equal weight to each
+            datapoint.
+
         Returns
         -------
         CPD: TabularCPD
@@ -132,7 +147,7 @@ class MaximumLikelihoodEstimator(ParameterEstimator):
         ╘══════╧══════╧══════╧══════╧══════╛
         """
 
-        state_counts = self.state_counts(node)
+        state_counts = self.state_counts(node, weighted=weighted)
 
         # if a column contains only `0`s (no states observed for some configuration
         # of parents' states) fill that column uniformly instead
diff --git a/pgmpy/estimators/StructureScore.py b/pgmpy/estimators/StructureScore.py
index e9fdfd69..30c029fb 100644
--- a/pgmpy/estimators/StructureScore.py
+++ b/pgmpy/estimators/StructureScore.py
@@ -79,6 +79,11 @@ class StructureScore(BaseEstimator):
         """A (log) prior distribution over models. Currently unused (= uniform)."""
         return 0
 
+    def structure_prior_ratio(self, operation):
+        """Return the log ratio of the prior probabilities for a given proposed change to the DAG.
+        Currently unused (=uniform)."""
+        return 0
+
 
 class K2Score(StructureScore):
     def __init__(self, data, **kwargs):
@@ -182,6 +187,9 @@ class BDeuScore(StructureScore):
         self.equivalent_sample_size = equivalent_sample_size
         super(BDeuScore, self).__init__(data, **kwargs)
 
+    def get_number_of_parent_states(self, state_counts):
+        return float(state_counts.shape[1])
+
     def local_score(self, variable, parents):
         'Computes a score that measures how much a \
         given variable is "influenced" by a given list of potential parents.'
@@ -189,14 +197,12 @@ class BDeuScore(StructureScore):
         var_states = self.state_names[variable]
         var_cardinality = len(var_states)
         state_counts = self.state_counts(variable, parents)
-        num_parents_states = float(state_counts.shape[1])
+        num_parents_states = self.get_number_of_parent_states(state_counts)
 
         counts = np.asarray(state_counts)
         log_gamma_counts = np.zeros_like(counts, dtype=np.float_)
-
         alpha = self.equivalent_sample_size / num_parents_states
         beta = self.equivalent_sample_size / counts.size
-
         # Compute log(gamma(counts + beta))
         gammaln(counts + beta, out=log_gamma_counts)
 
@@ -210,7 +216,71 @@ class BDeuScore(StructureScore):
             + num_parents_states * lgamma(alpha)
             - counts.size * lgamma(beta)
         )
+        return score
+
+
+class BDsScore(BDeuScore):
+    def __init__(self, data, equivalent_sample_size=10, **kwargs):
+        """
+        Class for Bayesian structure scoring for BayesianModels with Dirichlet priors.
+        The BDs score is the result of setting all Dirichlet hyperparameters/pseudo_counts to
+        `equivalent_sample_size/modified_variable_cardinality` where for the modified_variable_cardinality
+        only the number of parent configurations where there were observed variable counts are considered.
+        The `score`-method measures how well a model is able to describe the given data set.
+
+        Parameters
+        ----------
+        data: pandas DataFrame object
+            datafame 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`.)
+
+        equivalent_sample_size: int (default: 10)
+            The equivalent/imaginary sample size (of uniform pseudo samples) for the dirichlet
+            hyperparameters.
+            The score is sensitive to this value, runs with different values might be useful.
 
+        state_names: dict (optional)
+            A dict indicating, for each variable, the discrete set of states (or values)
+            that the variable can take. If unspecified, the observed values in the data set
+            are taken to be the only possible states.
+
+        complete_samples_only: bool (optional, default `True`)
+            Specifies how to deal with missing data, if present. If set to `True` all rows
+            that contain `np.Nan` somewhere are ignored. If `False` then, for each variable,
+            every row where neither the variable nor its parents are `np.NaN` is used.
+            This sets the behavior of the `state_count`-method.
+
+        References
+        ---------
+        [1] Scutari, Marco. An Empirical-Bayes Score for Discrete Bayesian Networks.
+        Journal of Machine Learning Research, 2016, pp. 438–48
+
+        """
+        super(BDsScore, self).__init__(data, equivalent_sample_size, **kwargs)
+
+    def get_number_of_parent_states(self, state_counts):
+        return float(len(np.where(state_counts.sum(axis=0) > 0)[0]))
+
+    def structure_prior_ratio(self, operation):
+        """Return the log ratio of the prior probabilities for a given proposed change to
+        the DAG.
+        """
+        if operation == "+":
+            return -log(2.0)
+        if operation == "-":
+            return log(2.0)
+        return 0
+
+    def structure_prior(self, model):
+        """
+        Implements the marginal uniform prior for the graph structure where each arc
+        is independent with the probability of an arc for any two nodes in either direction
+        is 1/4 and the probability of no arc between any two nodes is 1/2."""
+        nedges = float(len(model.edges()))
+        nnodes = float(len(model.nodes()))
+        possible_edges = nnodes * (nnodes - 1) / 2.0
+        score = -(nedges + possible_edges) * log(2.0)
         return score
 
 
diff --git a/pgmpy/estimators/__init__.py b/pgmpy/estimators/__init__.py
index d3b9e749..261215f4 100644
--- a/pgmpy/estimators/__init__.py
+++ b/pgmpy/estimators/__init__.py
@@ -1,13 +1,20 @@
 from pgmpy.estimators.base import BaseEstimator, ParameterEstimator, StructureEstimator
 from pgmpy.estimators.MLE import MaximumLikelihoodEstimator
 from pgmpy.estimators.BayesianEstimator import BayesianEstimator
-from pgmpy.estimators.StructureScore import StructureScore, K2Score, BDeuScore, BicScore
+from pgmpy.estimators.StructureScore import (
+    StructureScore,
+    K2Score,
+    BDeuScore,
+    BDsScore,
+    BicScore,
+)
 from pgmpy.estimators.ExhaustiveSearch import ExhaustiveSearch
 from pgmpy.estimators.HillClimbSearch import HillClimbSearch
 from pgmpy.estimators.TreeSearch import TreeSearch
 from pgmpy.estimators.SEMEstimator import SEMEstimator, IVEstimator
 from pgmpy.estimators.ScoreCache import ScoreCache
 from pgmpy.estimators.MmhcEstimator import MmhcEstimator
+from pgmpy.estimators.EM import ExpectationMaximization
 from pgmpy.estimators.PC import PC
 
 __all__ = [
@@ -22,10 +29,12 @@ __all__ = [
     "StructureScore",
     "K2Score",
     "BDeuScore",
+    "BDsScore",
     "BicScore",
     "ScoreCache",
     "SEMEstimator",
     "IVEstimator",
     "MmhcEstimator",
     "PC",
+    "ExpectationMaximization",
 ]
diff --git a/pgmpy/estimators/base.py b/pgmpy/estimators/base.py
index c2912b78..fdceac4e 100644
--- a/pgmpy/estimators/base.py
+++ b/pgmpy/estimators/base.py
@@ -1,10 +1,7 @@
 #!/usr/bin/env python
-from warnings import warn
 from functools import lru_cache
 
-import numpy as np
 import pandas as pd
-from scipy.stats import chisquare
 
 from pgmpy.utils.decorators import convert_args_tuple
 
@@ -17,7 +14,6 @@ class BaseEstimator(object):
 
         Parameters
         ----------
-
         data: pandas DataFrame object datafame 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`.)
@@ -61,13 +57,15 @@ class BaseEstimator(object):
                         self.state_names[var] = self._collect_state_names(var)
 
     def _collect_state_names(self, variable):
-        "Return a list of states that the variable takes in the data"
+        "Return a list of states that the variable takes in the data."
         states = sorted(list(self.data.loc[:, variable].dropna().unique()))
         return states
 
     @convert_args_tuple
     @lru_cache(maxsize=2048)
-    def state_counts(self, variable, parents=[], complete_samples_only=None):
+    def state_counts(
+        self, variable, parents=[], complete_samples_only=None, weighted=False
+    ):
         """
         Return counts how often each state of 'variable' occurred in the data.
         If a list of parents is provided, counting is done conditionally
@@ -88,6 +86,10 @@ class BaseEstimator(object):
             every row where neither the variable nor its parents are `np.NaN` is used.
             Desired default behavior can be passed to the class constructor.
 
+        weighted: bool
+            If True, data must have a `_weight` column specifying the weight of the
+            datapoint (row). If False, each datapoint has a weight of `1`.
+
         Returns
         -------
         state_counts: pandas.DataFrame
@@ -129,9 +131,16 @@ class BaseEstimator(object):
             else self.data.dropna(subset=[variable] + parents)
         )
 
+        if weighted and ("_weight" not in self.data.columns):
+            raise ValueError("data must contain a `_weight` column if weighted=True")
+
         if not parents:
             # count how often each state of 'variable' occured
-            state_count_data = data.loc[:, variable].value_counts()
+            if weighted:
+                state_count_data = data.groupby([variable]).sum()["_weight"]
+            else:
+                state_count_data = data.loc[:, variable].value_counts()
+
             state_counts = (
                 state_count_data.reindex(self.state_names[variable])
                 .fillna(0)
@@ -141,9 +150,16 @@ class BaseEstimator(object):
         else:
             parents_states = [self.state_names[parent] for parent in parents]
             # count how often each state of 'variable' occured, conditional on parents' states
-            state_count_data = (
-                data.groupby([variable] + parents).size().unstack(parents)
-            )
+            if weighted:
+                state_count_data = (
+                    data.groupby([variable] + parents).sum()["_weight"].unstack(parents)
+                )
+
+            else:
+                state_count_data = (
+                    data.groupby([variable] + parents).size().unstack(parents)
+                )
+
             if not isinstance(state_count_data.columns, pd.MultiIndex):
                 state_count_data.columns = pd.MultiIndex.from_arrays(
                     [state_count_data.columns]
@@ -189,7 +205,7 @@ class ParameterEstimator(BaseEstimator):
             This sets the behavior of the `state_count`-method.
         """
 
-        if not set(model.nodes()) <= set(data.columns.values):
+        if not (set(model.nodes()) - model.latents) <= set(data.columns.values):
             raise ValueError(
                 "variable names of the model must be identical to column names in data"
             )
@@ -197,7 +213,7 @@ class ParameterEstimator(BaseEstimator):
 
         super(ParameterEstimator, self).__init__(data, **kwargs)
 
-    def state_counts(self, variable, **kwargs):
+    def state_counts(self, variable, weighted=False, **kwargs):
         """
         Return counts how often each state of 'variable' occurred in the data.
         If the variable has parents, counting is done conditionally
@@ -243,12 +259,9 @@ class ParameterEstimator(BaseEstimator):
 
         parents = sorted(self.model.get_parents(variable))
         return super(ParameterEstimator, self).state_counts(
-            variable, parents=parents, **kwargs
+            variable, parents=parents, weighted=weighted, **kwargs
         )
 
-    def get_parameters(self):
-        pass
-
 
 class StructureEstimator(BaseEstimator):
     def __init__(self, data=None, independencies=None, **kwargs):
diff --git a/pgmpy/factors/discrete/DiscreteFactor.py b/pgmpy/factors/discrete/DiscreteFactor.py
index 6b07ca04..69a864ee 100644
--- a/pgmpy/factors/discrete/DiscreteFactor.py
+++ b/pgmpy/factors/discrete/DiscreteFactor.py
@@ -936,7 +936,15 @@ class DiscreteFactor(BaseFactor, StateNameMixin):
 
     __div__ = __truediv__
 
-    def __eq__(self, other):
+    def __eq__(self, other, atol=1e-08):
+        """
+        Method for checking if two factors are equal.
+
+        Parameters
+        ----------
+        atol: float
+            The maximum allowed difference in values to be considered equal.
+        """
         if not (isinstance(self, DiscreteFactor) and isinstance(other, DiscreteFactor)):
             return False
 
@@ -974,7 +982,7 @@ class DiscreteFactor(BaseFactor, StateNameMixin):
 
             if phi.values.shape != self.values.shape:
                 return False
-            elif not np.allclose(phi.values, self.values):
+            elif not np.allclose(phi.values, self.values, atol=atol):
                 return False
             elif not all(self.cardinality == phi.cardinality):
                 return False
diff --git a/pgmpy/inference/CausalInference.py b/pgmpy/inference/CausalInference.py
index 418c0311..5a3c1c85 100644
--- a/pgmpy/inference/CausalInference.py
+++ b/pgmpy/inference/CausalInference.py
@@ -379,6 +379,72 @@ class CausalInference(object):
         ]
         return np.mean(ate)
 
+    def get_proper_backdoor_graph(self, X, Y, inplace=False):
+        """
+        Returns a proper backdoor graph for the exposure `X` and outcome `Y`.
+        A proper backdoor graph is a graph which remove the first edge of every
+        proper causal path from `X` to `Y`.
+
+        Parameters
+        ----------
+        X: list (array-like)
+            A list of exposure variables.
+
+        Y: list (array-like)
+            A list of outcome variables
+
+        inplace: boolean
+            If inplace is True, modifies the object itself. Otherwise retuns
+            a modified copy of self.
+        """
+        for var in chain(X, Y):
+            if var not in self.model.nodes():
+                raise ValueError(f"{var} not found in the model.")
+
+        model = self.model if inplace else self.model.copy()
+        edges_to_remove = []
+        for source in X:
+            paths = nx.all_simple_edge_paths(model, source, Y)
+            for path in paths:
+                edges_to_remove.append(path[0])
+        model.remove_edges_from(edges_to_remove)
+        return model
+
+    def is_valid_adjustment_set(self, X, Y, adjustment_set):
+        """
+        Method to test whether `adjustment_set` is a valid adjustment set for
+        identifying the causal effect of `X` on `Y`.
+
+        Parameters
+        ----------
+        X: list (array-like)
+            The set of cause variables.
+
+        Y: list (array-like)
+            The set of predictor variables.
+
+        adjustment_set: list (array-like)
+            The set of variables for which to test whether they satisfy the
+            adjustment set criteria.
+
+        Returns
+        -------
+        boolean: Returns True if `adjustment_set` is a valid adjustment set for
+                identifying the effect of `X` on `Y`. Else returns False.
+
+        Examples
+        --------
+
+        References
+        ----------
+        [1] Perkovic, Emilija, et al. "Complete graphical characterization and construction of adjustment sets in Markov equivalence classes of ancestral graphs." The Journal of Machine Learning Research 18.1 (2017): 8132-8193.
+        """
+        backdoor_graph = self.get_proper_backdoor_graph(X, Y, inplace=False)
+        for (x, y) in zip(X, Y):
+            if backdoor_graph.is_dconnected(start=x, end=y, observed=adjustment_set):
+                return False
+        return True
+
     def query(
         self,
         variables,
diff --git a/pgmpy/inference/base.py b/pgmpy/inference/base.py
index 1b41e520..712322be 100644
--- a/pgmpy/inference/base.py
+++ b/pgmpy/inference/base.py
@@ -143,7 +143,7 @@ class Inference(object):
         # Step 1: Remove all the variables that are d-separated from `variables` when conditioned
         #         on `evidence`
         d_connected = bn.active_trail_nodes(
-            variables=variables, observed=list(evidence.keys())
+            variables=variables, observed=list(evidence.keys()), include_latents=True
         )
         d_connected = set.union(*d_connected.values()).union(evidence.keys())
         bn = bn.subgraph(d_connected)
diff --git a/pgmpy/utils/decorators.py b/pgmpy/utils/decorators.py
index 8257356b..17dcc117 100644
--- a/pgmpy/utils/decorators.py
+++ b/pgmpy/utils/decorators.py
@@ -1,8 +1,8 @@
 def convert_args_tuple(func):
     def _convert_param_to_tuples(
-        obj, variable, parents=tuple(), complete_samples_only=None
+        obj, variable, parents=tuple(), complete_samples_only=None, weighted=False
     ):
         parents = tuple(parents)
-        return func(obj, variable, parents, complete_samples_only)
+        return func(obj, variable, parents, complete_samples_only, weighted)
 
     return _convert_param_to_tuples


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_BDsScore.py b/pgmpy/tests/test_estimators/test_BDsScore.py
new file mode 100644
index 00000000..5f800873
--- /dev/null
+++ b/pgmpy/tests/test_estimators/test_BDsScore.py
@@ -0,0 +1,37 @@
+import unittest
+
+import pandas as pd
+
+from pgmpy.models import BayesianModel
+from pgmpy.estimators import BDsScore
+
+
+class TestBDsScore(unittest.TestCase):
+    def setUp(self):
+        """Example taken from https://arxiv.org/pdf/1708.00689.pdf"""
+        self.d1 = pd.DataFrame(
+            data={
+                "X": [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
+                "Y": [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1],
+                "Z": [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
+                "W": [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1],
+            }
+        )
+        self.m1 = BayesianModel([("W", "X"), ("Z", "X")])
+        self.m1.add_node("Y")
+        self.m2 = BayesianModel([("W", "X"), ("Z", "X"), ("Y", "X")])
+
+    def test_score(self):
+        self.assertAlmostEqual(
+            BDsScore(self.d1, equivalent_sample_size=1).score(self.m1),
+            -36.82311976667139,
+        )
+        self.assertEqual(
+            BDsScore(self.d1, equivalent_sample_size=1).score(self.m2),
+            -45.788991276221964,
+        )
+
+    def tearDown(self):
+        del self.d1
+        del self.m1
+        del self.m2
diff --git a/pgmpy/tests/test_estimators/test_EM.py b/pgmpy/tests/test_estimators/test_EM.py
new file mode 100644
index 00000000..05023b1a
--- /dev/null
+++ b/pgmpy/tests/test_estimators/test_EM.py
@@ -0,0 +1,40 @@
+import unittest
+
+import pandas as pd
+import numpy as np
+
+from pgmpy.models import BayesianModel
+from pgmpy.sampling import BayesianModelSampling
+from pgmpy.estimators import ExpectationMaximization as EM
+from pgmpy.utils import get_example_model
+
+
+class TestEMObserved(unittest.TestCase):
+    def setUp(self):
+        self.model1 = get_example_model("cancer")
+        s = BayesianModelSampling(self.model1)
+        self.data1 = s.forward_sample(int(1e4))
+
+        self.model2 = BayesianModel(self.model1.edges(), latents={"Smoker"})
+        self.model2.add_cpds(*self.model1.cpds)
+        s = BayesianModelSampling(self.model2)
+        self.data2 = s.forward_sample(int(1e4))
+
+    def test_get_parameters(self):
+        est = EM(self.model1, self.data1)
+        cpds = est.get_parameters(seed=42)
+        for est_cpd in cpds:
+            var = est_cpd.variables[0]
+            orig_cpd = self.model1.get_cpds(var)
+            self.assertTrue(orig_cpd.__eq__(est_cpd, atol=0.1))
+
+        est = EM(self.model2, self.data2)
+        cpds = est.get_parameters(seed=42)
+        for est_cpd in cpds:
+            var = est_cpd.variables[0]
+            orig_cpd = self.model2.get_cpds(var)
+
+            if "Smoker" in orig_cpd.variables:
+                orig_cpd.state_names["Smoker"] = [0, 1]
+
+            self.assertTrue(orig_cpd.__eq__(est_cpd, atol=0.3))
diff --git a/pgmpy/tests/test_estimators/test_HillClimbSearch.py b/pgmpy/tests/test_estimators/test_HillClimbSearch.py
index 3aa4a34e..118057fe 100644
--- a/pgmpy/tests/test_estimators/test_HillClimbSearch.py
+++ b/pgmpy/tests/test_estimators/test_HillClimbSearch.py
@@ -14,7 +14,9 @@ class TestHillClimbEstimator(unittest.TestCase):
         )
         self.rand_data["C"] = self.rand_data["B"]
         self.est_rand = HillClimbSearch(self.rand_data)
-        self.score_rand = K2Score(self.rand_data).local_score
+        k2score = K2Score(self.rand_data)
+        self.score_rand = k2score.local_score
+        self.score_structure_prior = k2score.structure_prior_ratio
 
         self.model1 = BayesianModel()
         self.model1.add_nodes_from(["A", "B", "C"])
@@ -47,6 +49,7 @@ class TestHillClimbEstimator(unittest.TestCase):
             self.est_rand._legal_operations(
                 model=self.model2,
                 score=self.score_rand,
+                structure_score=self.score_structure_prior,
                 tabu_list=set(),
                 max_indegree=float("inf"),
                 black_list=set(),
@@ -72,6 +75,7 @@ class TestHillClimbEstimator(unittest.TestCase):
             self.est_rand._legal_operations(
                 model=self.model2,
                 score=self.score_rand,
+                structure_score=self.score_structure_prior,
                 tabu_list=set(),
                 max_indegree=float("inf"),
                 black_list=set([("A", "B"), ("A", "C"), ("C", "A"), ("C", "B")]),
@@ -92,6 +96,7 @@ class TestHillClimbEstimator(unittest.TestCase):
             self.est_rand._legal_operations(
                 model=self.model2,
                 score=self.score_rand,
+                structure_score=self.score_structure_prior,
                 tabu_list=set(),
                 max_indegree=float("inf"),
                 black_list=set(),
@@ -118,6 +123,7 @@ class TestHillClimbEstimator(unittest.TestCase):
         legal_ops = self.est_titanic1._legal_operations(
             model=start_model,
             score=self.score_titanic1,
+            structure_score=self.score_structure_prior,
             tabu_list=[],
             max_indegree=float("inf"),
             black_list=set(),
@@ -134,6 +140,7 @@ class TestHillClimbEstimator(unittest.TestCase):
         legal_ops_tabu = self.est_titanic1._legal_operations(
             model=start_model,
             score=self.score_titanic1,
+            structure_score=self.score_structure_prior,
             tabu_list=tabu_list,
             max_indegree=float("inf"),
             black_list=set(),
@@ -145,6 +152,7 @@ class TestHillClimbEstimator(unittest.TestCase):
         legal_ops_indegree = self.est_titanic1._legal_operations(
             model=start_model,
             score=self.score_titanic1,
+            structure_score=self.score_structure_prior,
             tabu_list=[],
             max_indegree=1,
             black_list=set(),
@@ -156,6 +164,7 @@ class TestHillClimbEstimator(unittest.TestCase):
         legal_ops_both = self.est_titanic1._legal_operations(
             model=start_model,
             score=self.score_titanic1,
+            structure_score=self.score_structure_prior,
             tabu_list=tabu_list,
             max_indegree=1,
             black_list=set(),
diff --git a/pgmpy/tests/test_estimators/test_MaximumLikelihoodEstimator.py b/pgmpy/tests/test_estimators/test_MaximumLikelihoodEstimator.py
index 7977e7ff..92379213 100644
--- a/pgmpy/tests/test_estimators/test_MaximumLikelihoodEstimator.py
+++ b/pgmpy/tests/test_estimators/test_MaximumLikelihoodEstimator.py
@@ -12,6 +12,7 @@ class TestMLE(unittest.TestCase):
     def setUp(self):
         self.m1 = BayesianModel([("A", "C"), ("B", "C")])
         self.model_latents = BayesianModel([("A", "C"), ("B", "C")], latents=["C"])
+        self.data_latents = pd.DataFrame(data={"A": [0, 0, 1], "B": [0, 1, 0]})
 
         self.d1 = pd.DataFrame(data={"A": [0, 0, 1], "B": [0, 1, 0], "C": [1, 1, 0]})
         self.d2 = pd.DataFrame(
@@ -37,7 +38,10 @@ class TestMLE(unittest.TestCase):
 
     def test_error_latent_model(self):
         self.assertRaises(
-            ValueError, MaximumLikelihoodEstimator, self.model_latents, self.d1
+            ValueError,
+            MaximumLikelihoodEstimator,
+            self.model_latents,
+            self.data_latents,
         )
 
     def test_get_parameters_incomplete_data(self):
diff --git a/pgmpy/tests/test_inference/test_CausalInference.py b/pgmpy/tests/test_inference/test_CausalInference.py
index bb5c60fb..5de74ed0 100644
--- a/pgmpy/tests/test_inference/test_CausalInference.py
+++ b/pgmpy/tests/test_inference/test_CausalInference.py
@@ -26,6 +26,58 @@ class TestCausalGraphMethods(unittest.TestCase):
         self.inference.is_valid_backdoor_adjustment_set("X", "Y", Z="C")
 
 
+class TestAdjustmentSet(unittest.TestCase):
+    def setUp(self):
+        # Model example taken from Constructing Separators and Adjustment Sets
+        # in Ancestral Graphs UAI 2014.
+        self.model = BayesianModel(
+            [("x1", "y1"), ("x1", "z1"), ("z1", "z2"), ("z2", "x2"), ("y2", "z2")]
+        )
+        self.infer = CausalInference(self.model)
+
+    def test_proper_backdoor_graph_error(self):
+        self.assertRaises(
+            ValueError, self.infer.get_proper_backdoor_graph, X=["x3"], Y=["y1", "y2"]
+        )
+        self.assertRaises(
+            ValueError, self.infer.get_proper_backdoor_graph, X=["x2"], Y=["y1", "y3"]
+        )
+        self.assertRaises(
+            ValueError,
+            self.infer.get_proper_backdoor_graph,
+            X=["x3", "x2"],
+            Y=["y1", "y3"],
+        )
+
+    def test_proper_backdoor_graph(self):
+        bd_graph = self.infer.get_proper_backdoor_graph(X=["x1", "x2"], Y=["y1", "y2"])
+        self.assertTrue(("x1", "y1") not in bd_graph.edges())
+        self.assertEqual(len(bd_graph.edges()), 4)
+        self.assertTrue(
+            set(bd_graph.edges()),
+            set([("x1", "z1"), ("z1", "z2"), ("z2", "x2"), ("y2", "z2")]),
+        )
+
+    def test_is_valid_adjustment_set(self):
+        self.assertTrue(
+            self.infer.is_valid_adjustment_set(
+                X=["x1", "x2"], Y=["y1", "y2"], adjustment_set=["z1", "z2"]
+            )
+        )
+
+        self.assertFalse(
+            self.infer.is_valid_adjustment_set(
+                X=["x1", "x2"], Y=["y1", "y2"], adjustment_set=["z1"]
+            )
+        )
+
+        self.assertTrue(
+            self.infer.is_valid_adjustment_set(
+                X=["x1", "x2"], Y=["y1", "y2"], adjustment_set=["z2"]
+            )
+        )
+
+
 class TestBackdoorPaths(unittest.TestCase):
     """
     These tests are drawn from games presented in The Book of Why by Judea Pearl. See the Jupyter Notebook called