Repository: pgmpy/pgmpy
Pull request #707 (2016-06-28 09:46:30)

Canonical factor
This PR deals with the CanonicalFactor class.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6779844/badge)](https://coveralls.io/builds/6779844)

Coverage decreased (-0.5%) to 95.719% when pulling **8ef8f1f552a40556f16fe631c2cc66cfa0e1291b on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6779844/badge)](https://coveralls.io/builds/6779844)

Coverage decreased (-0.5%) to 95.719% when pulling **8ef8f1f552a40556f16fe631c2cc66cfa0e1291b on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6796198/badge)](https://coveralls.io/builds/6796198)

Coverage decreased (-0.8%) to 95.435% when pulling **b07c75f3e74aaf994c191b8371c11be169f20b4a on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6830391/badge)](https://coveralls.io/builds/6830391)

Coverage decreased (-1.03%) to 95.214% when pulling **82c80810fab3e9c20cd7c8bef45476abb7970e1b on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6830417/badge)](https://coveralls.io/builds/6830417)

Coverage decreased (-1.02%) to 95.225% when pulling **82c80810fab3e9c20cd7c8bef45476abb7970e1b on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6869480/badge)](https://coveralls.io/builds/6869480)

Coverage decreased (-1.04%) to 95.203% when pulling **6e3c164d1ff836e559b0c129c0053acec59c2f86 on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6885543/badge)](https://coveralls.io/builds/6885543)

Coverage decreased (-1.03%) to 95.214% when pulling **51a468c5ee4ca4c9197e56890804236bfc793a6e on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6885591/badge)](https://coveralls.io/builds/6885591)

Coverage decreased (-1.04%) to 95.203% when pulling **51a468c5ee4ca4c9197e56890804236bfc793a6e on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/6922136/badge)](https://coveralls.io/builds/6922136)

Coverage decreased (-0.8%) to 95.417% when pulling **9902999792b736be74fe6a9f5ddbd902e363570d on yashu-seth:canonical_factor** into **6401506f356190301691a3a8c3482d12d593fa29 on pgmpy:dev**.

ankurankan: @yashu-seth Yes, changing pdf to an internal attribute makes sense. But changing assignment to pdf won't be consistent with Discrete Factor. So, I would suggest you to leave that. Also you need to finish `to_canonical_factor` in this PR. I am merging this right now. Please create a new PR for these changes.
yashu-seth: @ankurankan But i still need to add more tests in this PR. What about the
JointGaussianDistribution PR?
On 9 Jul 2016 19:34, "Ankur Ankan" <notifications@github.com> wrote:

> @yashu-seth <https://github.com/yashu-seth> Yes, changing pdf to an
> internal attribute makes sense. But changing assignment to pdf won't be
> consistent with Discrete Factor. So, I would suggest you to leave that.
> Also you need to finish to_canonical_factor in this PR. I am merging this
> right now. Please create a new PR for these changes.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/pgmpy/pgmpy/pull/707#issuecomment-231535906>, or mute
> the thread
> <https://github.com/notifications/unsubscribe/ANikoc59pbcdjAYXRl7PodpwzOwvW_GEks5qT6pygaJpZM4I_7Jp>
> .
>

yashu-seth: @ankurankan There are still some things to be done on CanonicalFactor. I found some bugs during adding the tests. I was thinking that the PR #693 would get merged first then I will rebase its changes here.
ankurankan: @yashu-seth Please add the tests in a new PR. All the changes of that PR were in this PR. Right ?
ankurankan: @yashu-seth Okay I will merge that now. 
yashu-seth: @ankurankan Yes but I rebased all those commits into 1 commit thinking I will remove that commit when I rebase the changes from PR #693. The work on CanonicalFactor is incomplete. But you can go ahead with PR #693 .
yashu-seth: @ankurankan Okay.

Patch:
diff --git a/pgmpy/factors/__init__.py b/pgmpy/factors/__init__.py
index 5f580335..9402b76b 100644
--- a/pgmpy/factors/__init__.py
+++ b/pgmpy/factors/__init__.py
@@ -5,6 +5,7 @@ from .JointProbabilityDistribution import JointProbabilityDistribution
 from .continuous.ContinuousNode import ContinuousNode
 from .continuous.ContinuousFactor import ContinuousFactor
 from .continuous.JointGaussianDistribution import JointGaussianDistribution
+from .continuous.CanonicalFactor import CanonicalFactor
 
 __all__ = ['Factor',
            'State',
@@ -17,4 +18,5 @@ __all__ = ['Factor',
            'factorset_divide',
            'ContinuousNode',
            'ContinuousFactor',
-           'JointGaussianDistribution']
+           'JointGaussianDistribution',
+           'CanonicalFactor']
diff --git a/pgmpy/factors/continuous/CanonicalFactor.py b/pgmpy/factors/continuous/CanonicalFactor.py
new file mode 100644
index 00000000..acf367b8
--- /dev/null
+++ b/pgmpy/factors/continuous/CanonicalFactor.py
@@ -0,0 +1,477 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import division
+
+import numpy as np
+from scipy.stats import multivariate_normal
+
+from pgmpy.extern import six
+from pgmpy.factors import ContinuousFactor
+from pgmpy.factors import JointGaussianDistribution
+
+
+class CanonicalFactor(ContinuousFactor):
+    u"""
+    The intermediate factors in a Gaussian network can be described
+    compactly using a simple parametric representation called the
+    canonical form. This representation is closed under the basic
+    operations used in inference: factor product, factor division,
+    factor reduction, and marginalization. Thus, we define this
+    CanonicalFactor class that allows the inference process to be
+    performed on joint Gaussian networks.
+
+    A canonical form C (X; K,h,g) is defined as
+
+    C (X; K,h,g) = exp( ((-1/2) * X.T * K * X) + (h.T * X) + g)
+
+    Reference
+    ---------
+    Probabilistic Graphical Models, Principles and Techniques,
+    Daphne Koller and Nir Friedman, Section 14.2, Chapter 14.
+
+    """
+    def __init__(self, variables, K, h, g):
+        """
+        Parameters
+        ----------
+        variables: list or array-like
+        The variables for wich the distribution is defined.
+
+        K: n x n, 2-d array-like
+
+        h : n x 1, array-like
+
+        g : int, float
+
+        pdf: function
+        The probability density function of the distribution.
+
+        The terms K, h and g are defined parameters for canonical
+        factors representation.
+
+        Examples
+        --------
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi = CanonicalFactor(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
+                                  np.array([[1], [-1]]), -3)
+        >>> phi.variables
+        ['X', 'Y']
+
+        >>> phi.K
+        array([[1, -1],
+               [-1, 1]])
+
+        >>> phi.h
+        array([[1],
+               [-1]])
+
+        >>> phi.g
+        -3
+
+        """
+        no_of_var = len(variables)
+
+        if len(h) != no_of_var:
+            raise ValueError("Length of h parameter vector must be equal to the"
+                             "number of variables.")
+
+        self.h = np.asarray(np.reshape(h, (no_of_var, 1)), dtype=float)
+        self.g = g
+        self.K = np.asarray(K, dtype=float)
+
+        if self.K.shape != (no_of_var, no_of_var):
+            raise ValueError("The K matrix should be a square matrix with order equal to"
+                             "the number of variables. Got: {got_shape}, Expected: {exp_shape}".format
+                             (got_shape=self.K.shape, exp_shape=(no_of_var, no_of_var)))
+
+        super(CanonicalFactor, self).__init__(variables, None)
+
+    @property
+    def pdf(self):
+        def fun(*args):
+            x = np.array(args)
+            return np.exp(self.g + np.dot(x, self.h)[0] - 0.5 * np.dot(x.T, np.dot(self.K, x)))
+        return fun
+
+    def copy(self):
+        """
+        Makes a copy of the factor.
+
+        Returns
+        -------
+        CanonicalFactor object: Copy of the factor
+
+        Examples
+        --------
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi = CanonicalFactor(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
+                                  np.array([[1], [-1]]), -3)
+        >>> phi.variables
+        ['X', 'Y']
+
+        >>> phi.K
+        array([[1, -1],
+               [-1, 1]])
+
+        >>> phi.h
+        array([[1],
+               [-1]])
+
+        >>> phi.g
+        -3
+
+        >>> phi2 = phi.copy()
+
+        >>> phi2.variables
+        ['X', 'Y']
+
+        >>> phi2.K
+        array([[1, -1],
+               [-1, 1]])
+
+        >>> phi2.h
+        array([[1],
+               [-1]])
+
+        >>> phi2.g
+        -3
+
+        """
+        copy_factor = CanonicalFactor(self.scope(), self.K.copy(),
+                                      self.h.copy(), self.g)
+
+        return copy_factor
+
+    def to_joint_gaussian(self):
+        """
+        Return an equivalent Joint Gaussian Distribution.
+
+        Examples
+        --------
+
+        >>> import numpy as np
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi = CanonicalFactor(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
+                                  np.array([[5], [-1]]), 1)
+        >>> jgd = phi.to_joint_gaussian()
+        >>> jgd.variables
+        ['x1', 'x2']
+        >>> jgd.covariance
+        array([[ 0.5  ,  0.25 ],
+               [ 0.25 ,  0.375]])
+        >>> jgd.mean
+        array([[ 2.25 ],
+               [ 0.875]])
+
+        """
+        covariance = np.linalg.inv(self.K)
+        mean = np.dot(covariance, self.h)
+
+        return JointGaussianDistribution(self.scope(), mean, covariance)
+
+    def reduce(self, values, inplace=True):
+        """
+        Reduces the distribution to the context of the given variable values.
+
+        Let C(X,Y ; K, h, g) be some canonical form over X,Y where,
+
+        k = [[K_XX, K_XY],        ;       h = [[h_X],
+             [K_YX, K_YY]]                     [h_Y]]
+
+        The formula for the obtained conditional distribution for setting
+        Y = y is given by,
+
+        .. math:: K' = K_{XX}
+        .. math:: h' = h_X - K_{XY} * y
+        .. math:: g' = g + {h^T}_Y * y - 0.5 * y^T * K_{YY} * y
+
+
+        Parameters
+        ----------
+        values: list, array-like
+            A list of tuples of the form (variable name, variable value).
+
+        inplace: boolean
+            If inplace=True it will modify the factor itself, else would return
+            a new CaninicalFactor object.
+
+        Returns
+        -------
+        CanonicalFactor or None:
+                if inplace=True (default) returns None
+                if inplace=False returns a new CanonicalFactor instance.
+
+        Examples
+        --------
+        >>> imort numpy as np
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi = CanonicalFactor(['X1', 'X2', 'X3'],
+        ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+        ...                       np.array([[1], [4], [-1]]), -2)
+        >>> phi.variables
+        ['X1', 'X2', 'X3']
+
+        >>> phi.K
+        array([[ 1., -1.],
+                [-1.,  3.]])
+
+        >>> phi.h
+        array([[ 1. ],
+                [ 3.5]])
+
+        >>> phi.g
+        -2
+
+        >>> phi.reduce([('X3', 0.25)])
+
+        >>> phi.variables
+        ['X1', 'X2']
+
+        >>> phi.K
+        array([[ 1, -1],
+                [-1,  4]])
+
+        >>> phi.h
+        array([[ 1. ],
+                [ 4.5]])
+
+        >>> phi.g
+        -2.375
+        """
+        if not isinstance(values, (list, tuple, np.ndarray)):
+            raise TypeError("variables: Expected type list or array-like, "
+                            "got type {var_type}".format(var_type=type(values)))
+
+        if not all([var in self.variables for var, value in values]):
+            raise ValueError("Variable not in scope.")
+
+        phi = self if inplace else self.copy()
+
+        var_to_reduce = [var for var, value in values]
+
+        # index_to_keep -> j vector
+        index_to_keep = [self.variables.index(var) for var in self.variables
+                         if var not in var_to_reduce]
+        # index_to_reduce -> i vector
+        index_to_reduce = [self.variables.index(var) for var in var_to_reduce]
+
+        K_i_i = self.K[np.ix_(index_to_keep, index_to_keep)]
+        K_i_j = self.K[np.ix_(index_to_keep, index_to_reduce)]
+        K_j_j = self.K[np.ix_(index_to_reduce, index_to_reduce)]
+        h_i = self.h[index_to_keep]
+        h_j = self.h[index_to_reduce]
+
+        # The values for the reduced variables.
+        y = np.array([value for index,
+                      value in sorted([(self.variables.index(var), value) for var,
+                                      value in values])]).reshape(len(index_to_reduce), 1)
+
+        phi.variables = [self.variables[index] for index in index_to_keep]
+        phi.K = K_i_i
+        phi.h = h_i - np.dot(K_i_j, y)
+        phi.g = self.g + (np.dot(h_j.T, y) - (0.5 * np.dot(np.dot(y.T, K_j_j), y)))[0][0]
+
+        if not inplace:
+            return phi
+
+    def marginalize(self, variables, inplace=True):
+        u"""
+        Modifies the factor with marginalized values.
+
+        Let C(X,Y ; K, h, g) be some canonical form over X,Y where,
+
+        k = [[K_XX, K_XY],        ;       h = [[h_X],
+             [K_YX, K_YY]]                     [h_Y]]
+
+        In this case, the result of the integration operation is a canonical
+        from C (K', h', g') given by,
+
+        .. math:: K' = K_{XX} - K_{XY} * {K^{-1}}_{YY} * K_YX
+        .. math:: h' = h_X - K_{XY} * {K^{-1}}_{YY} * h_Y
+        .. math:: g' = g + 0.5 * (|Y| * log(2*pi) - log(|K_{YY}|) + {h^T}_Y * K_{YY} * h_Y)
+
+        Parameters
+        ----------
+
+        variables: list or array-like
+                List of variables over which to marginalize.
+
+        inplace: boolean
+                If inplace=True it will modify the distribution itself,
+                else would return a new distribution.
+
+        Returns
+        -------
+        CanonicalFactor or None :
+                if inplace=True (default) returns None
+                if inplace=False return a new CanonicalFactor instance
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi = CanonicalFactor(['X1', 'X2', 'X3'],
+        ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+        ...                       np.array([[1], [4], [-1]]), -2)
+        >>> phi.K
+        array([[ 1, -1,  0],
+                [-1,  4, -2],
+                [ 0, -2,  4]])
+
+        >>> phi.h
+        array([[ 1],
+                [ 4],
+                [-1]])
+
+        >>> phi.g
+        -2
+
+        >>> phi.marginalize(['X3'])
+
+        >>> phi.K
+        array([[ 1., -1.],
+                [-1.,  3.]])
+
+        >>> phi.h
+        array([[ 1. ],
+                [ 3.5]])
+
+        >>> phi.g
+        0.22579135
+        """
+        if not isinstance(variables, (list, tuple, np.ndarray)):
+            raise TypeError("variables: Expected type list or array-like, "
+                            "got type {var_type}".format(var_type=type(variables)))
+
+        if not all([var in self.variables for var in variables]):
+            raise ValueError("Variable not in scope.")
+
+        phi = self if inplace else self.copy()
+
+        # index_to_keep -> i vector
+        index_to_keep = [self.variables.index(var) for var in self.variables
+                         if var not in variables]
+        # index_to_marginalize -> j vector
+        index_to_marginalize = [self.variables.index(var) for var in variables]
+
+        K_i_i = self.K[np.ix_(index_to_keep, index_to_keep)]
+        K_i_j = self.K[np.ix_(index_to_keep, index_to_marginalize)]
+        K_j_i = self.K[np.ix_(index_to_marginalize, index_to_keep)]
+        K_j_j = self.K[np.ix_(index_to_marginalize, index_to_marginalize)]
+        K_j_j_inv = np.linalg.inv(K_j_j)
+        h_i = self.h[index_to_keep]
+        h_j = self.h[index_to_marginalize]
+
+        phi.variables = [self.variables[index] for index in index_to_keep]
+
+        phi.K = K_i_i - np.dot(np.dot(K_i_j, K_j_j_inv), K_j_i)
+        phi.h = h_i - np.dot(np.dot(K_i_j, K_j_j_inv), h_j)
+        phi.g = self.g + 0.5 * (len(variables) * np.log(2 * np.pi) -
+                                np.log(abs(np.linalg.det(K_j_j))) + np.dot(np.dot(h_j.T, K_j_j), h_j))[0][0]
+
+        if not inplace:
+            return phi
+
+    def _operate(self, other, operation, inplace=True):
+        """
+        Gives the CanonicalFactor operation (product or divide) with
+        the other factor.
+
+        The product of two canonical factors over the same scope
+        X is simply:
+
+        C(K1, h1, g1) * C(K2, h2, g2) = C(K1+K2, h1+h2, g1+g2)
+
+        The division of canonical forms is defined analogously:
+
+        C(K1, h1, g1) / C(K2, h2, g2) = C(K1-K2, h1-h2, g1- g2)
+
+        When we have two canonical factors over different scopes X and Y,
+        we simply extend the scope of both to make their scopes match and
+        then perform the operation of the above equation. The extension of
+        the scope is performed by simply adding zero entries to both the K
+        matrices and the h vectors.
+
+        Parameters
+        ----------
+        other: CanonicalFactor
+            The CanonicalFactor to be multiplied.
+
+        operation: String
+            'product' for multiplication operation and
+            'divide' for division operation.
+
+        Returns
+        -------
+        CanonicalFactor or None:
+                        if inplace=True (default) returns None
+                        if inplace=False returns a new CanonicalFactor instance.
+
+        Example
+        -------
+        >>> import numpy as np
+        >>> from pgmpy.factors import CanonicalFactor
+        >>> phi1 = CanonicalFactor(['x1', 'x2', 'x3'],
+                                   np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                                   np.array([[1], [4], [-1]]), -2)
+        >>> phi2 = CanonicalFactor(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
+                                   np.array([[5], [-1]]), 1)
+
+        >>> phi3 = phi1 * phi2
+        >>> phi3.K
+        array([[ 4., -3.,  0.],
+               [-3.,  8., -2.],
+               [ 0., -2.,  4.]])
+        >>> phi3.h
+        array([ 6.,  3., -1.])
+        >>> phi3.g
+        -1
+
+        >>> phi4 = phi1 / phi2
+        >>> phi4.K
+        array([[-2.,  1.,  0.],
+               [ 1.,  0., -2.],
+               [ 0., -2.,  4.]])
+        >>> phi4.h
+        array([-4.,  5., -1.])
+        >>> phi4.g
+        -3
+
+        """
+        if not isinstance(other, CanonicalFactor):
+            raise TypeError("CanonicalFactor object can only be multiplied or divided with "
+                            "an another CanonicalFactor object. Got {other_type}, expected "
+                            "CanonicalFactor.".format(other_type=type(other)))
+
+        phi = self if inplace else self.copy()
+
+        all_vars = self.variables + [var for var in other.variables if var not in self.variables]
+        no_of_var = len(all_vars)
+
+        self_var_index = [all_vars.index(var) for var in self.variables]
+        other_var_index = [all_vars.index(var) for var in other.variables]
+
+        def _extend_K_scope(K, index):
+            ext_K = np.zeros([no_of_var, no_of_var])
+            ext_K[np.ix_(index, index)] = K
+            return ext_K
+
+        def _extend_h_scope(h, index):
+            ext_h = np.zeros(no_of_var).reshape(no_of_var, 1)
+            ext_h[index] = h
+            return ext_h
+
+        phi.variables = all_vars
+
+        if operation == 'product':
+            phi.K = _extend_K_scope(self.K, self_var_index) + _extend_K_scope(other.K, other_var_index)
+            phi.h = _extend_h_scope(self.h, self_var_index) + _extend_h_scope(other.h, other_var_index)
+            phi.g = self.g + other.g
+
+        else:
+            phi.K = _extend_K_scope(self.K, self_var_index) - _extend_K_scope(other.K, other_var_index)
+            phi.h = _extend_h_scope(self.h, self_var_index) - _extend_h_scope(other.h, other_var_index)
+            phi.g = self.g - other.g
+
+        if not inplace:
+            return phi
diff --git a/pgmpy/factors/continuous/ContinuousFactor.py b/pgmpy/factors/continuous/ContinuousFactor.py
index 089b8bbe..7cb28cda 100644
--- a/pgmpy/factors/continuous/ContinuousFactor.py
+++ b/pgmpy/factors/continuous/ContinuousFactor.py
@@ -43,7 +43,14 @@ class ContinuousFactor(object):
             raise ValueError("Variable names cannot be same.")
 
         self.variables = list(variables)
-        self.pdf = pdf
+        self._pdf = pdf
+
+    @property
+    def pdf(self):
+        """
+        Returns the pdf of the ContinuousFactor.
+        """
+        return self._pdf
 
     def scope(self):
         """
@@ -209,7 +216,7 @@ class ContinuousFactor(object):
             return pdf(*reduced_args, **reduced_kwargs)
 
         phi.variables = var_to_keep
-        phi.pdf = reduced_pdf
+        phi._pdf = reduced_pdf
 
         if not inplace:
             return phi
@@ -274,7 +281,7 @@ class ContinuousFactor(object):
         def marginalized_pdf(*args):
             return integrate.nquad(reordered_pdf, [[-np.inf, np.inf] for i in range(len(variables))], args=args)[0]
 
-        phi.pdf = marginalized_pdf
+        phi._pdf = marginalized_pdf
         phi.variables = var_to_keep
 
         if not inplace:
@@ -315,7 +322,7 @@ class ContinuousFactor(object):
 
         pdf_mod = integrate.nquad(pdf, [[-np.inf, np.inf] for var in self.variables])[0]
 
-        phi.pdf = lambda *args: pdf(*args) / pdf_mod
+        phi._pdf = lambda *args: pdf(*args) / pdf_mod
 
         if not inplace:
             return phi
@@ -366,7 +373,7 @@ class ContinuousFactor(object):
                 return pdf(*self_pdf_args) / other.pdf(*other_pdf_args)
 
         phi.variables = modified_pdf_var
-        phi.pdf = modified_pdf
+        phi._pdf = modified_pdf
 
         if not inplace:
             return phi
diff --git a/pgmpy/factors/continuous/JointGaussianDistribution.py b/pgmpy/factors/continuous/JointGaussianDistribution.py
index 3959660c..91b6adba 100644
--- a/pgmpy/factors/continuous/JointGaussianDistribution.py
+++ b/pgmpy/factors/continuous/JointGaussianDistribution.py
@@ -55,8 +55,6 @@ class JointGaussianDistribution(ContinuousFactor):
                                  number of variables.")
 
         self.mean = np.asarray(np.reshape(mean, (no_of_var, 1)), dtype=float)
-        normal_pdf = lambda *args: multivariate_normal.pdf(args, self.mean.reshape(1, no_of_var)[0],
-                                                           covariance)
         self.covariance = np.asarray(covariance, dtype=float)
         self._precision_matrix = None
 
@@ -65,7 +63,12 @@ class JointGaussianDistribution(ContinuousFactor):
                               the number of variables. Got: {got_shape}, Expected: {exp_shape}".format
                              (got_shape=self.covariance.shape, exp_shape=(no_of_var, no_of_var)))
 
-        super(JointGaussianDistribution, self).__init__(variables, normal_pdf)
+        super(JointGaussianDistribution, self).__init__(variables, None)
+
+    @property
+    def pdf(self):
+        return lambda *args: multivariate_normal.pdf(args, self.mean.reshape(1, len(self.variables))[0],
+                                                     self.covariance)
 
     @property
     def precision_matrix(self):
@@ -150,7 +153,6 @@ class JointGaussianDistribution(ContinuousFactor):
         phi.mean = phi.mean[index_to_keep]
         phi.covariance = phi.covariance[np.ix_(index_to_keep, index_to_keep)]
         phi._precision_matrix = None
-        phi.pdf = lambda *args: multivariate_normal(args, phi.mean.reshape(1, len(phi.variables)), phi.covariance)
 
         if not inplace:
             return phi
@@ -241,8 +243,6 @@ class JointGaussianDistribution(ContinuousFactor):
         phi.mean = mu_j + np.dot(np.dot(sig_j_i, sig_i_i_inv), x_i - mu_i)
         phi.covariance = sig_j_j - np.dot(np.dot(sig_j_i, sig_i_i_inv), sig_i_j)
         phi._precision_matrix = None
-        phi.pdf = lambda *args: multivariate_normal(args, phi.mean.reshape(1, len(phi.variables)),
-                                                    phi.covariance)
 
         if not inplace:
             return phi
@@ -314,6 +314,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Example
         -------
 
+        >>> import numpy as np
         >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
@@ -331,10 +332,8 @@ class JointGaussianDistribution(ContinuousFactor):
         >>> phi.g
         -6.51533
         """
-        # TODO: This method will return a CanonicalFactor object.
-        # Currently this cannot be used until the CanonicalFactor class is
-        # created. For now the method returns a tuple of the CanonicalFactor
-        # parameters - (K, h, g)
+        from pgmpy.factors import CanonicalFactor
+        
         mu = self.mean
         sigma = self.covariance
 
@@ -343,6 +342,53 @@ class JointGaussianDistribution(ContinuousFactor):
         h = np.dot(K, mu)
 
         g = -(0.5) * np.dot(mu.T, h)[0, 0] - np.log(
-            np.power(2 * np.pi, len(self.variables)/2) * np.power(np.linalg.det(sigma), 0.5))
+            np.power(2 * np.pi, len(self.variables)/2) * np.power(abs(np.linalg.det(sigma)), 0.5))
+
+        return CanonicalFactor(self.scope(), K, h, g)
+
+    def _operate(self, other, operation, inplace=True):
+        """
+        Gives the CanonicalFactor operation (product or divide) with
+        the other factor.
 
-        return K, h, g
+        Parameters
+        ----------
+        other: CanonicalFactor
+            The CanonicalFactor to be multiplied.
+
+        operation: String
+            'product' for multiplication operation and
+            'divide' for division operation.
+
+        Returns
+        -------
+        CanonicalFactor or None:
+                        if inplace=True (default) returns None
+                        if inplace=False returns a new CanonicalFactor instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
+        >>> dis1 = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
+        >>> dis2 = JGD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
+        >>> dis3 = dis1 * dis2
+        >>> dis3.covariance
+        array([[ 3.6,  1. , -0.4, -0.6],
+               [ 1. ,  2.5, -1. , -1.5],
+               [-0.4, -1. ,  1.6,  2.4],
+               [-1. , -2.5,  4. ,  4.5]])
+        >>> dis3.mean
+        array([[ 1.6],
+               [-1.5],
+               [ 1.6],
+               [ 3.5]])
+        """
+        phi = self if inplace else self.copy()
+
+        phi = self.to_canonical_factor()._operate(other.to_canonical_factor(),
+                                                 operation, inplace=False).to_joint_gaussian()
+
+        if not inplace:
+            return phi
diff --git a/pgmpy/inference/continuous/__init__.py b/pgmpy/inference/continuous/__init__.py
index 412db0e8..b5114e61 100644
--- a/pgmpy/inference/continuous/__init__.py
+++ b/pgmpy/inference/continuous/__init__.py
@@ -1,6 +1,6 @@
 from .base import (LeapFrog, ModifiedEuler, BaseSimulateHamiltonianDynamics,
                    BaseGradLogPDF, GradLogPDFGaussian)
-from .sampling import HamiltonianMCda, HamiltonianMC, NoUTurnSampler, NoUTurnSamplerDA
+from .sampling import HamiltonianMCda, HamiltonianMC
 
 
 __all__ = ['LeapFrog',
@@ -9,6 +9,4 @@ __all__ = ['LeapFrog',
            'BaseGradLogPDF',
            'GradLogPDFGaussian',
            'HamiltonianMC',
-           'HamiltonianMCda',
-           'NoUTurnSampler',
-           'NoUTurnSamplerDA']
+           'HamiltonianMCda']
diff --git a/pgmpy/inference/continuous/sampling.py b/pgmpy/inference/continuous/sampling.py
index 788acf6d..9dcbcc37 100644
--- a/pgmpy/inference/continuous/sampling.py
+++ b/pgmpy/inference/continuous/sampling.py
@@ -1,8 +1,6 @@
-# -*- coding: UTF-8 -*-
 """
     A collection of methods for sampling from continuous models in pgmpy
 """
-from __future__ import division
 from math import sqrt
 
 import numpy as np
@@ -51,12 +49,12 @@ class HamiltonianMC(object):
     >>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
     ...                          trajectory_length=2, stepsize=0.4)
     >>> samples
-    rec.array([(1.0, 1.0), (-3.1861687131079086, 3.7940994520145654),
-     (-1.6920542547310844, 6.347410703806017), ...,
-     (-1.8093621120575312, 5.940860883943261),
-     (0.3933248026088032, 6.3853098838119235),
-     (-0.8654072934719572, 6.023803629334816)],
-              dtype=[('x', '<f8'), ('y', '<f8')])
+    array([(5e-324, 5e-324), (-2.2348941964735225, 4.43066330647519),
+           (-2.316454719617516, 7.430291195678112), ...,
+           (-1.1443831048872348, 3.573135519428842),
+           (-0.2325915892988598, 4.155961788010201),
+           (-0.7582492446601238, 3.5416519297297056)], 
+          dtype=[('x', '<f8'), ('y', '<f8')])
 
     >>> samples = np.array([samples[var_name] for var_name in model.variables])
     >>> np.cov(samples)
@@ -263,7 +261,7 @@ class HamiltonianMC(object):
 
         self.acceptance_rate = self.accepted_proposals / num_samples
 
-        if HAS_PANDAS is True:
+        if HAS_PANDAS == True:
             return pd.DataFrame.from_records(samples)
 
         return samples
@@ -311,7 +309,7 @@ class HamiltonianMC(object):
         array([[ 0.1467264 ,  0.27143857],
                [ 4.0371448 ,  0.15871274],
                [ 3.24656208, -1.03742621],
-               ...,
+               ..., 
                [ 6.45975905,  1.97941306],
                [ 4.89007171,  0.15413156],
                [ 5.9528083 ,  1.92983158]])
@@ -396,7 +394,7 @@ class HamiltonianMCda(HamiltonianMC):
     def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog, delta=0.65):
 
         if not isinstance(delta, float) or delta > 1.0 or delta < 0.0:
-            raise ValueError(
+            raise AttributeError(
                 "delta should be a floating value in between 0 and 1")
 
         self.delta = delta
@@ -404,7 +402,7 @@ class HamiltonianMCda(HamiltonianMC):
         super(HamiltonianMCda, self).__init__(model=model, grad_log_pdf=grad_log_pdf,
                                               simulate_dynamics=simulate_dynamics)
 
-    def _adapt_params(self, stepsize, stepsize_bar, h_bar, mu, index_i, alpha, n_alpha=1):
+    def _adapt_params(self, stepsize, stepsize_bar, h_bar, mu, index_i, alpha):
         """
         Run tha adaptation for stepsize for better proposals of position
         """
@@ -414,7 +412,7 @@ class HamiltonianMCda(HamiltonianMC):
         # See equation (6) section 3.2.1 for details
 
         estimate = 1.0 / (index_i + t0)
-        h_bar = (1 - estimate) * h_bar + estimate * (self.delta - alpha / n_alpha)
+        h_bar = (1 - estimate) * h_bar + estimate * (self.delta - alpha)
 
         stepsize = np.exp(mu - sqrt(index_i) / gamma * h_bar)
         i_kappa = index_i ** (-kappa)
@@ -482,7 +480,7 @@ class HamiltonianMCda(HamiltonianMC):
             stepsize = self._find_reasonable_stepsize(initial_pos)
 
         if num_adapt <= 1:  # Return samples genrated using Simple HMC algorithm
-            return HamiltonianMC.sample(self, initial_pos, num_samples, trajectory_length, stepsize)
+            return HamiltonianMC.sample(self, initial_pos, num_samples, stepsize)
 
         # stepsize is epsilon
         mu = np.log(10.0 * stepsize)  # freely chosen point, after each iteration xt(/position) is shrunk towards it
@@ -511,7 +509,7 @@ class HamiltonianMCda(HamiltonianMC):
 
         self.acceptance_rate = self.accepted_proposals / num_samples
 
-        if HAS_PANDAS is True:
+        if HAS_PANDAS == True:
             return pd.DataFrame.from_records(samples)
 
         return samples
@@ -570,7 +568,7 @@ class HamiltonianMCda(HamiltonianMC):
             stepsize = self._find_reasonable_stepsize(initial_pos)
 
         if num_adapt <= 1:  # return sample generated using Simple HMC algorithm
-            for sample in HamiltonianMC.generate_sample(self, initial_pos, num_samples, trajectory_length, stepsize):
+            for sample in HamiltonianMC.generate_sample(self, initial_pos, num_samples, stepsize, trajectory_length):
                 yield sample
             return
         mu = np.log(10.0 * stepsize)
@@ -593,629 +591,3 @@ class HamiltonianMCda(HamiltonianMC):
             yield position_m
 
         self.acceptance_rate = self.accepted_proposals / num_samples
-
-
-class NoUTurnSampler(HamiltonianMCda):
-    """
-    Class for performing sampling in Continuous model
-    using No U Turn Sampler (a variant of Hamiltonian Monte Carlo)
-
-    Parameters:
-    -----------
-    model: An instance pgmpy.models
-        Model from which sampling has to be done
-
-    grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
-        Class to compute the log and gradient log of distribution
-
-    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
-        Class to propose future states of position and momentum in time by simulating
-        HamiltonianDynamics
-
-    Public Methods:
-    ---------------
-    sample()
-    generate_sample()
-
-    Example:
-    --------
-    >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, LeapFrog
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
-    >>> import numpy as np
-    >>> mean = np.array([1, 2, 3])
-    >>> covariance = np.array([[4, 0.1, 0.2], [0.1, 1, 0.3], [0.2, 0.3, 8]])
-    >>> model = JGD(['x', 'y', 'z'], mean, covariance)
-    >>> sampler = NUTS(model=model, grad_log_pdf=None, simulate_dynamics=LeapFrog)
-    >>> samples = sampler.sample(initial_pos=np.array([0.1, 0.9, 0.3]), num_samples=20000,stepsize=0.4)
-    >>> samples
-    rec.array([(0.1, 0.9, 0.3),
-     (-0.27303886844752756, 0.5028580705249155, 0.2895768065049909),
-     (1.7139810571103862, 2.809135711303245, 5.690811523613858), ...,
-     (-0.7742669710786649, 2.092867703984895, 6.139480724333439),
-     (1.3916152816323692, 1.394952482021687, 3.446906546649354),
-     (-0.2726336476939125, 2.6230854954595357, 2.923948403903159)],
-              dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])
-
-    References
-    ----------
-    Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively
-    Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
-    Machine Learning Research 15 (2014) 1351-1381
-    Algorithm 3 : Efficient No-U-Turn Sampler
-    """
-
-    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog):
-
-        super(NoUTurnSampler, self).__init__(model=model, grad_log_pdf=grad_log_pdf,
-                                             simulate_dynamics=simulate_dynamics)
-
-    def _initalize_tree(self, position, momentum, slice_var, stepsize):
-        """
-        Initalizes root node of the tree, i.e depth = 0
-        """
-
-        position_bar, momentum_bar, _ = self.simulate_dynamics(self.model, position, momentum, stepsize,
-                                                               self.grad_log_pdf).get_proposed_values()
-
-        _, logp_bar = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
-
-        hamiltonian = logp_bar - 0.5 * np.dot(momentum_bar, momentum_bar)
-
-        candidate_set_size = slice_var < np.exp(hamiltonian)
-        accept_set_bool = hamiltonian > np.log(slice_var) - 10000  # delta_max = 10000
-
-        return position_bar, momentum_bar, candidate_set_size, accept_set_bool
-
-    def _update_acceptance_criteria(self, position_forward, position_backward, momentum_forward, momentum_backward,
-                                    accept_set_bool, candidate_set_size, candidate_set_size2):
-
-        # criteria1 = I[(θ+ − θ−)·r− ≥ 0]
-        criteria1 = np.dot((position_forward - position_backward), momentum_backward) >= 0
-
-        # criteira2 = I[(θ+ − θ− )·r+ ≥ 0]
-        criteria2 = np.dot((position_forward - position_backward), momentum_forward) >= 0
-
-        accept_set_bool = accept_set_bool and criteria1 and criteria2
-        candidate_set_size += candidate_set_size2
-
-        return accept_set_bool, candidate_set_size
-
-    def _build_tree(self, position, momentum, slice_var, direction, depth, stepsize):
-        """
-        Recursively builds a tree for proposing new position and momentum
-        """
-
-        # Parameter names in algorithm (here -> representation in algorithm)
-        # position -> theta, momentum -> r, slice_var -> u, direction -> v, depth ->j, stepsize -> epsilon
-        # candidate_set_size -> n, accept_set_bool -> s
-        if depth == 0:
-            # Take single leapfrog step in the given direction (direction * stepsize)
-            position_bar, momentum_bar, candidate_set_size, accept_set_bool =\
-                self._initalize_tree(position, momentum, slice_var, direction * stepsize)
-
-            return (position_bar, momentum_bar, position_bar, momentum_bar, position_bar,
-                    candidate_set_size, accept_set_bool)
-
-        else:
-            # Build left and right subtrees
-            (position_backward, momentum_backward, position_forward, momentum_forward, position_bar,
-             candidate_set_size, accept_set_bool) = self._build_tree(position, momentum,
-                                                                     slice_var, direction, depth - 1, stepsize)
-            if accept_set_bool == 1:
-                if direction == -1:
-                    # Build tree in backward direction
-                    (position_backward, momentum_backward, _, _, position_bar2, candidate_set_size2,
-                     accept_set_bool2) = self._build_tree(position_backward, momentum_backward,
-                                                          slice_var, direction, depth - 1, stepsize)
-                else:
-                    # Build tree in forward direction
-                    (_, _, position_forward, momentum_forward, position_bar2, candidate_set_size2,
-                     accept_set_bool2) = self._build_tree(position_forward, momentum_forward,
-                                                          slice_var, direction, depth - 1, stepsize)
-
-                if np.random.rand() < candidate_set_size2 / (candidate_set_size2 + candidate_set_size):
-                    position_bar = position_bar2
-
-                accept_set_bool, candidate_set_size =\
-                    self._update_acceptance_criteria(position_forward, position_backward, momentum_forward,
-                                                     momentum_backward, accept_set_bool2, candidate_set_size,
-                                                     candidate_set_size2)
-
-            return (position_backward, momentum_backward, position_forward, momentum_forward,
-                    position_bar, candidate_set_size, accept_set_bool)
-
-    def _sample(self, position, stepsize):
-        """
-        Returns a sample using a single iteration of NUTS
-        """
-
-        # Re-sampling momentum
-        momentum = np.random.normal(0, 1, len(position))
-
-        # Initializations
-        depth = 0
-        position_backward, position_forward = position, position
-        momentum_backward, momentum_forward = momentum, momentum
-        candidate_set_size = accept_set_bool = 1
-        _, log_pdf = self.model.get_gradient_log_pdf(position, self.grad_log_pdf)
-
-        # Resample slice variable `u`
-        slice_var = np.random.uniform(0, np.exp(log_pdf - 0.5 * np.dot(momentum, momentum)))
-
-        while accept_set_bool == 1:
-            direction = np.random.choice([-1, 1], p=[0.5, 0.5])
-            if direction == -1:
-                # Build a tree in backward direction
-                (position_backward, momentum_backward, _, _, position_bar,
-                 candidate_set_size2, accept_set_bool2) = self._build_tree(position_backward, momentum_backward,
-                                                                           slice_var, direction, depth, stepsize)
-            else:
-                # Build tree in forward direction
-                (_, _, position_forward, momentum_forward, position_bar,
-                 candidate_set_size2, accept_set_bool2) = self._build_tree(position_forward, momentum_forward,
-                                                                           slice_var, direction, depth, stepsize)
-            if accept_set_bool2 == 1:
-                if np.random.rand() < candidate_set_size2 / candidate_set_size:
-                    position = position_bar.copy()
-
-            accept_set_bool, candidate_set_size = self._update_acceptance_criteria(position_forward, position_backward,
-                                                                                   momentum_forward, momentum_backward,
-                                                                                   accept_set_bool2, candidate_set_size,
-                                                                                   candidate_set_size2)
-            depth += 1
-
-        return position
-
-    def sample(self, initial_pos, num_samples, stepsize=None):
-        """
-        Method to return samples using No U Turn Sampler
-
-        Parameters
-        ----------
-        initial_pos: A 1d array like object
-            Vector representing values of parameter position, the starting
-            state in markov chain.
-
-        num_samples: int
-            Number of samples to be generated
-
-        stepsize: float , defaults to None
-            The stepsize for proposing new values of position and momentum in simulate_dynamics
-            If None, then will be choosen suitably
-
-        Returns
-        -------
-        Returns two different types (based on installations)
-
-        pandas.DataFrame: Returns samples as pandas.DataFrame if environment has a installation of pandas
-
-        numpy.recarray: Returns samples in form of numpy recorded arrays (numpy.recarray)
-
-        Examples
-        ---------
-        >>> # If environment has a installation of pandas
-        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
-        >>> import numpy as np
-        >>> mean = np.array([0, 0, 0])
-        >>> covariance = np.array([[6, 0.7, 0.2], [0.7, 3, 0.9], [0.2, 0.9, 1]])
-        >>> model = JGD(['x', 'y', 'z'], mean, covariance)
-        >>> sampler = NUTS(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
-        >>> samples = sampler.sample(initial_pos=np.array([1, 1, 1]), num_samples=10, stepsize=0.4)
-        >>> samples
-                  x         y         z
-        0  1.000000  1.000000  1.000000
-        1  1.760756  0.271543 -0.613309
-        2  1.883387  0.990745 -0.611720
-        3  0.980812  0.340336 -0.916283
-        4  0.781338  0.647220 -0.948640
-        5  0.040308 -1.391406  0.412201
-        6  1.179549 -1.450552  1.105216
-        7  1.100320 -1.313926  1.207815
-        8  1.484520 -1.349247  0.768599
-        9  0.934942 -1.894589  0.471772
-        """
-        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
-        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
-
-        if stepsize is None:
-            stepsize = self._find_reasonable_stepsize(initial_pos)
-
-        types = [(var_name, 'float') for var_name in self.model.variables]
-        samples = np.zeros(num_samples, dtype=types).view(np.recarray)
-
-        samples[0] = tuple(initial_pos)
-        position_m = initial_pos
-
-        for i in range(1, num_samples):
-            # Genrating sample
-            position_m = self._sample(position_m, stepsize)
-            samples[i] = position_m
-
-        if HAS_PANDAS is True:
-            return pd.DataFrame.from_records(samples)
-
-        return samples
-
-    def generate_sample(self, initial_pos, num_samples, stepsize=None):
-        """
-        Returns a generator type object whose each iteration yields a sample
-
-        Parameters
-        ----------
-        initial_pos: A 1d array like object
-            Vector representing values of parameter position, the starting
-            state in markov chain.
-
-        num_samples: int
-            Number of samples to be generated
-
-        stepsize: float , defaults to None
-            The stepsize for proposing new values of position and momentum in simulate_dynamics
-            If None, then will be choosen suitably
-
-        Returns
-        -------
-        genrator: yielding a numpy.array type object for a sample
-
-        Examples
-        ---------
-        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
-        >>> import numpy as np
-        >>> mean = np.array([11, -6])
-        >>> covariance = np.array([[0.7, 0.2], [0.2, 14]])
-        >>> model = JGD(['x', 'y'], mean, covariance)
-        >>> sampler = NUTS(model=model)
-        >>> samples = sampler.generate_sample(initial_pos=np.array([1, 1]), num_samples=10, stepsize=0.4)
-        >>> samples = np.array([sample for sample in samples])
-        >>> samples
-        array([[ 10.26357538,   0.10062725],
-               [ 12.70600336,   0.63392499],
-               [ 10.95523217,  -0.62079273],
-               [ 10.66263031,  -4.08135962],
-               [ 10.59255762,  -8.48085076],
-               [  9.99860242,  -9.47096032],
-               [ 10.5733564 ,  -9.83504745],
-               [ 11.51302059,  -9.49919523],
-               [ 11.31892143,  -8.5873259 ],
-               [ 11.29008667,  -0.43809674]])
-        """
-        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
-        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
-
-        if stepsize is None:
-            stepsize = self._find_reasonable_stepsize(initial_pos)
-
-        position_m = initial_pos
-
-        for _ in range(0, num_samples):
-
-            position_m = self._sample(position_m, stepsize)
-
-            yield position_m
-
-
-class NoUTurnSamplerDA(NoUTurnSampler):
-    """
-    Class for performing sampling in Continuous model
-    using No U Turn sampler with dual averaging for
-    adaptaion of parameter stepsize.
-
-    Parameters:
-    -----------
-    model: An instance pgmpy.models
-        Model from which sampling has to be done
-
-    grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
-        Class to compute the log and gradient log of distribution
-
-    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
-        Class to propose future states of position and momentum in time by simulating
-        HamiltonianDynamics
-
-    delta: float (in between 0 and 1), defaults to 0.65
-        The target HMC acceptance probability
-
-    Public Methods:
-    ---------------
-    sample()
-    generate_sample()
-
-    Example:
-    --------
-    >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
-    >>> import numpy as np
-    >>> mean = np.array([-1, 12, -3])
-    >>> covariance = np.array([[-2, 7, 2], [7, 14, 4], [2, 4, -1]])
-    >>> model = JGD(['x', 'v', 't'], mean, covariance)
-    >>> sampler = NUTSda(model=model)
-    >>> samples = sampler.sample(initial_pos=np.array([0, 0, 0]), num_adapt=10, num_samples=10, stepsize=0.25)
-    >>> samples
-    rec.array([(0.0, 0.0, 0.0),
-     (0.06100992691638076, -0.17118088764170125, 0.14048470935160887),
-     (0.06100992691638076, -0.17118088764170125, 0.14048470935160887),
-     (-0.7451883138013118, 1.7975387358691155, 2.3090698721374436),
-     (-0.6207457594500309, 1.4611049498441024, 2.5890867012835574),
-     (0.24043604780911487, 1.8660976216530618, 3.2508715592645347),
-     (0.21509819341468212, 2.157760225367607, 3.5749582768731476),
-     (0.20699150582681913, 2.0605044285377305, 3.8588980251618135),
-     (0.20699150582681913, 2.0605044285377305, 3.8588980251618135),
-     (0.085332419611991, 1.7556171374575567, 4.49985082288814)],
-              dtype=[('x', '<f8'), ('v', '<f8'), ('t', '<f8')])
-
-    References
-    ----------
-    Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively
-    Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
-    Machine Learning Research 15 (2014) 1351-1381
-    Algorithm 6 : No-U-Turn Sampler with Dual Averaging
-    """
-
-    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog, delta=0.65):
-
-        if not isinstance(delta, float) or delta > 1.0 or delta < 0.0:
-            raise ValueError(
-                "delta should be a floating value in between 0 and 1")
-
-        self.delta = delta
-
-        super(NoUTurnSamplerDA, self).__init__(model=model, grad_log_pdf=grad_log_pdf,
-                                               simulate_dynamics=simulate_dynamics)
-
-    def _build_tree(self, position, momentum, slice_var, direction, depth, stepsize, position0, momentum0):
-        """
-        Recursively builds a tree for proposing new position and momentum
-        """
-        if depth == 0:
-
-            position_bar, momentum_bar, candidate_set_size, accept_set_bool =\
-                self._initalize_tree(position, momentum, slice_var, direction * stepsize)
-
-            alpha = min(1, self._acceptance_prob(position, position_bar, momentum, momentum_bar))
-
-            return (position_bar, momentum_bar, position_bar, momentum_bar, position_bar,
-                    candidate_set_size, accept_set_bool, alpha, 1)
-
-        else:
-            (position_backward, momentum_backward, position_forward, momentum_forward, position_bar,
-             candidate_set_size, accept_set_bool, alpha, n_alpha) =\
-                self._build_tree(position, momentum, slice_var,
-                                 direction, depth - 1, stepsize, position0, momentum0)
-
-            if accept_set_bool == 1:
-                if direction == -1:
-                    # Build tree in backward direction
-                    (position_backward, momentum_backward, _, _, position_bar2, candidate_set_size2, accept_set_bool2,
-                     alpha2, n_alpha2) = self._build_tree(position_backward, momentum_backward, slice_var, direction,
-                                                          depth - 1, stepsize, position0, momentum0)
-                else:
-                    # Build tree in forward direction
-                    (_, _, position_forward, momentum_forward, position_bar2, candidate_set_size2, accept_set_bool2,
-                     alpha2, n_alpha2) = self._build_tree(position_forward, momentum_forward, slice_var, direction,
-                                                          depth - 1, stepsize, position0, momentum0)
-
-                if np.random.rand() < candidate_set_size2 / (candidate_set_size2 + candidate_set_size):
-                    position_bar = position_bar2
-
-                alpha += alpha2
-                n_alpha += n_alpha2
-                accept_set_bool, candidate_set_size =\
-                    self._update_acceptance_criteria(position_forward, position_backward, momentum_forward,
-                                                     momentum_backward, accept_set_bool2, candidate_set_size,
-                                                     candidate_set_size2)
-
-            return (position_backward, momentum_backward, position_forward, momentum_forward, position_bar,
-                    candidate_set_size, accept_set_bool, alpha, n_alpha)
-
-    def _sample(self, position, stepsize):
-        """
-        Returns a sample using a single iteration of NUTS with dual averaging
-        """
-
-        # Re-sampling momentum
-        momentum = np.random.normal(0, 1, len(position))
-
-        # Initializations
-        depth = 0
-        position_backward, position_forward = position, position
-        momentum_backward, momentum_forward = momentum, momentum
-        candidate_set_size = accept_set_bool = 1
-        position_m_1 = position
-        _, log_pdf = self.model.get_gradient_log_pdf(position, self.grad_log_pdf)
-
-        # Resample slice variable `u`
-        slice_var = np.random.uniform(0, np.exp(log_pdf - 0.5 * np.dot(momentum, momentum)))
-
-        while accept_set_bool == 1:
-            direction = np.random.choice([-1, 1], p=[0.5, 0.5])
-            if direction == -1:
-                # Build a tree in backward direction
-                (position_backward, momentum_backward, _, _, position_bar, candidate_set_size2, accept_set_bool2,
-                 alpha, n_alpha) = self._build_tree(position_backward, momentum_backward, slice_var, direction,
-                                                    depth, stepsize, position_m_1, momentum)
-            else:
-                # Build tree in forward direction
-                (_, _, position_forward, momentum_forward, position_bar, candidate_set_size2, accept_set_bool2,
-                 alpha, n_alpha) = self._build_tree(position_forward, momentum_forward, slice_var, direction,
-                                                    depth, stepsize, position_m_1, momentum)
-
-            if accept_set_bool2 == 1:
-                if np.random.rand() < candidate_set_size2 / candidate_set_size:
-                    position = position_bar
-
-            accept_set_bool, candidate_set_size = self._update_acceptance_criteria(position_forward, position_backward,
-                                                                                   momentum_forward, momentum_backward,
-                                                                                   accept_set_bool2, candidate_set_size,
-                                                                                   candidate_set_size2)
-
-            depth += 1
-
-        return position, alpha, n_alpha
-
-    def sample(self, initial_pos, num_adapt, num_samples, stepsize=None):
-        """
-        Returns samples using No U Turn Sampler with dual averaging
-
-        Parameters
-        ----------
-        initial_pos: A 1d array like object
-            Vector representing values of parameter position, the starting
-            state in markov chain.
-
-        num_adapt: int
-            The number of interations to run the adaptation of stepsize
-
-        num_samples: int
-            Number of samples to be generated
-
-        stepsize: float , defaults to None
-            The stepsize for proposing new values of position and momentum in simulate_dynamics
-            If None, then will be choosen suitably
-
-        Returns
-        -------
-        Returns two different types (based on installations)
-
-        pandas.DataFrame: Returns samples as pandas.DataFrame if environment has a installation of pandas
-
-        numpy.recarray: Returns samples in form of numpy recorded arrays (numpy.recarray)
-
-        Examples
-        ---------
-        >>> # If environment has a installation of pandas
-        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
-        >>> import numpy as np
-        >>> mean = np.array([10, -13])
-        >>> covariance = np.array([[16, -3], [-3, 13]])
-        >>> model = JGD(['x', 'y'], mean, covariance)
-        >>> sampler = NUTSda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
-        >>> samples = sampler.sample(initial_pos=np.array([12, -4]), num_adapt=10, num_samples=10, stepsize=0.1)
-        >>> samples
-                   x          y
-        0  12.000000  -4.000000
-        1  11.864821  -3.696109
-        2  10.546986  -4.892169
-        3   8.526596 -21.555793
-        4   8.526596 -21.555793
-        5  11.343194  -6.353789
-        6  -1.583269 -12.802931
-        7  12.411957 -11.704859
-        8  13.253336 -20.169492
-        9  11.295901  -7.665058
-        """
-        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
-        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
-
-        if stepsize is None:
-            stepsize = self._find_reasonable_stepsize(initial_pos)
-
-        if num_adapt <= 1:
-            return NoUTurnSampler(self.model, self.grad_log_pdf,
-                                  self.simulate_dynamics).sample(initial_pos, num_samples, stepsize)
-
-        mu = np.log(10.0 * stepsize)
-        stepsize_bar = 1.0
-        h_bar = 0.0
-
-        types = [(var_name, 'float') for var_name in self.model.variables]
-        samples = np.zeros(num_samples, dtype=types).view(np.recarray)
-        samples[0] = tuple(initial_pos)
-        position_m = initial_pos
-
-        for i in range(1, num_samples):
-
-            position_m, alpha, n_alpha = self._sample(position_m, stepsize)
-            samples[i] = position_m
-
-            if i <= num_adapt:
-                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha, n_alpha)
-            else:
-                stepsize = stepsize_bar
-
-        if HAS_PANDAS is True:
-            return pd.DataFrame.from_records(samples)
-
-        return samples
-
-    def generate_sample(self, initial_pos, num_adapt, num_samples, stepsize=None):
-        """
-        Returns a generator type object whose each iteration yields a sample
-
-        Parameters
-        ----------
-        initial_pos: A 1d array like object
-            Vector representing values of parameter position, the starting
-            state in markov chain.
-
-        num_adapt: int
-            The number of interations to run the adaptation of stepsize
-
-        num_samples: int
-            Number of samples to be generated
-
-        stepsize: float , defaults to None
-            The stepsize for proposing new values of position and momentum in simulate_dynamics
-            If None, then will be choosen suitably
-
-        Returns
-        -------
-        genrator: yielding a numpy.array type object for a sample
-
-        Examples
-        --------
-        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
-        >>> import numpy as np
-        >>> mean = np.array([1, -100])
-        >>> covariance = np.array([[-12, 45], [45, -10]])
-        >>> model = JGD(['a', 'b'], mean, covariance)
-        >>> sampler = NUTSda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
-        >>> samples = sampler.generate_sample(initial_pos=np.array([12, -4]), num_adapt=10,
-        ...                                   num_samples=10, stepsize=0.1)
-        >>> samples
-        <generator object NoUTurnSamplerDA.generate_sample at 0x7f4fed46a4c0>
-        >>> samples_array = np.array([sample for sample in samples])
-        >>> samples_array
-        array([[ 11.89963386,  -4.06572636],
-               [ 10.3453755 ,  -7.5700289 ],
-               [-26.56899659, -15.3920684 ],
-               [-29.97143077, -12.0801625 ],
-               [-29.97143077, -12.0801625 ],
-               [-33.07960829,  -8.90440347],
-               [-55.28263496, -17.31718524],
-               [-55.28263496, -17.31718524],
-               [-56.63440044, -16.03309364],
-               [-63.880094  , -19.19981944]])
-        """
-        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
-        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
-
-        if stepsize is None:
-            stepsize = self._find_reasonable_stepsize(initial_pos)
-
-        if num_adapt <= 1:  # return sample generated using Simple HMC algorithm
-            for sample in NoUTurnSampler(self.model, self.grad_log_pdf,
-                                         self.simulate_dynamics).generate_sample(initial_pos, num_samples, stepsize):
-                yield sample
-            return
-        mu = np.log(10.0 * stepsize)
-
-        stepsize_bar = 1.0
-        h_bar = 0.0
-
-        position_m = initial_pos.copy()
-        num_adapt += 1
-
-        for i in range(1, num_samples + 1):
-
-            position_m, alpha, n_alpha = self._sample(position_m, stepsize)
-
-            if i <= num_adapt:
-                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha, n_alpha)
-            else:
-                stepsize = stepsize_bar
-
-            yield position_m


Test patch:
diff --git a/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
new file mode 100644
index 00000000..92e0d7bc
--- /dev/null
+++ b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
@@ -0,0 +1,187 @@
+import unittest
+
+import numpy as np
+import numpy.testing as np_test
+from scipy.stats import multivariate_normal
+
+from pgmpy.factors import JointGaussianDistribution as JGD
+from pgmpy.factors import CanonicalFactor
+
+
+class TestCanonicalFactor(unittest.TestCase):
+    def test_class_init(self):
+        phi = CanonicalFactor(['x1', ('y', 'z'), 'x3'],
+                              np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                              np.array([[1], [4.7], [-1]]), -2)
+        self.assertEqual(phi.variables, ['x1', ('y', 'z'), 'x3'])
+        np_test.assert_array_equal(phi.K, np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]], dtype=float))
+        np_test.assert_array_equal(phi.h, np.array([[1], [4.7], [-1]], dtype=float))
+        self.assertEqual(phi.g, -2)
+
+        phi = CanonicalFactor(['x1', ('y', 'z'), 'x3'],
+                              np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                              np.array([1, 4.7, -1]), -2)
+        self.assertEqual(phi.variables, ['x1', ('y', 'z'), 'x3'])
+        np_test.assert_array_equal(phi.K, np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]], dtype=float))
+        np_test.assert_array_equal(phi.h, np.array([[1], [4.7], [-1]], dtype=float))
+        self.assertEqual(phi.g, -2)
+
+        phi = CanonicalFactor(['x'], [[1]], [0], 1)
+        self.assertEqual(phi.variables, ['x'])
+        np_test.assert_array_equal(phi.K, np.array([[1]], dtype=float))
+        np_test.assert_array_equal(phi.h, np.array([[0]], dtype=float))
+        self.assertEqual(phi.g, 1)
+
+    def test_class_init_valueerror(self):
+        self.assertRaises(ValueError, CanonicalFactor, ['x1', 'x2', 'x3'],
+                         np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                         np.array([1, 2]), 7)
+        self.assertRaises(ValueError, CanonicalFactor, ['x1', 'x2', 'x3'],
+                         np.array([[1.1, -1, 0], [-1, 4], [0, -2, 4]]),
+                         np.array([1, 2, 3]), 7)
+        self.assertRaises(ValueError, CanonicalFactor, ['x1', 'x2', 'x3'],
+                         np.array([[1.1, -1, 0], [0, -2, 4]]),
+                         np.array([1, 2, 3]), 7)
+        self.assertRaises(ValueError, CanonicalFactor, ['x1', 'x3'],
+                         np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                         np.array([1, 2, 3]), 7)
+
+
+class TestJGDMethods(unittest.TestCase):
+    def setUp(self):
+        self.phi1 = CanonicalFactor(['x1', 'x2', 'x3'],
+                                    np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                                    np.array([[1], [4.7], [-1]]), -2)
+        self.phi2 = CanonicalFactor(['x'], [[1]], [0], 1)
+        self.phi3 = self.phi1.copy()
+
+        self.gauss_phi1 = JGD(['x1', 'x2', 'x3'],
+                              np.array([[3.13043478], [2.44347826], [0.97173913]]),
+                              np.array([[ 1.30434783,  0.43478261,  0.2173913 ],
+                                       [ 0.43478261,  0.47826087,  0.23913043],
+                                       [ 0.2173913 ,  0.23913043,  0.36956522]], dtype=float))
+        self.gauss_phi2 = JGD(['x'], np.array([0]), np.array([[1]]))
+
+    def test_assignment(self):
+        np_test.assert_almost_equal(self.phi1.assignment(1,2,3), 0.0007848640)
+        np_test.assert_almost_equal(self.phi2.assignment(1.2), 1.323129812337)
+
+    def test_to_joint_gaussian(self):
+        jgd1 = self.phi1.to_joint_gaussian()
+        jgd2 = self.phi2.to_joint_gaussian()
+
+        self.assertEqual(jgd1.variables, self.gauss_phi1.variables)
+        np_test.assert_almost_equal(jgd1.covariance, self.gauss_phi1.covariance)
+        np_test.assert_almost_equal(jgd1.mean, self.gauss_phi1.mean)
+
+        self.assertEqual(jgd2.variables, self.gauss_phi2.variables)
+        np_test.assert_almost_equal(jgd2.covariance, self.gauss_phi2.covariance)
+        np_test.assert_almost_equal(jgd2.mean, self.gauss_phi2.mean)
+
+    def test_reduce(self):
+        phi = self.phi1.reduce([('x1', 7)], inplace=False)
+        self.assertEqual(phi.variables, ['x2', 'x3'])
+        np_test.assert_almost_equal(phi.K, np.array([[ 4., -2.], [-2.,  4.]]))
+        np_test.assert_almost_equal(phi.h, np.array([[ 11.7], [ -1. ]]))
+        np_test.assert_almost_equal(phi.g, -21.95)
+
+        phi = self.phi1.reduce([('x1', 4), ('x2', 1.23)], inplace=False)
+        self.assertEqual(phi.variables, ['x3'])
+        np_test.assert_almost_equal(phi.K, np.array([[ 4.]]))
+        np_test.assert_almost_equal(phi.h, np.array([[ 1.46]]))
+        np_test.assert_almost_equal(phi.g, 0.8752)
+
+        self.phi1.reduce([('x1', 7)])
+        self.assertEqual(self.phi1.variables, ['x2', 'x3'])
+        np_test.assert_almost_equal(self.phi1.K, np.array([[ 4., -2.], [-2.,  4.]]))
+        np_test.assert_almost_equal(self.phi1.h, np.array([[ 11.7], [ -1. ]]))
+        np_test.assert_almost_equal(self.phi1.g, -21.95)
+
+        self.phi1 = self.phi3
+        self.phi1.reduce([('x1', 4), ('x2', 1.23)])
+        self.assertEqual(self.phi1.variables, ['x3'])
+        np_test.assert_almost_equal(self.phi1.K, np.array([[ 4.]]))
+        np_test.assert_almost_equal(self.phi1.h, np.array([[ 1.46]]))
+        np_test.assert_almost_equal(self.phi1.g, 0.8752)
+
+        self.phi1 = self.phi3
+
+    def test_marginalize(self):
+        phi = self.phi1.marginalize(['x1'], inplace=False)
+        self.assertEqual(phi.variables, ['x2', 'x3'])
+        np_test.assert_almost_equal(phi.K, np.array([[ 3.090909, -2.], [-2.,  4.]]))
+        np_test.assert_almost_equal(phi.h, np.array([[ 5.6090909], [ -1. ]]))
+        np_test.assert_almost_equal(phi.g, -0.5787165566)
+
+        phi = self.phi1.marginalize(['x1', 'x2'], inplace=False)
+        self.assertEqual(phi.variables, ['x3'])
+        np_test.assert_almost_equal(phi.K, np.array([[ 2.70588235]]))
+        np_test.assert_almost_equal(phi.h, np.array([[ 2.62941176]]))
+        np_test.assert_almost_equal(phi.g, 39.25598935059)
+
+        self.phi1.marginalize(['x1'])
+        self.assertEqual(self.phi1.variables, ['x2', 'x3'])
+        np_test.assert_almost_equal(self.phi1.K, np.array([[ 3.090909, -2.], [-2.,  4.]]))
+        np_test.assert_almost_equal(self.phi1.h, np.array([[ 5.6090909], [ -1. ]]))
+        np_test.assert_almost_equal(self.phi1.g, -0.5787165566)
+
+        self.phi1 = self.phi3
+        self.phi1.marginalize(['x1', 'x2'])
+        self.assertEqual(self.phi1.variables, ['x3'])
+        np_test.assert_almost_equal(self.phi1.K, np.array([[ 2.70588235]]))
+        np_test.assert_almost_equal(self.phi1.h, np.array([[ 2.62941176]]))
+        np_test.assert_almost_equal(self.phi1.g, 39.25598935059)
+
+        self.phi1 = self.phi3
+
+    def test_operate(self):
+        phi1 = self.phi1 * CanonicalFactor(['x2', 'x4'], [[1, 2], [3, 4]], [0, 4.56], -6.78)
+        phi2 = self.phi1 / CanonicalFactor(['x2', 'x3'], [[1, 2], [3, 4]], [0, 4.56], -6.78)
+
+        self.assertEqual(phi1.variables, ['x1', 'x2', 'x3', 'x4'])
+        np_test.assert_almost_equal(phi1.K, np.array([[ 1.1, -1. ,  0. ,  0. ],
+                                                      [-1. ,  5. , -2. ,  2. ],
+                                                      [ 0. , -2. ,  4. ,  0. ],
+                                                      [ 0. ,  3. ,  0. ,  4. ]]))
+        np_test.assert_almost_equal(phi1.h, np.array([[ 1.  ],
+                                                      [ 4.7 ],
+                                                      [-1.  ],
+                                                      [ 4.56]]))
+        np_test.assert_almost_equal(phi1.g, -8.78)
+
+        self.assertEqual(phi2.variables, ['x1', 'x2', 'x3'])
+        np_test.assert_almost_equal(phi2.K, np.array([[ 1.1, -1. ,  0. ],
+                                                      [-1. ,  3. , -4. ],
+                                                      [ 0. , -5. ,  0. ]]))
+        np_test.assert_almost_equal(phi2.h, np.array([[ 1.  ],
+                                                      [ 4.7 ],
+                                                      [-5.56]]))
+        np_test.assert_almost_equal(phi2.g, 4.78)        
+
+    def test_copy(self):
+        copy_phi1 = self.phi1.copy()
+        self.assertEqual(copy_phi1.variables, self.phi1.variables)
+        np_test.assert_array_equal(copy_phi1.K, self.phi1.K)
+        np_test.assert_array_equal(copy_phi1.h, self.phi1.h)
+        np_test.assert_array_equal(copy_phi1.g, self.phi1.g)
+
+        copy_phi1.marginalize(['x1'])
+        self.assertEqual(self.phi1.variables, ['x1', 'x2', 'x3'])
+        np_test.assert_array_equal(self.phi1.K, np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]], dtype=float))
+        np_test.assert_array_equal(self.phi1.h, np.array([[1], [4.7], [-1]], dtype=float))
+        self.assertEqual(self.phi1.g, -2)
+
+        self.phi1.marginalize(['x2'])
+        self.assertEqual(copy_phi1.variables, ['x2', 'x3'])
+        np_test.assert_almost_equal(copy_phi1.K, np.array([[ 3.090909, -2.], [-2.,  4.]]))
+        np_test.assert_almost_equal(copy_phi1.h, np.array([[ 5.6090909], [ -1. ]]))
+        np_test.assert_almost_equal(copy_phi1.g, -0.5787165566)
+
+        self.phi1 = self.phi3
+
+    def tearDown(self):
+        del self.phi1
+        del self.phi2
+        del self.phi3
+        del self.gauss_phi1
+        del self.gauss_phi2
diff --git a/pgmpy/tests/test_factors/test_continuous/test_ContinuousFactor.py b/pgmpy/tests/test_factors/test_continuous/test_ContinuousFactor.py
index ad226b1b..3d0284cc 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_ContinuousFactor.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_ContinuousFactor.py
@@ -13,17 +13,17 @@ class TestContinuousFactor(unittest.TestCase):
         pdf1 = lambda x, y: (np.power(x, 1)*np.power(y, 2))/beta(x, y)
         phi1 = ContinuousFactor(['x', 'y'], pdf1)
         self.assertEqual(phi1.variables, ['x', 'y'])
-        self.assertEqual(phi1.pdf, pdf1)
+        self.assertEqual(phi1._pdf, pdf1)
 
         pdf2 = lambda *args: multivariate_normal.pdf(args, [0, 0], [[1, 0], [0, 1]])
         phi2 = ContinuousFactor(['x1', 'x2'], pdf2)
         self.assertEqual(phi2.variables, ['x1', 'x2'])
-        self.assertEqual(phi2.pdf, pdf2)
+        self.assertEqual(phi2._pdf, pdf2)
 
         pdf3 = lambda x, y, z: z*(np.power(x, 1)*np.power(y, 2))/beta(x, y)
         phi3 = ContinuousFactor(['x', 'y', 'z'], pdf3)
         self.assertEqual(phi3.variables, ['x', 'y', 'z'])
-        self.assertEqual(phi3.pdf, pdf3)
+        self.assertEqual(phi3._pdf, pdf3)
 
     def test_class_init_typeerror(self):
         pdf1 = lambda x, y: (np.power(x, 1)*np.power(y, 2))/beta(x, y)
@@ -366,12 +366,12 @@ class TestContinuousFactorMethods(unittest.TestCase):
 
         copy1.variables = ['A', 'B']
         self.assertEqual(copy4.variables, self.phi1.variables)
-        copy1.pdf = lambda a, b: (a + b) / (a * a + b * b)
+        copy1._pdf = lambda a, b: (a + b) / (a * a + b * b)
         copy1_pdf = lambda a, b: (a + b) / (a * a + b * b)
         self.assertEqual(copy4.pdf, self.phi1.pdf)
         copy4.variables = ['X', 'Y']
         self.assertEqual(copy1.variables, ['A', 'B'])
-        copy4.pdf = lambda a, b: a + b
+        copy4._pdf = lambda a, b: a + b
         for inp in np.random.rand(4, 2):
             self.assertEqual(copy1.pdf(inp[0], inp[1]), copy1_pdf(inp[0], inp[1]))
 
diff --git a/pgmpy/tests/test_inference/test_continuous_sampling.py b/pgmpy/tests/test_inference/test_continuous_sampling.py
index 2d900f46..8951f5a5 100644
--- a/pgmpy/tests/test_inference/test_continuous_sampling.py
+++ b/pgmpy/tests/test_inference/test_continuous_sampling.py
@@ -3,8 +3,7 @@ import unittest
 import numpy as np
 
 from pgmpy.inference.continuous import (HamiltonianMC as HMC, HamiltonianMCda as HMCda, ModifiedEuler,
-                                        GradLogPDFGaussian as GradGaussian, NoUTurnSampler as NUTS,
-                                        NoUTurnSamplerDA as NUTSda)
+                                        GradLogPDFGaussian as GradGaussian)
 from pgmpy.models import JointGaussianDistribution as JGD
 
 
@@ -21,14 +20,12 @@ class TestHMCInference(unittest.TestCase):
             HMCda(model=self.test_model, grad_log_pdf=1)
         with self.assertRaises(TypeError):
             HMCda(model=self.test_model, grad_log_pdf=GradGaussian, simulate_dynamics=1)
-        with self.assertRaises(ValueError):
-            HMCda(model=self.test_model, delta=-1)
         with self.assertRaises(TypeError):
             self.hmc_sampler.sample(initial_pos=1, num_adapt=1, num_samples=1, trajectory_length=1)
         with self.assertRaises(TypeError):
             self.hmc_sampler.generate_sample(1, 1, 1, 1).send(None)
         with self.assertRaises(TypeError):
-            HMC(model=self.test_model).sample(initial_pos=1, num_samples=1, trajectory_length=1)
+            HMC(model=self.test_model).sample(position0=1, num_samples=1, trajectory_length=1)
         with self.assertRaises(TypeError):
             HMC(model=self.test_model).generate_sample(1, 1, 1).send(None)
         # TODO: Remove them after implementation of JointGaussianDistribution
@@ -92,65 +89,3 @@ class TestHMCInference(unittest.TestCase):
     def tearDown(self):
         del self.hmc_sampler
         del self.test_model
-
-
-class TestNUTSInference(unittest.TestCase):
-
-    def setUp(self):
-        mean = np.array([-1, 1, 0])
-        covariance = np.array([[6, 0.7, 0.2], [0.7, 3, 0.9], [0.2, 0.9, 1]])
-        self.test_model = JGD(['x', 'y', 'z'], mean, covariance)
-        self.nuts_sampler = NUTSda(model=self.test_model, grad_log_pdf=GradGaussian)
-
-    def test_errors(self):
-        with self.assertRaises(TypeError):
-            NUTS(model=self.test_model, grad_log_pdf=JGD)
-        with self.assertRaises(TypeError):
-            NUTS(model=self.test_model, grad_log_pdf=None, simulate_dynamics=GradGaussian)
-        with self.assertRaises(ValueError):
-            NUTSda(model=self.test_model, delta=-0.2)
-        with self.assertRaises(ValueError):
-            NUTSda(model=self.test_model, delta=1.1)
-        with self.assertRaises(TypeError):
-            NUTS(self.test_model).sample(initial_pos={1, 1, 1}, num_samples=1)
-        with self.assertRaises(ValueError):
-            NUTS(self.test_model).sample(initial_pos=[1, 1], num_samples=1)
-        with self.assertRaises(TypeError):
-            NUTSda(self.test_model).sample(initial_pos=1, num_samples=1, num_adapt=1)
-        with self.assertRaises(ValueError):
-            NUTSda(self.test_model).sample(initial_pos=[1, 1, 1, 1], num_samples=1, num_adapt=1)
-        with self.assertRaises(TypeError):
-            NUTS(self.test_model).generate_sample(initial_pos=0.1, num_samples=1).send(None)
-        with self.assertRaises(ValueError):
-            NUTS(self.test_model).generate_sample(initial_pos=(0, 1, 1, 1), num_samples=1).send(None)
-        with self.assertRaises(TypeError):
-            NUTSda(self.test_model).generate_sample(initial_pos=[[1, 2, 3]], num_samples=1, num_adapt=1).send(None)
-        with self.assertRaises(ValueError):
-            NUTSda(self.test_model).generate_sample(initial_pos=[1], num_samples=1, num_adapt=1).send(None)
-
-    def test_sampling(self):
-        np.random.seed(1010101)
-        samples = self.nuts_sampler.sample(initial_pos=[-0.4, 1, 3.6], num_adapt=0, num_samples=10000)
-        sample_covariance = np.cov(samples.values.T)
-        self.assertTrue(np.linalg.norm(sample_covariance - self.test_model.covariance) < 3)
-
-        np.random.seed(1210161)
-        samples = self.nuts_sampler.generate_sample(initial_pos=[-0.4, 1, 3.6], num_adapt=0, num_samples=10000)
-        samples_array = np.array([sample for sample in samples])
-        sample_covariance = np.cov(samples_array.T)
-        self.assertTrue(np.linalg.norm(sample_covariance - self.test_model.covariance) < 3)
-
-        np.random.seed(12313131)
-        samples = self.nuts_sampler.sample(initial_pos=[0.2, 0.4, 2.2], num_adapt=10000, num_samples=10000)
-        sample_covariance = np.cov(samples.values.T)
-        self.assertTrue(np.linalg.norm(sample_covariance - self.test_model.covariance) < 0.4)
-
-        np.random.seed(921312312)
-        samples = self.nuts_sampler.generate_sample(initial_pos=[0.2, 0.4, 2.2], num_adapt=10000, num_samples=10000)
-        samples_array = np.array([sample for sample in samples])
-        sample_covariance = np.cov(samples_array.T)
-        self.assertTrue(np.linalg.norm(sample_covariance - self.test_model.covariance) < 0.4)
-
-    def tearDown(self):
-        del self.test_model
-        del self.nuts_sampler