Repository: pgmpy/pgmpy
Pull request #866 (2017-04-01 11:13:00)

Refactor JointGaussianDistribution


Patch:
diff --git a/docs/factors.rst b/docs/factors.rst
index cfdf9c20..32d0acd1 100644
--- a/docs/factors.rst
+++ b/docs/factors.rst
@@ -29,7 +29,7 @@ Continuous
 Canonical Factor
 ^^^^^^^^^^^^^^^^
 
-.. automodule:: pgmpy.factors.continuous.CanonicalFactor
+.. automodule:: pgmpy.factors.continuous.CanonicalDistribution
    :members:
 
 Continuous Factor
@@ -42,7 +42,7 @@ Continuous Factor
 Joint Gaussian Distribution
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. automodule:: pgmpy.factors.continuous.JointGaussianDistribution
+.. automodule:: pgmpy.factors.continuous.GaussianDistribution
    :members:
 
 Linear Gaussian CPD
diff --git a/pgmpy/factors/continuous/CanonicalFactor.py b/pgmpy/factors/continuous/CanonicalDistribution.py
similarity index 66%
rename from pgmpy/factors/continuous/CanonicalFactor.py
rename to pgmpy/factors/continuous/CanonicalDistribution.py
index a40296c9..8c26b163 100644
--- a/pgmpy/factors/continuous/CanonicalFactor.py
+++ b/pgmpy/factors/continuous/CanonicalDistribution.py
@@ -4,18 +4,18 @@ from __future__ import division
 
 import numpy as np
 
-from pgmpy.factors.continuous import ContinuousFactor
-from pgmpy.factors.distributions import JointGaussianDistribution
+from pgmpy.factors.distributions import BaseDistribution
+from pgmpy.factors.distributions import GaussianDistribution
 
 
-class CanonicalFactor(ContinuousFactor):
+class CanonicalDistribution(BaseDistribution):
     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
+    CanonicalDistribution class that allows the inference process to be
     performed on joint Gaussian networks.
 
     A canonical form C (X; K,h,g) is defined as
@@ -49,8 +49,8 @@ class CanonicalFactor(ContinuousFactor):
 
         Examples
         --------
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi = CanonicalFactor(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi = CanonicalDistribution(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
                                   np.array([[1], [-1]]), -3)
         >>> phi.variables
         ['X', 'Y']
@@ -70,19 +70,20 @@ class CanonicalFactor(ContinuousFactor):
         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.")
+            raise ValueError("Length of h parameter vector must be equal to "
+                             "the number of variables.")
 
+        self.variables = 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)
+            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)))
 
     @property
     def pdf(self):
@@ -91,18 +92,39 @@ class CanonicalFactor(ContinuousFactor):
             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 assignment(self, *x):
+        """
+        Returns the probability value of the PDF at the given parameter values.
+        Parameters
+        ----------
+        *x: values of all variables of this distribution,
+            collective defining a point at which the probability value is to be computed.
+        Returns
+        -------
+        float: The probability value at the point.
+        Examples
+        --------
+        >>> from pgmpy.factors.distributions import GaussianDistribution
+        >>> dist = GaussianDistribution(variables=['x1', 'x2'],
+        ...                             mean=[[0], [0]],
+        ...                             covariance=[[1, 0], [0, 1]])
+        >>> dist.assignment(0, 0)
+        0.15915494309189535
+        """
+        return self.pdf(*x)
+
     def copy(self):
         """
         Makes a copy of the factor.
 
         Returns
         -------
-        CanonicalFactor object: Copy of the factor
+        CanonicalDistribution object: Copy of the factor
 
         Examples
         --------
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi = CanonicalFactor(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi = CanonicalDistribution(['X', 'Y'], np.array([[1, -1], [-1, 1]]),
                                   np.array([[1], [-1]]), -3)
         >>> phi.variables
         ['X', 'Y']
@@ -135,8 +157,8 @@ class CanonicalFactor(ContinuousFactor):
         -3
 
         """
-        copy_factor = CanonicalFactor(self.scope(), self.K.copy(),
-                                      self.h.copy(), self.g)
+        copy_factor = CanonicalDistribution(self.variables, self.K.copy(),
+                                            self.h.copy(), self.g)
 
         return copy_factor
 
@@ -148,8 +170,8 @@ class CanonicalFactor(ContinuousFactor):
         --------
 
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi = CanonicalFactor(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi = CanonicalDistribution(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
                                   np.array([[5], [-1]]), 1)
         >>> jgd = phi.to_joint_gaussian()
         >>> jgd.variables
@@ -165,7 +187,7 @@ class CanonicalFactor(ContinuousFactor):
         covariance = np.linalg.inv(self.K)
         mean = np.dot(covariance, self.h)
 
-        return JointGaussianDistribution(self.scope(), mean, covariance)
+        return GaussianDistribution(self.variables, mean, covariance)
 
     def reduce(self, values, inplace=True):
         """
@@ -195,15 +217,15 @@ class CanonicalFactor(ContinuousFactor):
 
         Returns
         -------
-        CanonicalFactor or None:
+        CanonicalDistribution or None:
                 if inplace=True (default) returns None
-                if inplace=False returns a new CanonicalFactor instance.
+                if inplace=False returns a new CanonicalDistribution instance.
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi = CanonicalFactor(['X1', 'X2', 'X3'],
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi = CanonicalDistribution(['X1', 'X2', 'X3'],
         ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
         ...                       np.array([[1], [4], [-1]]), -2)
         >>> phi.variables
@@ -298,15 +320,15 @@ class CanonicalFactor(ContinuousFactor):
 
         Returns
         -------
-        CanonicalFactor or None :
+        CanonicalDistribution or None :
                 if inplace=True (default) returns None
-                if inplace=False return a new CanonicalFactor instance
+                if inplace=False return a new CanonicalDistribution instance
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi = CanonicalFactor(['X1', 'X2', 'X3'],
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi = CanonicalDistribution(['X1', 'X2', 'X3'],
         ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
         ...                       np.array([[1], [4], [-1]]), -2)
         >>> phi.K
@@ -370,7 +392,7 @@ class CanonicalFactor(ContinuousFactor):
 
     def _operate(self, other, operation, inplace=True):
         """
-        Gives the CanonicalFactor operation (product or divide) with
+        Gives the CanonicalDistribution operation (product or divide) with
         the other factor.
 
         The product of two canonical factors over the same scope
@@ -391,7 +413,7 @@ class CanonicalFactor(ContinuousFactor):
         Parameters
         ----------
         other: CanonicalFactor
-            The CanonicalFactor to be multiplied.
+            The CanonicalDistribution to be multiplied.
 
         operation: String
             'product' for multiplication operation and
@@ -399,18 +421,18 @@ class CanonicalFactor(ContinuousFactor):
 
         Returns
         -------
-        CanonicalFactor or None:
+        CanonicalDistribution or None:
                         if inplace=True (default) returns None
-                        if inplace=False returns a new CanonicalFactor instance.
+                        if inplace=False returns a new CanonicalDistribution instance.
 
         Example
         -------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import CanonicalFactor
-        >>> phi1 = CanonicalFactor(['x1', 'x2', 'x3'],
+        >>> from pgmpy.factors.continuous import CanonicalDistribution
+        >>> phi1 = CanonicalDistribution(['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]]),
+        >>> phi2 = CanonicalDistribution(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
                                    np.array([[5], [-1]]), 1)
 
         >>> phi3 = phi1 * phi2
@@ -434,10 +456,11 @@ class CanonicalFactor(ContinuousFactor):
         -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)))
+        if not isinstance(other, CanonicalDistribution):
+            raise TypeError(
+                "CanonicalDistribution object can only be multiplied or divided "
+                "with an another CanonicalDistribution object. Got {other_type}, "
+                "expected CanonicalDistribution.".format(other_type=type(other)))
 
         phi = self if inplace else self.copy()
 
@@ -471,3 +494,95 @@ class CanonicalFactor(ContinuousFactor):
 
         if not inplace:
             return phi
+
+    def product(self, other, inplace=True):
+        """
+        Returns the product of two gaussian distributions.
+
+        Parameters
+        ----------
+        other: CanonicalFactor
+            The GaussianDistribution to be multiplied.
+
+        inplace: boolean
+            If True, modifies the distribution itself, otherwise returns a new
+            CanonicalDistribution object.
+
+        Returns
+        -------
+        CanonicalDistribution or None:
+                    if inplace=True (default) returns None.
+                    if inplace=False returns a new CanonicalDistribution instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis1 = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        ...            np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
+        >>> dis2 = GD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
+        >>> dis3 = dis1.product(dis2, inplace=False)
+        >>> 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]])
+        """
+        return self._operate(other, operation='product', inplace=inplace)
+
+    def divide(self, other, inplace=True):
+        """
+        Returns the division of two gaussian distributions.
+
+        Parameters
+        ----------
+        other: GaussianDistribution
+            The GaussianDistribution to be divided.
+
+        inplace: boolean
+            If True, modifies the distribution itself, otherwise returns a new
+            GaussianDistribution object.
+
+        Returns
+        -------
+        CanonicalDistribution or None:
+                    if inplace=True (default) returns None.
+                    if inplace=False returns a new CanonicalDistribution instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis1 = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        ...            np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
+        >>> dis2 = GD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
+        >>> dis3 = dis1.divide(dis2, inplace=False)
+        >>> 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]])
+        """
+        return self._operate(other, operation='divide', inplace=inplace)
+
+    def __mul__(self, other):
+        return self.product(other, inplace=False)
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+    def __truediv__(self, other):
+        return self.divide(other, inplace=False)
+
+    __div__ = __truediv__
+
diff --git a/pgmpy/factors/continuous/__init__.py b/pgmpy/factors/continuous/__init__.py
index bf76fdc3..4528e4f1 100644
--- a/pgmpy/factors/continuous/__init__.py
+++ b/pgmpy/factors/continuous/__init__.py
@@ -1,10 +1,10 @@
 from .discretize import BaseDiscretizer, RoundingDiscretizer, UnbiasedDiscretizer
 from .ContinuousFactor import ContinuousFactor
-from .CanonicalFactor import CanonicalFactor
+from .CanonicalDistribution import CanonicalDistribution
 from .LinearGaussianCPD import LinearGaussianCPD
 
 
-__all__ = ['CanonicalFactor',
+__all__ = ['CanonicalDistribution',
            'ContinuousFactor',
            'LinearGaussianCPD'
            'BaseDiscretizer',
diff --git a/pgmpy/factors/distributions/JointGaussianDistribution.py b/pgmpy/factors/distributions/GaussianDistribution.py
similarity index 66%
rename from pgmpy/factors/distributions/JointGaussianDistribution.py
rename to pgmpy/factors/distributions/GaussianDistribution.py
index 6a14f89c..ee3ddc5d 100644
--- a/pgmpy/factors/distributions/JointGaussianDistribution.py
+++ b/pgmpy/factors/distributions/GaussianDistribution.py
@@ -5,10 +5,10 @@ from __future__ import division
 import numpy as np
 from scipy.stats import multivariate_normal
 
-from pgmpy.factors.distributions import CustomDistribution
+from pgmpy.factors.distributions import BaseDistribution
 
 
-class JointGaussianDistribution(CustomDistribution):
+class GaussianDistribution(BaseDistribution):
     u"""
     In its most common representation, a multivariate Gaussian distribution
     over X1...........Xn is characterized by an n-dimensional mean vector μ,
@@ -33,8 +33,8 @@ class JointGaussianDistribution(CustomDistribution):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> dis.variables
         ['x1', 'x2', 'x3']
@@ -55,6 +55,7 @@ class JointGaussianDistribution(CustomDistribution):
             raise ValueError("Length of mean_vector must be equal to the\
                                  number of variables.")
 
+        self.variables = variables
         self.mean = np.asarray(np.reshape(mean, (no_of_var, 1)), dtype=float)
         self.covariance = np.asarray(covariance, dtype=float)
         self._precision_matrix = None
@@ -64,13 +65,34 @@ class JointGaussianDistribution(CustomDistribution):
                               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, None)
+#        super(GaussianDistribution, 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)
 
+    def assignment(self, *x):
+        """
+        Returns the probability value of the PDF at the given parameter values.
+        Parameters
+        ----------
+        *x: values of all variables of this distribution,
+            collective defining a point at which the probability value is to be computed.
+        Returns
+        -------
+        float: The probability value at the point.
+        Examples
+        --------
+        >>> from pgmpy.factors.distributions import GaussianDistribution
+        >>> dist = GaussianDistribution(variables=['x1', 'x2'],
+        ...                             mean=[[0], [0]],
+        ...                             covariance=[[1, 0], [0, 1]])
+        >>> dist.assignment(0, 0)
+        0.15915494309189535
+        """
+        return self.pdf(*x)
+
     @property
     def precision_matrix(self):
         """
@@ -79,8 +101,8 @@ class JointGaussianDistribution(CustomDistribution):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> dis.precision_matrix
         array([[ 0.3125    , -0.125     ,  0.        ],
@@ -108,15 +130,15 @@ class JointGaussianDistribution(CustomDistribution):
 
         Returns
         -------
-        JointGaussianDistribution or None :
+        GaussianDistribution or None :
                 if inplace=True (default) returns None
-                if inplace=False return a new JointGaussianDistribution instance
+                if inplace=False return a new GaussianDistribution instance
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> dis.variables
         ['x1', 'x2', 'x3']
@@ -182,15 +204,15 @@ class JointGaussianDistribution(CustomDistribution):
 
         Returns
         -------
-        JointGaussianDistribution or None:
+        GaussianDistribution or None:
                 if inplace=True (default) returns None
-                if inplace=False returns a new JointGaussianDistribution instance.
+                if inplace=False returns a new GaussianDistribution instance.
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> dis.variables
         ['x1', 'x2', 'x3']
@@ -261,13 +283,13 @@ class JointGaussianDistribution(CustomDistribution):
 
         Returns
         -------
-        JointGaussianDistribution: copy of the distribution
+        GaussianDistribution: copy of the distribution
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> gauss_dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> gauss_dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...                 np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> copy_dis = gauss_dis.copy()
         >>> copy_dis.variables
@@ -285,8 +307,8 @@ class JointGaussianDistribution(CustomDistribution):
                 [-0.125     ,  0.58333333,  0.33333333],
                 [ 0.        ,  0.33333333,  0.33333333]])
         """
-        copy_distribution = JointGaussianDistribution(self.variables, self.mean.copy(),
-                                                      self.covariance.copy())
+        copy_distribution = GaussianDistribution(self.variables, self.mean.copy(),
+                                                 self.covariance.copy())
         if self._precision_matrix is not None:
             copy_distribution._precision_matrix = self._precision_matrix.copy()
 
@@ -294,7 +316,7 @@ class JointGaussianDistribution(CustomDistribution):
 
     def to_canonical_factor(self):
         u"""
-        Returns an equivalent CanonicalFactor object.
+        Returns an equivalent CanonicalDistribution object.
 
         The formulas for calculating the cannonical factor parameters
         for N(μ; Σ) = C(K; h; g) are as follows -
@@ -315,8 +337,8 @@ class JointGaussianDistribution(CustomDistribution):
         -------
 
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...             np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> phi = dis.to_canonical_factor()
         >>> phi.variables
@@ -332,7 +354,7 @@ class JointGaussianDistribution(CustomDistribution):
         >>> phi.g
         -6.51533
         """
-        from pgmpy.factors.continuous import CanonicalFactor
+        from pgmpy.factors.continuous import CanonicalDistribution
 
         mu = self.mean
         sigma = self.covariance
@@ -344,17 +366,17 @@ class JointGaussianDistribution(CustomDistribution):
         g = -(0.5) * np.dot(mu.T, h)[0, 0] - np.log(
             np.power(2 * np.pi, len(self.variables)/2) * np.power(abs(np.linalg.det(sigma)), 0.5))
 
-        return CanonicalFactor(self.variables, K, h, g)
+        return CanonicalDistribution(self.variables, K, h, g)
 
     def _operate(self, other, operation, inplace=True):
         """
-        Gives the CanonicalFactor operation (product or divide) with
+        Gives the CanonicalDistribution operation (product or divide) with
         the other factor.
 
         Parameters
         ----------
-        other: CanonicalFactor
-            The CanonicalFactor to be multiplied.
+        other: CanonicalDistribution
+            The CanonicalDistribution to be multiplied.
 
         operation: String
             'product' for multiplication operation and
@@ -362,17 +384,17 @@ class JointGaussianDistribution(CustomDistribution):
 
         Returns
         -------
-        CanonicalFactor or None:
+        CanonicalDistribution or None:
                         if inplace=True (default) returns None
-                        if inplace=False returns a new CanonicalFactor instance.
+                        if inplace=False returns a new CanonicalDistribution instance.
 
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-        >>> dis1 = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis1 = GD(['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]])
+        >>> dis2 = GD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
         >>> dis3 = dis1 * dis2
         >>> dis3.covariance
         array([[ 3.6,  1. , -0.4, -0.6],
@@ -391,3 +413,82 @@ class JointGaussianDistribution(CustomDistribution):
         if not inplace:
             return phi
 
+    def product(self, other, inplace=True):
+        """
+        Returns the product of two gaussian distributions.
+
+        Parameters
+        ----------
+        other: GaussianDistribution
+            The GaussianDistribution to be multiplied.
+
+        inplace: boolean
+            If True, modifies the distribution itself, otherwise returns a new
+            GaussianDistribution object.
+
+        Returns
+        -------
+        CanonicalDistribution or None:
+                    if inplace=True (default) returns None.
+                    if inplace=False returns a new CanonicalDistribution instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis1 = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        ...            np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
+        >>> dis2 = GD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
+        >>> dis3 = dis1.product(dis2, inplace=False)
+        >>> 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]])
+        """
+        return self._operate(other, operation='product', inplace=inplace)
+
+    def divide(self, other, inplace=True):
+        """
+        Returns the division of two gaussian distributions.
+
+        Parameters
+        ----------
+        other: GaussianDistribution
+            The GaussianDistribution to be divided.
+
+        inplace: boolean
+            If True, modifies the distribution itself, otherwise returns a new
+            GaussianDistribution object.
+
+        Returns
+        -------
+        CanonicalDistribution or None:
+                    if inplace=True (default) returns None.
+                    if inplace=False returns a new CanonicalDistribution instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from pgmpy.factors.distributions import GaussianDistribution as GD
+        >>> dis1 = GD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
+        ...            np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
+        >>> dis2 = GD(['x3', 'x4'], [1, 2], [[2, 3], [5, 6]])
+        >>> dis3 = dis1.divide(dis2, inplace=False)
+        >>> 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]])
+        """
+        return self._operate(other, operation='divide', inplace=inplace)
diff --git a/pgmpy/factors/distributions/__init__.py b/pgmpy/factors/distributions/__init__.py
index ce6d8c90..b80b0435 100644
--- a/pgmpy/factors/distributions/__init__.py
+++ b/pgmpy/factors/distributions/__init__.py
@@ -1,9 +1,9 @@
 from .base import BaseDistribution
 from .CustomDistribution import CustomDistribution
-from .JointGaussianDistribution import JointGaussianDistribution
+from .GaussianDistribution import GaussianDistribution
 
 
 __all__ = ['BaseDistribution',
            'CustomDistribution',
-           'JointGaussianDistribution'
+           'GaussianDistribution'
            ]
diff --git a/pgmpy/factors/distributions/base.py b/pgmpy/factors/distributions/base.py
index f85ba5d1..7d632e67 100644
--- a/pgmpy/factors/distributions/base.py
+++ b/pgmpy/factors/distributions/base.py
@@ -2,42 +2,42 @@ from abc import ABCMeta, abstractmethod, abstractproperty
 
 
 class BaseDistribution(object):
-    @abstractproperty
-    def pdf(self):
-        pass
-    
-    @abstractproperty
-    def variables(self):
-        pass
-
-    @abstractmethod
-    def assignment(self, *args, **kwargs):
-        pass
-
-    @abstractmethod
-    def copy(self):
-        pass
-    @abstractmethod
-    def discretize(self, method, *args, **kwargs):
-        pass
-
-    @abstractmethod
-    def reduce(self, values, inplace=True):
-        pass
-
-    @abstractmethod
-    def marginalize(self, variables, inplace=True):
-        pass
-
-    @abstractmethod
-    def normalize(self, inplace=True):
-         pass
-
-    @abstractmethod
-    def product(self, other, inplace=True):
-        pass
-
-    @abstractmethod
-    def divide(self, other, inplace=True):
-        pass
-
+#    @abstractproperty
+#    def pdf(self):
+#        pass
+#
+#    @abstractproperty
+#    def variables(self):
+#        pass
+#
+#    @abstractmethod
+#    def assignment(self, *args, **kwargs):
+#        pass
+#
+#    @abstractmethod
+#    def copy(self):
+#        pass
+#    @abstractmethod
+#    def discretize(self, method, *args, **kwargs):
+#        pass
+#
+#    @abstractmethod
+#    def reduce(self, values, inplace=True):
+#        pass
+#
+#    @abstractmethod
+#    def marginalize(self, variables, inplace=True):
+#        pass
+#
+#    @abstractmethod
+#    def normalize(self, inplace=True):
+#         pass
+#
+#    @abstractmethod
+#    def product(self, other, inplace=True):
+#        pass
+#
+#    @abstractmethod
+#    def divide(self, other, inplace=True):
+#        pass
+    pass
diff --git a/pgmpy/models/LinearGaussianBayesianNetwork.py b/pgmpy/models/LinearGaussianBayesianNetwork.py
index 9db3c4b3..57c16667 100644
--- a/pgmpy/models/LinearGaussianBayesianNetwork.py
+++ b/pgmpy/models/LinearGaussianBayesianNetwork.py
@@ -5,7 +5,7 @@ import networkx as nx
 
 from pgmpy.models import BayesianModel
 from pgmpy.factors.continuous import LinearGaussianCPD
-from pgmpy.factors.distributions import JointGaussianDistribution
+from pgmpy.factors.distributions import GaussianDistribution
 
 
 class LinearGaussianBayesianNetwork(BayesianModel):
@@ -134,7 +134,7 @@ class LinearGaussianBayesianNetwork(BayesianModel):
 
         Returns
         -------
-        JointGaussianDistribution: An equivalent joint Gaussian
+        GaussianDistribution: An equivalent joint Gaussian
                                    distribution for the network.
 
         Reference
@@ -186,7 +186,7 @@ class LinearGaussianBayesianNetwork(BayesianModel):
                         [coeff * covariance[node_i_idx, variables.index(parent)]
                          for coeff, parent in zip(cpd_j.beta_vector, cpd_j.evidence)])
 
-        return JointGaussianDistribution(variables, mean, covariance)
+        return GaussianDistribution(variables, mean, covariance)
 
     def check_model(self):
         """
diff --git a/pgmpy/sampling/HMC.py b/pgmpy/sampling/HMC.py
index 1061ae18..d4c7a47e 100644
--- a/pgmpy/sampling/HMC.py
+++ b/pgmpy/sampling/HMC.py
@@ -37,7 +37,7 @@ class HamiltonianMC(object):
     Example:
     --------
     >>> from pgmpy.sampling import HamiltonianMC as HMC, LeapFrog, GradLogPDFGaussian
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+    >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([-3, 4])
     >>> covariance = np.array([[3, 0.7], [0.7, 5]])
@@ -215,7 +215,7 @@ class HamiltonianMC(object):
         Examples
         --------
         >>> from pgmpy.sampling import HamiltonianMC as HMC, GradLogPDFGaussian, ModifiedEuler
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, -1])
         >>> covariance = np.array([[1, 0.2], [0.2, 1]])
@@ -298,7 +298,7 @@ class HamiltonianMC(object):
         Examples
         --------
         >>> from pgmpy.sampling import HamiltonianMC as HMC, GradLogPDFGaussian as GLPG
-        >>> from pgmpy.factors import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([4, -1])
         >>> covariance = np.array([[3, 0.4], [0.4, 3]])
@@ -370,7 +370,7 @@ class HamiltonianMCDA(HamiltonianMC):
     Example:
     --------
     >>> from pgmpy.sampling import HamiltonianMCDA as HMCda, LeapFrog, GradLogPDFGaussian as GLPG
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+    >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([1, 2, 3])
     >>> covariance = np.array([[2, 0.4, 0.5], [0.4, 3, 0.6], [0.5, 0.6, 4]])
@@ -459,7 +459,7 @@ class HamiltonianMCDA(HamiltonianMC):
         Examples
         ---------
         >>> from pgmpy.sampling import HamiltonianMCDA as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, 0.7], [0.7, 3]])
@@ -548,7 +548,7 @@ class HamiltonianMCDA(HamiltonianMC):
         Examples
         --------
         >>> from pgmpy.sampling import HamiltonianMCDA as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, 0.7], [0.7, 3]])
diff --git a/pgmpy/sampling/NUTS.py b/pgmpy/sampling/NUTS.py
index c941c944..5cdeec8b 100644
--- a/pgmpy/sampling/NUTS.py
+++ b/pgmpy/sampling/NUTS.py
@@ -33,7 +33,7 @@ class NoUTurnSampler(HamiltonianMCDA):
     Example:
     --------
     >>> from pgmpy.sampling import NoUTurnSampler as NUTS, LeapFrog, GradLogPDFGaussian
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+    >>> from pgmpy.factors.continuous import GaussianDistribution 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]])
@@ -209,7 +209,7 @@ class NoUTurnSampler(HamiltonianMCDA):
         Examples
         ---------
         >>> from pgmpy.sampling import NoUTurnSampler as NUTS, GradLogPDFGaussian, LeapFrog
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution 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]])
@@ -273,7 +273,7 @@ class NoUTurnSampler(HamiltonianMCDA):
         Examples
         ---------
         >>> from pgmpy.sampling import NoUTurnSampler as NUTS, GradLogPDFGaussian
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([11, -6])
         >>> covariance = np.array([[0.7, 0.2], [0.2, 14]])
@@ -337,7 +337,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
     Example:
     --------
     >>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+    >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([-1, 12, -3])
     >>> covariance = np.array([[-2, 7, 2], [7, 14, 4], [2, 4, -1]])
@@ -499,7 +499,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
         Examples
         ---------
         >>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian, LeapFrog
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([10, -13])
         >>> covariance = np.array([[16, -3], [-3, 13]])
@@ -579,7 +579,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
         Examples
         --------
         >>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.continuous import GaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, -100])
         >>> covariance = np.array([[-12, 45], [45, -10]])
diff --git a/pgmpy/sampling/base.py b/pgmpy/sampling/base.py
index 5ff9f442..3296fc9c 100644
--- a/pgmpy/sampling/base.py
+++ b/pgmpy/sampling/base.py
@@ -27,7 +27,7 @@ class BaseGradLogPDF(object):
 
     Examples
     --------
-    >>> from pgmpy.factors import JointGaussianDistribution
+    >>> from pgmpy.factors import GaussianDistribution
     >>> from pgmpy.inference.continuous import BaseGradLogPDF
     >>> import numpy as np
     >>> class GradLogGaussian(BaseGradLogPDF):
@@ -41,7 +41,7 @@ class BaseGradLogPDF(object):
     ...         return grad, log_pdf
     >>> mean = np.array([1, 1])
     >>> covariance = np.array([[1, 0.2], [0.2, 7]])
-    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
     >>> dist_param = np.array([0.1, 0.9])
     >>> grad_logp, logp = GradLogGaussian(dist_param, model).get_gradient_log_pdf()
     >>> logp
@@ -79,11 +79,11 @@ class BaseGradLogPDF(object):
         --------
         >>> # Using implementation of GradLogPDFGaussian
         >>> from pgmpy.sampling.base import GradLogPDFGaussian
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution
+        >>> from pgmpy.factors.continuous import GaussianDistribution
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, -5], [-5, 2]])
-        >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+        >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
         >>> dist_param = np.array([0.6, 0.8])
         >>> grad_logp, logp = GradLogPDFGaussian(dist_param, model).get_gradient_log_pdf()
         >>> logp
@@ -98,23 +98,23 @@ class GradLogPDFGaussian(BaseGradLogPDF):
     """
     Class for finding gradient and gradient log of Joint Gaussian Distribution
     Inherits pgmpy.inference.base_continuous.BaseGradLogPDF
-    Here model must be an instance of JointGaussianDistribution
+    Here model must be an instance of GaussianDistribution
 
     Parameters
     ----------
     variable_assignments : A 1d array like object (numpy.ndarray or list)
         Vector representing values of variables at which we want to find gradient and log
 
-    model : An instance of pgmpy.models.JointGaussianDistribution
+    model : An instance of pgmpy.models.GaussianDistribution
 
     Example
     -------
     >>> from pgmpy.sampling import GradLogPDFGaussian
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution
+    >>> from pgmpy.factors.continuous import GaussianDistribution
     >>> import numpy as np
     >>> mean = np.array([3, 4])
     >>> covariance = np.array([[5, 4], [4, 5]])
-    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
     >>> dist_param = np.array([12, 21])
     >>> grad_logp, logp = GradLogPDFGaussian(dist_param, model).get_gradient_log_pdf()
     >>> logp
@@ -169,7 +169,7 @@ class BaseSimulateHamiltonianDynamics(object):
     Examples
     --------
     >>> from pgmpy.sampling import BaseSimulateHamiltonianDynamics
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution
+    >>> from pgmpy.factors.continuous import GaussianDistribution
     >>> from pgmpy.sampling import GradLogPDFGaussian
     >>> import numpy as np
     >>> # Class should initalize self.new_position, self.new_momentum and self.new_grad_logp
@@ -188,7 +188,7 @@ class BaseSimulateHamiltonianDynamics(object):
     >>> momentum = np.array([0, 0])
     >>> mean = np.array([0, 0])
     >>> covariance = np.eye(2)
-    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
     >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum,
     ...                                                 0.25, GradLogPDFGaussian).get_proposed_values()
     >>> new_pos
@@ -248,13 +248,13 @@ class BaseSimulateHamiltonianDynamics(object):
         -------
         >>> # Using implementation of ModifiedEuler
         >>> from pgmpy.inference.continuous import ModifiedEuler, GradLogPDFGaussian as GLPG
-        >>> from pgmpy.factors import JointGaussianDistribution
+        >>> from pgmpy.factors import GaussianDistribution
         >>> import numpy as np
         >>> pos = np.array([3, 4])
         >>> momentum = np.array([1, 1])
         >>> mean = np.array([-1, 1])
         >>> covariance = 3*np.eye(2)
-        >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+        >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
         >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum, 0.70, GLPG).get_proposed_values()
         >>> new_pos
         array([ 3.04666667,  4.21      ])
@@ -295,14 +295,14 @@ class LeapFrog(BaseSimulateHamiltonianDynamics):
 
     Example
     --------
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution
+    >>> from pgmpy.factors.continuous import GaussianDistribution
     >>> from pgmpy.sampling import LeapFrog, GradLogPDFGaussian as GLPG
     >>> import numpy as np
     >>> pos = np.array([2, 1])
     >>> momentum = np.array([7, 7])
     >>> mean = np.array([-5, 5])
     >>> covariance = np.array([[1, 2], [2, 1]])
-    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
     >>> new_pos, new_momentum, new_grad = LeapFrog(model, pos, momentum, 4.0, GLPG).get_proposed_values()
     >>> new_pos
     array([ 70., -19.])
@@ -366,14 +366,14 @@ class ModifiedEuler(BaseSimulateHamiltonianDynamics):
 
     Example
     --------
-    >>> from pgmpy.factors.continuous import JointGaussianDistribution
+    >>> from pgmpy.factors.continuous import GaussianDistribution
     >>> from pgmpy.sampling import GradLogPDFGaussian, ModifiedEuler
     >>> import numpy as np
     >>> pos = np.array([2, 1])
     >>> momentum = np.array([1, 1])
     >>> mean = np.array([0, 0])
     >>> covariance = np.eye(2)
-    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> model = GaussianDistribution(['x', 'y'], mean, covariance)
     >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum,
     ...                                                 0.25, GradLogPDFGaussian).get_proposed_values()
     >>> new_pos


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
index f8bd3918..5f71c6be 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
@@ -3,55 +3,55 @@ import unittest
 import numpy as np
 import numpy.testing as np_test
 
-from pgmpy.factors.distributions import JointGaussianDistribution as JGD
-from pgmpy.factors.continuous import CanonicalFactor
+from pgmpy.factors.distributions import GaussianDistribution as JGD
+from pgmpy.factors.continuous import CanonicalDistribution
 
 
 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)
+        phi = CanonicalDistribution(['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)
+        phi = CanonicalDistribution(['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)
+        phi = CanonicalDistribution(['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'],
+        self.assertRaises(ValueError, CanonicalDistribution, ['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'],
+        self.assertRaises(ValueError, CanonicalDistribution, ['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'],
+        self.assertRaises(ValueError, CanonicalDistribution, ['x1', 'x2', 'x3'],
                           np.array([[1.1, -1, 0], [0, -2, 4]]),
                           np.array([1, 2, 3]), 7)
-        self.assertRaises(ValueError, CanonicalFactor, ['x1', 'x3'],
+        self.assertRaises(ValueError, CanonicalDistribution, ['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.phi1 = CanonicalDistribution(['x1', 'x2', 'x3'],
+                                          np.array([[1.1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
+                                          np.array([[1], [4.7], [-1]]), -2)
+        self.phi2 = CanonicalDistribution(['x'], [[1]], [0], 1)
         self.phi3 = self.phi1.copy()
 
         self.gauss_phi1 = JGD(['x1', 'x2', 'x3'],
@@ -139,8 +139,8 @@ class TestJGDMethods(unittest.TestCase):
         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)
+        phi1 = self.phi1 * CanonicalDistribution(['x2', 'x4'], [[1, 2], [3, 4]], [0, 4.56], -6.78)
+        phi2 = self.phi1 / CanonicalDistribution(['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.0,  0.0],
diff --git a/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py b/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
index 62918756..c5b57deb 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
@@ -3,7 +3,7 @@ import unittest
 import numpy as np
 import numpy.testing as np_test
 
-from pgmpy.factors.distributions import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import GaussianDistribution as JGD
 
 
 class TestJGDInit(unittest.TestCase):
diff --git a/pgmpy/tests/test_readwrite/test_ProbModelXML.py b/pgmpy/tests/test_readwrite/test_ProbModelXML.py
index 39a8fb40..a94c01eb 100644
--- a/pgmpy/tests/test_readwrite/test_ProbModelXML.py
+++ b/pgmpy/tests/test_readwrite/test_ProbModelXML.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
+import os
 import unittest
 import warnings
 import json
@@ -1071,6 +1072,11 @@ class TestProbModelXMLWriter(unittest.TestCase):
         self.assertEqual(str(self.writer.__str__()[:-1]), str(etree.tostring(self.expected_xml)))
         self.assertEqual(str(data), str(etree.tostring(self.expected_xml).decode('utf-8')))
 
+    def tearDown(self):
+        try:
+            os.remove('test_xml.pgmx')
+        except OSError:
+            pass
 
 class TestProbModelXMLmethods(unittest.TestCase):
     def setUp(self):
diff --git a/pgmpy/tests/test_sampling/test_base_continuous.py b/pgmpy/tests/test_sampling/test_base_continuous.py
index b3220644..2d19585e 100644
--- a/pgmpy/tests/test_sampling/test_base_continuous.py
+++ b/pgmpy/tests/test_sampling/test_base_continuous.py
@@ -2,7 +2,7 @@ import unittest
 
 import numpy as np
 
-from pgmpy.factors.distributions import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import GaussianDistribution as JGD
 from pgmpy.sampling import LeapFrog, ModifiedEuler, GradLogPDFGaussian
 
 
diff --git a/pgmpy/tests/test_sampling/test_continuous_sampling.py b/pgmpy/tests/test_sampling/test_continuous_sampling.py
index de159ba5..a32b6d2b 100644
--- a/pgmpy/tests/test_sampling/test_continuous_sampling.py
+++ b/pgmpy/tests/test_sampling/test_continuous_sampling.py
@@ -2,7 +2,7 @@ import unittest
 
 import numpy as np
 
-from pgmpy.factors.distributions import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import GaussianDistribution as JGD
 from pgmpy.sampling import (HamiltonianMC as HMC, HamiltonianMCDA as HMCda, GradLogPDFGaussian, NoUTurnSampler as NUTS,
                             NoUTurnSamplerDA as NUTSda)