Repository: pgmpy/pgmpy
Pull request #862 (2017-03-24 16:10:00)

Custom distribution
Fixed issues in #837 

@ankurankan @yashu-seth Please review.


Patch:
diff --git a/pgmpy/factors/continuous/CanonicalFactor.py b/pgmpy/factors/continuous/CanonicalFactor.py
index 993696ec..a40296c9 100644
--- a/pgmpy/factors/continuous/CanonicalFactor.py
+++ b/pgmpy/factors/continuous/CanonicalFactor.py
@@ -5,7 +5,7 @@ from __future__ import division
 import numpy as np
 
 from pgmpy.factors.continuous import ContinuousFactor
-from pgmpy.factors.continuous import JointGaussianDistribution
+from pgmpy.factors.distributions import JointGaussianDistribution
 
 
 class CanonicalFactor(ContinuousFactor):
diff --git a/pgmpy/factors/continuous/__init__.py b/pgmpy/factors/continuous/__init__.py
index ff2bdb9d..bf76fdc3 100644
--- a/pgmpy/factors/continuous/__init__.py
+++ b/pgmpy/factors/continuous/__init__.py
@@ -1,15 +1,13 @@
 from .discretize import BaseDiscretizer, RoundingDiscretizer, UnbiasedDiscretizer
 from .ContinuousFactor import ContinuousFactor
-from .JointGaussianDistribution import JointGaussianDistribution
 from .CanonicalFactor import CanonicalFactor
 from .LinearGaussianCPD import LinearGaussianCPD
 
 
 __all__ = ['CanonicalFactor',
            'ContinuousFactor',
-           'JointGaussianDistribution',
            'LinearGaussianCPD'
            'BaseDiscretizer',
            'RoundingDiscretizer',
            'UnbiasedDiscretizer'
-           ]
\ No newline at end of file
+           ]
diff --git a/pgmpy/factors/distributions/CustomDistribution.py b/pgmpy/factors/distributions/CustomDistribution.py
new file mode 100644
index 00000000..864a13a3
--- /dev/null
+++ b/pgmpy/factors/distributions/CustomDistribution.py
@@ -0,0 +1,506 @@
+import numpy as np
+from scipy import integrate
+
+from .base import BaseDistribution
+
+class CustomDistribution(BaseDistribution):
+    def __init__(self, variables, distribution, *args, **kwargs):
+        """
+        Class for representing custom continuous distributions.
+
+        Parameters
+        ----------
+        variables: list or array-like
+            The variables for which the distribution is defined.
+
+        distribution: function
+            The probability density function of the distribution.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        # Two variable dirichlet distribution with alpha = (1, 2)
+        >>> def dirichlet_pdf(x, y):
+        ...     return (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> dirichlet_dist = CustomDistribution(variables=['x', 'y'],distribution=dirichlet_pdf)
+        >>> dirichlet_dist.variables
+        ['x', 'y']
+        """
+        if not isinstance(variables, (list, tuple, np.ndarray)):
+            raise TypeError(
+                "variables: Expected type: iterable, got: {type}".format(
+                    type=type(variables)))
+
+        if len(set(variables)) != len(variables):
+            raise ValueError("Multiple variables can't have the same name")
+
+        self._variables = list(variables)
+        self._pdf = distribution
+
+    @property
+    def pdf(self):
+        """
+        Returns the Probability Density Function of the distribution.
+
+        Returns
+        -------
+        function: The probability density function of the distribution
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        # Two variable dirichlet distribution with alpha = (1, 2)
+        >>> def dirichlet_pdf(x, y):
+        ...     return (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> dirichlet_dist = CustomDistribution(variables=['x', 'y'],
+        ...                                     distribution=dirichlet_pdf)
+        >>> dirichlet_dist.pdf()
+        <function __main__.diri_pdf>
+        """
+        return self._pdf
+
+    @property
+    def variables(self):
+        """
+        Returns the scope of the distribution.
+
+        Returns
+        -------
+        list: List of variables on which the distribution is defined.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        # Two variable dirichlet distribution with alpha = (1, 2)
+        >>> def dirichlet_pdf(x, y):
+        ...     return (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> dirichlet_dist = CustomDistribution(variables=['x', 'y'],
+        ...                                     distribution=dirichlet_pdf)
+        >>> dirichlet_dist.variables
+        ['x', 'y']
+        """
+        return self._variables
+
+    @variables.setter
+    def variables(self, value):
+        self._variables = value
+
+    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 CustomDistribution
+        >>> from scipy.stats import multivariate_normal
+        >>> normal_pdf = lambda x1, x2: multivariate_normal.pdf(
+        ...     x=(x1, x2), mean=[0, 0], cov=[[1, 0], [0, 1]])
+        >>> normal_dist = CustomDistribution(variables=['x1', 'x2'],
+        ...                                  distribution=normal_pdf)
+        >>> normal_dist.assignment(0, 0)
+        0.15915494309189535
+        0.15915494309189535
+        """
+        return self.pdf(*x)
+
+    def copy(self):
+        """
+        Returns a copy of the CustomDistribution instance.
+
+        Returns
+        -------
+        CustomDistribution object: copy of the instance
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        # Two variable drichlet ditribution with alpha = (1,2)
+        >>> def dirichlet_pdf(x, y):
+        ...     return (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> dirichlet_dist = CustomDistribution(variables=['x', 'y'],
+        ...                                     distribution=dirichlet_pdf)
+        >>> copy_dist = dirichlet_dist.copy()
+        >>> copy_dist.variables
+        ['x', 'y']
+        """
+        return CustomDistribution(self.variables, self.pdf)
+
+    # TODO: Discretize methods need to be fixed for this to work
+    def discretize(self, method, *args, **kwargs):
+        """
+        Discretizes the continuous distribution into discrete
+        probability masses using specified method.
+
+        Parameters
+        ----------
+        method: string, BaseDiscretizer instance
+            A Discretizer Class from pgmpy.factors.discretize
+
+        *args, **kwargs: values
+            The parameters to be given to the Discretizer Class.
+
+        Returns
+        -------
+        An n-D array or a DiscreteFactor object according to the discretiztion
+        method used.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.continuous import ContinuousFactor
+        >>> from pgmpy.factors.continuous import RoundingDiscretizer
+        >>> def dirichlet_pdf(x, y):
+        ...     return (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> dirichlet_factor = ContinuousFactor(['x', 'y'], dirichlet_pdf)
+        >>> dirichlet_factor.discretize(RoundingDiscretizer,
+        ...                             low=1, high=2, cardinality=5)
+        # TODO: finish this
+        """
+        super(CustomDistribution, self).discretize(method, *args, **kwargs)
+
+    def reduce(self, values, inplace=True):
+        """
+        Reduces the factor to the context of the given variable values.
+
+        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 CustomDistribution object.
+
+        Returns
+        -------
+        CustomDistribution or None:
+                    if inplace=True (default) returns None
+                    if inplace=False returns a new CustomDistribution instance.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from scipy.special import beta
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        >>> def custom_pdf(x, y, z):
+        ...     return z*(np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+        >>> custom_dist = CustomDistribution(['x', 'y', 'z'], custom_pdf)
+        >>> custom_dist.variables
+        ['x', 'y', 'z']
+        >>> custom_dist.assignment(1, 2, 3)
+        24.0
+
+        >>> custom_dist.reduce([('y', 2)])
+        >>> custom_dist.variables
+        ['x', 'z']
+        >>> custom_dist.assignment(1, 3)
+        24.0
+        """
+        if not isinstance(values, (list, tuple, np.ndarray)):
+            raise TypeError("variables: Expected type: iterable, "
+                            "got: {var_type}".format(var_type=type(values)))
+
+        for var, value in values:
+            if var not in self.variables:
+                raise ValueError("{var} not in scope.".format(var=var))
+
+        phi = self if inplace else self.copy()
+
+        var_to_remove = [var for var, value in values]
+        var_to_keep = [var for var in self.variables if var not in var_to_remove]
+
+        reduced_var_index = [(self.variables.index(var), value) for var, value in values]
+        pdf = self.pdf
+
+        def reduced_pdf(*args, **kwargs):
+            reduced_args = list(args)
+            reduced_kwargs = kwargs.copy()
+
+            if reduced_args:
+                for index, val in reduced_var_index:
+                    reduced_args.insert(index, val)
+            if reduced_kwargs:
+                for variable, val in values:
+                    reduced_kwargs[variable] = val
+            if reduced_args and reduced_kwargs:
+                reduced_args = [arg for arg in reduced_args if arg not in
+                                reduced_kwargs.values()]
+
+            return pdf(*reduced_args, **reduced_kwargs)
+
+        phi.variables = var_to_keep
+        phi._pdf = reduced_pdf
+
+        if not inplace:
+            return phi
+
+    def marginalize(self, variables, inplace=True):
+        """
+        Marginalize the distribution with respect to the given variables.
+
+        Parameters
+        ----------
+        variables: list, array-like
+            List of variables to be removed from the marginalized distribution.
+
+        inplace: boolean
+            If inplace=True it will modify the factor itself, else would return
+            a new CustomDistribution instance.
+
+        Returns
+        -------
+        Marginalized distribution or None:
+                if inplace=True (default) returns None
+                if inplace=False returns a new CustomDistribution instance.
+
+        Examples
+        --------
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        >>> from scipy.stats import multivariate_normal
+        >>> normal_pdf = lambda x1, x2: multivariate_normal.pdf(
+        ...                x=[x1, x2], mean=[0, 0], cov=[[1, 0], [0, 1]])
+        >>> normal_dist = CustomDistribution(variables=['x1', 'x2'],
+        ...                                  distribution=normal_pdf)
+        >>> normal_dist.variables
+        ['x1', 'x2']
+        >>> normal_dist.assignment(1, 1)
+        0.058549831524319168
+        >>> normal_dist.marginalize(['x2'])
+        >>> normal_dist.variables
+        ['x1']
+        >>> normal_dist.assignment(1)
+        0.24197072451914328
+        """
+        if len(variables) == 0:
+            raise ValueError("Shouldn't be calling marginalize over no variable.")
+
+        if not isinstance(variables, (list, tuple, np.ndarray)):
+            raise TypeError("variables: Expected type iterable, "
+                            "got: {var_type}".format(var_type=type(variables)))
+
+        for var in variables:
+            if var not in self.variables:
+                raise ValueError("{var} not in scope.".format(var=var))
+
+        phi = self if inplace else self.copy()
+
+        all_var = [var for var in self.variables]
+        var_to_keep = [var for var in self.variables if var not in variables]
+        reordered_var_index = [all_var.index(var) for var in variables + var_to_keep]
+        pdf = phi._pdf
+
+        # The arguments need to be reordered because integrate.nquad
+        # integrates the first n-arguments of the function passed.
+
+        def reordered_pdf(*args):
+            # ordered_args restores the original order as it was in self.variables
+            ordered_args = [args[reordered_var_index.index(index_id)] for index_id in range(len(all_var))]
+            return pdf(*ordered_args)
+
+        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.variables = var_to_keep
+
+        if not inplace:
+            return phi
+        
+    def normalize(self, inplace=True):
+        """
+        Normalizes the pdf of the distribution so that it
+        integrates to 1 over all the variables.
+
+        Parameters
+        ----------
+        inplace: boolean
+            If inplace=True it will modify the distribution itself, else would return
+            a new distribution.
+
+        Returns
+        -------
+        CustomDistribution or None:
+             if inplace=True (default) returns None
+             if inplace=False returns a new CustomDistribution instance.
+
+        Examples
+        --------
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        >>> from scipy.stats import multivariate_normal
+        >>> normal_pdf_x2 = lambda x1, x2: 2 * multivariate_normal.pdf(
+        ...                     x=[x1, x2], mean=[0, 0], cov=[[1, 0], [0, 1]])
+        >>> normal_dist_x2 = CustomDistribution(variables=['x1', 'x2'],
+        ...                                     distribution=normal_pdf_x2)
+        >>> normal_dist_x2.assignment(1, 1)
+        0.117099663049
+        >>> normal_dist = normal_dist_x2.normalize(inplace=False))
+        >>> normal_dist.assignment(1, 1)
+        0.0585498315243
+        """
+        phi = self if inplace else self.copy()
+        pdf = self.pdf
+
+        pdf_mod = integrate.nquad(pdf, [[-np.inf, np.inf] for var in self.variables])[0]
+
+        phi._pdf = lambda *args: pdf(*args) / pdf_mod
+
+        if not inplace:
+            return phi
+
+    def _operate(self, other, operation, inplace=True):
+        """
+        Gives the CustomDistribution operation (product or divide) with
+        the other distribution.
+
+        Parameters
+        ----------
+        other: CustomDistribution
+            The CustomDistribution to be multiplied.
+
+        operation: str
+            'product' for multiplication operation and 'divide' for
+            division operation.
+
+        inplace: boolean
+            If inplace=True it will modify the factor itself, else would return
+            a new distribution.
+
+        Returns
+        -------
+        CustomDistribution or None:
+                        if inplace=True (default) returns None
+                        if inplace=False returns a new `CustomDistribution` instance.
+
+        """
+        if not isinstance(other, CustomDistribution):
+            raise TypeError("CustomDistribution objects can only be multiplied "
+                            "or divided with another CustomDistribution  "
+                            "object. Got {other_type}, expected: "
+                            "CustomDistribution.".format(other_type=type(other)))
+
+        phi = self if inplace else self.copy()
+        pdf = self.pdf
+        self_var = [var for var in self.variables]
+
+        modified_pdf_var = self_var + [var for var in other.variables if var not in self_var]
+
+        def modified_pdf(*args):
+            self_pdf_args = list(args[:len(self_var)])
+            other_pdf_args = [args[modified_pdf_var.index(var)] for var in other.variables]
+
+            if operation == 'product':
+                return pdf(*self_pdf_args) * other._pdf(*other_pdf_args)
+            if operation == 'divide':
+                return pdf(*self_pdf_args) / other._pdf(*other_pdf_args)
+
+        phi.variables = modified_pdf_var
+        phi._pdf = modified_pdf
+
+        if not inplace:
+            return phi
+
+    def product(self, other, inplace=True):
+        """
+        Gives the CustomDistribution product with the other distribution.
+
+        Parameters
+        ----------
+        other: CustomDistribution
+            The CustomDistribution to be multiplied.
+
+        Returns
+        -------
+        CustomDistribution or None:
+                        if inplace=True (default) returns None
+                        if inplace=False returns a new `CustomDistribution` instance.
+
+        Example
+        -------
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        >>> from scipy.stats import multivariate_normal
+        >>> sn_pdf1 = lambda x: multivariate_normal.pdf(
+        ...                                 x=[x], mean=[0], cov=[[1]])
+        >>> sn_pdf2 = lambda x1,x2: multivariate_normal.pdf(
+        ...                     x=[x1, x2], mean=[0, 0], cov=[[1, 0], [0, 1]])
+        >>> sn1 = CustomDistribution(variables=['x2'], distribution=sn_pdf1)
+        >>> sn2 = CustomDistribution(variables=['x1', 'x2'],
+        ...                          distribution=sn_pdf2)
+
+        >>> sn3 = sn1.product(sn2, inplace=False)
+        >>> sn3.assignment(0, 0)
+        0.063493635934240983
+
+        >>> sn3 = sn1 * sn2
+        >>> sn3.assignment(0, 0)
+        0.063493635934240983
+        """
+        return self._operate(other, 'product', inplace)
+
+    def divide(self, other, inplace=True):
+        """
+        Gives the CustomDistribution divide with the other factor.
+
+        Parameters
+        ----------
+        other: ContinuousCustomDistributionFactor
+            The CustomDistribution to be multiplied.
+
+        Returns
+        -------
+        CustomDistribution or None:
+                        if inplace=True (default) returns None
+                        if inplace=False returns a new `CustomDistribution` instance.
+
+        Example
+        -------
+        >>> from pgmpy.factors.distributions import CustomDistribution
+        >>> from scipy.stats import multivariate_normal
+        >>> sn_pdf1 = lambda x: multivariate_normal.pdf(
+        ...                                     x=[x], mean=[0], cov=[[1]])
+        >>> sn_pdf2 = lambda x1, x2: multivariate_normal.pdf(
+        ...                 x=[x1, x2], mean=[0, 0], cov=[[1, 0], [0, 1]])
+        >>> sn1 = CustomDistribution(variables=['x2'], distribution=sn_pdf1)
+        >>> sn2 = CustomDistribution(variables=['x1', 'x2'],
+        ...                          distribution=sn_pdf2)
+
+        >>> sn3 = sn2.divide(sn1, inplace=False)
+        >>> sn3.assignment(0, 0)
+        0.3989422804014327
+
+        >>> sn3 = sn2 / sn1
+        >>> sn3.assignment(0, 0)
+        0.3989422804014327
+        """
+        if set(other.variables) - set(self.variables):
+            raise ValueError("Scope of divisor should be a subset of dividend")
+
+        return self._operate(other, 'divide', 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/JointGaussianDistribution.py b/pgmpy/factors/distributions/JointGaussianDistribution.py
similarity index 93%
rename from pgmpy/factors/continuous/JointGaussianDistribution.py
rename to pgmpy/factors/distributions/JointGaussianDistribution.py
index 4aaabc7c..6a14f89c 100644
--- a/pgmpy/factors/continuous/JointGaussianDistribution.py
+++ b/pgmpy/factors/distributions/JointGaussianDistribution.py
@@ -5,15 +5,17 @@ from __future__ import division
 import numpy as np
 from scipy.stats import multivariate_normal
 
-from pgmpy.factors.continuous import ContinuousFactor
+from pgmpy.factors.distributions import CustomDistribution
 
 
-class JointGaussianDistribution(ContinuousFactor):
+class JointGaussianDistribution(CustomDistribution):
     u"""
     In its most common representation, a multivariate Gaussian distribution
     over X1...........Xn is characterized by an n-dimensional mean vector μ,
     and a symmetric n x n covariance matrix Σ.
     This is the base class for its representation.
+
+    #TODO: Add more examples after refactoring ContinuousFactor class
     """
     def __init__(self, variables, mean, covariance):
         """
@@ -31,7 +33,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]]))
         >>> dis.variables
@@ -77,7 +79,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]]))
         >>> dis.precision_matrix
@@ -113,7 +115,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]]))
         >>> dis.variables
@@ -187,7 +189,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]]))
         >>> dis.variables
@@ -264,7 +266,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions import JointGaussianDistribution as JGD
         >>> gauss_dis = JGD(['x1', 'x2', 'x3'], np.array([[1], [-3], [4]]),
         ...                 np.array([[4, 2, -2], [2, 5, -5], [-2, -5, 8]]))
         >>> copy_dis = gauss_dis.copy()
@@ -283,7 +285,7 @@ class JointGaussianDistribution(ContinuousFactor):
                 [-0.125     ,  0.58333333,  0.33333333],
                 [ 0.        ,  0.33333333,  0.33333333]])
         """
-        copy_distribution = JointGaussianDistribution(self.scope(), self.mean.copy(),
+        copy_distribution = JointGaussianDistribution(self.variables, self.mean.copy(),
                                                       self.covariance.copy())
         if self._precision_matrix is not None:
             copy_distribution._precision_matrix = self._precision_matrix.copy()
@@ -313,7 +315,7 @@ class JointGaussianDistribution(ContinuousFactor):
         -------
 
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]]))
         >>> phi = dis.to_canonical_factor()
@@ -342,7 +344,7 @@ class JointGaussianDistribution(ContinuousFactor):
         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.scope(), K, h, g)
+        return CanonicalFactor(self.variables, K, h, g)
 
     def _operate(self, other, operation, inplace=True):
         """
@@ -367,7 +369,7 @@ class JointGaussianDistribution(ContinuousFactor):
         Examples
         --------
         >>> import numpy as np
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors.distributions 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]])
diff --git a/pgmpy/factors/distributions/__init__.py b/pgmpy/factors/distributions/__init__.py
new file mode 100644
index 00000000..ce6d8c90
--- /dev/null
+++ b/pgmpy/factors/distributions/__init__.py
@@ -0,0 +1,9 @@
+from .base import BaseDistribution
+from .CustomDistribution import CustomDistribution
+from .JointGaussianDistribution import JointGaussianDistribution
+
+
+__all__ = ['BaseDistribution',
+           'CustomDistribution',
+           'JointGaussianDistribution'
+           ]
diff --git a/pgmpy/factors/distributions/base.py b/pgmpy/factors/distributions/base.py
new file mode 100644
index 00000000..f85ba5d1
--- /dev/null
+++ b/pgmpy/factors/distributions/base.py
@@ -0,0 +1,43 @@
+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
+
diff --git a/pgmpy/models/LinearGaussianBayesianNetwork.py b/pgmpy/models/LinearGaussianBayesianNetwork.py
index fab63cff..9db3c4b3 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.continuous import JointGaussianDistribution
+from pgmpy.factors.distributions import JointGaussianDistribution
 
 
 class LinearGaussianBayesianNetwork(BayesianModel):


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 b1a5cef0..f8bd3918 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
@@ -3,7 +3,7 @@ import unittest
 import numpy as np
 import numpy.testing as np_test
 
-from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import JointGaussianDistribution as JGD
 from pgmpy.factors.continuous import CanonicalFactor
 
 
diff --git a/pgmpy/tests/test_factors/test_continuous/test_CustomDistribution.py b/pgmpy/tests/test_factors/test_continuous/test_CustomDistribution.py
new file mode 100644
index 00000000..9ab56839
--- /dev/null
+++ b/pgmpy/tests/test_factors/test_continuous/test_CustomDistribution.py
@@ -0,0 +1,436 @@
+import unittest
+
+import numpy as np
+import numpy.testing as np_test
+from scipy.special import beta
+from scipy.stats import multivariate_normal
+
+from pgmpy.factors.distributions import CustomDistribution
+
+
+class TestCustomDistribution(unittest.TestCase):
+    def pdf1(self, x, y):
+        return np.power(x, 1) * np.power(y, 2) / beta(x, y)
+
+    def pdf2(self, *args):
+        return multivariate_normal.pdf(args, [0, 0], [[1, 0], [0, 1]])
+
+    def pdf3(self, x, y, z):
+        return z * (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+
+    def test_class_init(self):
+        phi1 = CustomDistribution(['x', 'y'], self.pdf1)
+        self.assertEqual(phi1.variables, ['x', 'y'])
+        self.assertEqual(phi1._pdf, self.pdf1)
+
+        phi2 = CustomDistribution(['x1', 'x2'], self.pdf2)
+        self.assertEqual(phi2.variables, ['x1', 'x2'])
+        self.assertEqual(phi2._pdf, self.pdf2)
+
+        phi3 = CustomDistribution(['x', 'y', 'z'], self.pdf3)
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        self.assertEqual(phi3._pdf, self.pdf3)
+
+    def test_class_init_typeerror(self):
+        self.assertRaises(TypeError, CustomDistribution, 'x y', self.pdf1)
+        self.assertRaises(TypeError, CustomDistribution, 'x', self.pdf1)
+
+        self.assertRaises(TypeError, CustomDistribution, 'x1 x2', self.pdf2)
+        self.assertRaises(TypeError, CustomDistribution, 'x1', self.pdf1)
+
+        self.assertRaises(TypeError, CustomDistribution, 'x y z', self.pdf3)
+        self.assertRaises(TypeError, CustomDistribution, 'x', self.pdf3)
+
+        self.assertRaises(TypeError, CustomDistribution, set(['x', 'y']), self.pdf1)
+        self.assertRaises(TypeError, CustomDistribution, {'x': 1, 'y': 2}, self.pdf1)
+
+        self.assertRaises(TypeError, CustomDistribution, set(['x1', 'x2']), self.pdf2)
+        self.assertRaises(TypeError, CustomDistribution, {'x1': 1, 'x2': 2}, self.pdf1)
+
+        self.assertRaises(TypeError, CustomDistribution, set(['x', 'y', 'z']), self.pdf3)
+        self.assertRaises(TypeError, CustomDistribution, {'x': 1, 'y': 2, 'z': 3}, self.pdf3)
+
+    def test_class_init_valueerror(self):
+        self.assertRaises(ValueError, CustomDistribution, ['x', 'x'], self.pdf1)
+        self.assertRaises(ValueError, CustomDistribution, ['x', 'y', 'y'], self.pdf1)
+
+        self.assertRaises(ValueError, CustomDistribution, ['x1', 'x1'], self.pdf2)
+        self.assertRaises(ValueError, CustomDistribution, ['x1', 'x2', 'x2'], self.pdf2)
+
+        self.assertRaises(ValueError, CustomDistribution, ['x', 'x'], self.pdf1)
+        self.assertRaises(ValueError, CustomDistribution, ['x', 'y', 'y', 'z', 'z'], self.pdf1)
+
+
+class TestCustomDistributionMethods(unittest.TestCase):
+    def pdf1(self, x, y):
+        return np.power(x, 1) * np.power(y, 2) / beta(x, y)
+
+    def pdf2(self, x1, x2):
+        return multivariate_normal.pdf([x1, x2], [0, 0], [[1, 0], [0, 1]])
+
+    def pdf3(self, x, y, z):
+        return z * (np.power(x, 1) * np.power(y, 2)) / beta(x, y)
+
+    def pdf4(self, x1, x2, x3):
+        return multivariate_normal.pdf([x1, x2, x3], [0, 0, 0],
+                                       [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+
+    def setUp(self):
+        self.phi1 = CustomDistribution(['x', 'y'], self.pdf1)
+        self.phi2 = CustomDistribution(['x1', 'x2'], self.pdf2)
+        self.phi3 = CustomDistribution(['x', 'y', 'z'], self.pdf3)
+        self.phi4 = CustomDistribution(['x1', 'x2', 'x3'], self.pdf4)
+
+    def test_variables(self):
+        self.assertEqual(self.phi1.variables, self.phi1._variables)
+        self.assertEqual(self.phi2.variables, self.phi2._variables)
+        self.assertEqual(self.phi3.variables, self.phi3._variables)
+
+    def test_assignment(self):
+        self.assertEqual(self.phi1.assignment(1.212, 2), self.pdf1(1.212, 2))
+        self.assertEqual(self.phi2.assignment(1, -2.231), self.pdf2(1, -2.231))
+        self.assertEqual(self.phi3.assignment(1.212, 2.213, -3), self.pdf3(1.212, 2.213, -3))
+
+    def test_reduce(self):
+        phi1 = self.phi1.copy()
+        phi1.reduce([('x', 1)])
+        reduced_pdf1 = lambda y: (np.power(1, 1) * np.power(y, 2))/beta(1, y)
+        self.assertEqual(phi1.variables, ['y'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi1._pdf(inp), reduced_pdf1(inp))
+            self.assertEqual(phi1._pdf(y=inp), reduced_pdf1(inp))
+
+        phi1 = self.phi1.reduce([('x', 1)], inplace=False)
+        self.assertEqual(phi1.variables, ['y'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi1._pdf(inp), reduced_pdf1(inp))
+            self.assertEqual(phi1._pdf(y=inp), reduced_pdf1(inp))
+
+        phi2 = self.phi2.copy()
+        phi2.reduce([('x2', 7.213)])
+        reduced_pdf2 = lambda x1: multivariate_normal.pdf([x1, 7.213], [0, 0], [[1, 0], [0, 1]])
+        self.assertEqual(phi2.variables, ['x1'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi2._pdf(inp), reduced_pdf2(inp))
+            self.assertEqual(phi2._pdf(x1=inp), reduced_pdf2(inp))
+
+        phi2 = self.phi2.reduce([('x2', 7.213)], inplace=False)
+        self.assertEqual(phi2.variables, ['x1'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi2._pdf(inp), reduced_pdf2(inp))
+            self.assertEqual(phi2._pdf(x1=inp), reduced_pdf2(inp))
+
+        phi3 = self.phi3.copy()
+        phi3.reduce([('y', 0.112), ('z', 23)])
+        reduced_pdf4 = lambda x: 23*(np.power(x, 1)*np.power(0.112, 2))/beta(x, 0.112)
+        self.assertEqual(phi3.variables, ['x'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi3._pdf(inp), reduced_pdf4(inp))
+            self.assertEqual(phi3._pdf(x=inp), reduced_pdf4(inp))
+
+        phi3 = self.phi3.copy()
+        phi3.reduce([('y', 0.112)])
+        reduced_pdf3 = lambda x, z: z*(np.power(x, 1)*np.power(0.112, 2))/beta(x, 0.112)
+        self.assertEqual(phi3.variables, ['x', 'z'])
+        for inp in np.random.rand(4, 2):
+            self.assertEqual(phi3._pdf(inp[0], inp[1]),
+                             reduced_pdf3(inp[0], inp[1]))
+            self.assertEqual(phi3._pdf(x=inp[0], z=inp[1]),
+                             reduced_pdf3(inp[0], inp[1]))
+
+        phi3 = self.phi3.reduce([('y', 0.112)], inplace=False)
+        self.assertEqual(phi3.variables, ['x', 'z'])
+        for inp in np.random.rand(4, 2):
+            self.assertEqual(phi3._pdf(inp[0], inp[1]),
+                             reduced_pdf3(inp[0], inp[1]))
+            self.assertEqual(phi3._pdf(x=inp[0], z=inp[1]),
+                             reduced_pdf3(inp[0], inp[1]))
+            self.assertEqual(phi3._pdf(inp[0], z=inp[1]),
+                             reduced_pdf3(inp[0], inp[1]))
+
+        phi3 = self.phi3.reduce([('y', 0.112), ('z', 23)], inplace=False)
+        self.assertEqual(phi3.variables, ['x'])
+        for inp in np.random.rand(4):
+            self.assertEqual(phi3._pdf(inp), reduced_pdf4(inp))
+            self.assertEqual(phi3._pdf(x=inp), reduced_pdf4(inp))
+
+    def test_reduce_error(self):
+        self.assertRaises(TypeError, self.phi1.reduce, 'x1')
+        self.assertRaises(TypeError, self.phi1.reduce, set(['x', 'y']))
+        self.assertRaises(TypeError, self.phi1.reduce, {'x': 1, 'y': 1})
+
+        self.assertRaises(TypeError, self.phi4.reduce, 'x4')
+        self.assertRaises(TypeError, self.phi4.reduce, set(['x1', 'x2', 'x3']))
+        self.assertRaises(TypeError, self.phi4.reduce, {'x1': 1, 'x2': 1, 'x3': 1})
+
+        self.assertRaises(ValueError, self.phi1.reduce, [('z', 3)])
+        self.assertRaises(ValueError, self.phi1.reduce, [('x', 0), ('y', 1), ('z', 4)])
+
+        self.assertRaises(ValueError, self.phi4.reduce, [('x4', 7)])
+        self.assertRaises(ValueError, self.phi4.reduce, [('x1', 1), ('x2', 2), ('x3', 3), ('x4', 4)])
+
+    def test_marginalize(self):
+        phi2 = self.phi2.copy()
+        phi2.marginalize(['x2'])
+        self.assertEqual(phi2.variables, ['x1'])
+        for inp in np.random.rand(4):
+            np_test.assert_almost_equal(phi2._pdf(inp),
+                                        multivariate_normal.pdf([inp], [0], [[1]]))
+
+        phi2 = self.phi2.marginalize(['x2'], inplace=False)
+        self.assertEqual(phi2.variables, ['x1'])
+        for inp in np.random.rand(4):
+            np_test.assert_almost_equal(phi2._pdf(inp),
+                                        multivariate_normal.pdf([inp], [0], [[1]]))
+
+        phi4 = self.phi4.copy()
+        phi4.marginalize(['x2'])
+
+        self.assertEqual(phi4.variables, ['x1', 'x3'])
+        for inp in np.random.rand(4, 2):
+            np_test.assert_almost_equal(
+                phi4._pdf(inp[0], inp[1]), multivariate_normal.pdf(
+                    [inp[0], inp[1]], [0, 0], [[1, 0], [0, 1]]))
+
+        phi4.marginalize(['x3'])
+        self.assertEqual(phi4.variables, ['x1'])
+        for inp in np.random.rand(1):
+            np_test.assert_almost_equal(phi4._pdf(inp),
+                                        multivariate_normal.pdf([inp], [0], [[1]]))
+
+        phi4 = self.phi4.marginalize(['x2'], inplace=False)
+        self.assertEqual(phi4.variables, ['x1', 'x3'])
+        for inp in np.random.rand(4, 2):
+            np_test.assert_almost_equal(phi4._pdf(inp[0], inp[1]),
+                                        multivariate_normal.pdf([inp[0], inp[1]], [0, 0], [[1, 0], [0, 1]]))
+
+        phi4 = phi4.marginalize(['x3'], inplace=False)
+        self.assertEqual(phi4.variables, ['x1'])
+        for inp in np.random.rand(1):
+            np_test.assert_almost_equal(phi4._pdf(inp),
+                                        multivariate_normal.pdf([inp], [0], [[1]]))
+
+    def test_marginalize_error(self):
+        self.assertRaises(TypeError, self.phi1.marginalize, 'x1')
+        self.assertRaises(TypeError, self.phi1.marginalize, set(['x', 'y']))
+        self.assertRaises(TypeError, self.phi1.marginalize, {'x': 1, 'y': 1})
+
+        self.assertRaises(TypeError, self.phi4.marginalize, 'x4')
+        self.assertRaises(TypeError, self.phi4.marginalize, set(['x1', 'x2', 'x3']))
+        self.assertRaises(TypeError, self.phi4.marginalize, {'x1': 1, 'x2': 1, 'x3': 1})
+
+        self.assertRaises(ValueError, self.phi1.marginalize, ['z'])
+        self.assertRaises(ValueError, self.phi1.marginalize, ['x', 'y', 'z'])
+
+        self.assertRaises(ValueError, self.phi4.marginalize, ['x4'])
+        self.assertRaises(ValueError, self.phi4.marginalize, ['x1', 'x2', 'x3', 'x4'])
+
+    def test_normalize(self):
+        def pdf2(x1, x2):
+            return 2 * self.pdf2(x1, x2)
+
+        phi2 = CustomDistribution(['x1', 'x2'], pdf2)
+        phi4 = phi2.copy()
+
+        phi4.normalize()
+        self.assertEqual(phi4.variables, phi2.variables)
+        for inp in np.random.rand(1, 2):
+            np_test.assert_almost_equal(phi4._pdf(inp[0], inp[1]),
+                                        self.pdf2(inp[0], inp[1]))
+
+        phi4 = phi2.normalize(inplace=False)
+        self.assertEqual(phi4.variables, phi4.variables)
+        for inp in np.random.rand(1, 2):
+            np_test.assert_almost_equal(phi4._pdf(inp[0], inp[1]),
+                                        self.pdf2(inp[0], inp[1]))
+
+    def test_operate(self):
+        phi1 = self.phi1.copy()
+        phi1._operate(self.phi2, 'product')
+        self.assertEqual(phi1.variables, ['x', 'y', 'x1', 'x2'])
+        for inp in np.random.rand(4, 4):
+            self.assertEqual(phi1._pdf(*inp), self.phi1._pdf(inp[0], inp[1])
+                             * self.phi2._pdf(inp[2], inp[3]))
+
+        phi1 = self.phi1._operate(self.phi2, 'product', inplace=False)
+        self.assertEqual(phi1.variables, ['x', 'y', 'x1', 'x2'])
+        for inp in np.random.rand(4, 4):
+            self.assertEqual(phi1._pdf(*inp), self.phi1._pdf(inp[0], inp[1])
+                             * self.phi2._pdf(inp[2], inp[3]))
+
+        phi1 = self.phi1 * self.phi2
+        self.assertEqual(phi1.variables, ['x', 'y', 'x1', 'x2'])
+        for inp in np.random.rand(4, 4):
+            self.assertEqual(phi1._pdf(*inp), self.phi1._pdf(inp[0], inp[1])
+                             * self.phi2._pdf(inp[2], inp[3]))
+
+        phi3 = self.phi3.copy()
+        phi3._operate(self.phi1, 'product')
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) *
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi3 = self.phi3._operate(self.phi1, 'product', inplace=False)
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) *
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi3 = self.phi3 * self.phi1
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) *
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi3 = self.phi3.copy()
+        phi3._operate(self.phi1, 'divide')
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) /
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi3 = self.phi3._operate(self.phi1, 'divide', inplace=False)
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) /
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi3 = self.phi3 / self.phi1
+        self.assertEqual(phi3.variables, ['x', 'y', 'z'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi3._pdf(*inp), self.phi3._pdf(*inp) /
+                             self.phi1._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4.copy()
+        phi4._operate(self.phi2, 'product')
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) *
+                             self.phi2._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4._operate(self.phi2, 'product', inplace=False)
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) *
+                             self.phi2._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4 * self.phi2
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) *
+                             self.phi2._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4.copy()
+        phi4._operate(self.phi2, 'divide')
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) /
+                             self.phi2._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4._operate(self.phi2, 'divide', inplace=False)
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) /
+                             self.phi2._pdf(inp[0], inp[1]))
+
+        phi4 = self.phi4 / self.phi2
+        self.assertEqual(phi4.variables, ['x1', 'x2', 'x3'])
+        for inp in np.random.rand(4, 3):
+            self.assertEqual(phi4._pdf(*inp), self.phi4._pdf(*inp) /
+                             self.phi2._pdf(inp[0], inp[1]))
+
+    def test_operate_error(self):
+        self.assertRaises(TypeError, self.phi1._operate, 1, 'product')
+        self.assertRaises(TypeError, self.phi1._operate, 1, 'divide')
+        self.assertRaises(TypeError, self.phi1._operate, '1', 'product')
+        self.assertRaises(TypeError, self.phi1._operate, '1', 'divide')
+        self.assertRaises(TypeError, self.phi1._operate,
+                          self.phi2._pdf, 'product')
+        self.assertRaises(TypeError, self.phi1._operate,
+                          self.phi2._pdf, 'divide')
+        self.assertRaises(TypeError, self.phi1._operate, [1], 'product')
+        self.assertRaises(TypeError, self.phi1._operate, [1], 'divide')
+
+        self.assertRaises(TypeError, self.phi4._operate, 1, 'product')
+        self.assertRaises(TypeError, self.phi4._operate, 1, 'divide')
+        self.assertRaises(TypeError, self.phi4._operate, '1', 'product')
+        self.assertRaises(TypeError, self.phi4._operate, '1', 'divide')
+        self.assertRaises(TypeError, self.phi4._operate,
+                          self.phi2._pdf, 'product')
+        self.assertRaises(TypeError, self.phi4._operate,
+                          self.phi2._pdf, 'divide')
+        self.assertRaises(TypeError, self.phi4._operate, [1], 'product')
+        self.assertRaises(TypeError, self.phi4._operate, [1], 'divide')
+
+        self.assertRaises(TypeError, self.phi1._operate, 1, 'product', False)
+        self.assertRaises(TypeError, self.phi1._operate, 1, 'divide', False)
+        self.assertRaises(TypeError, self.phi1._operate, '1', 'product', False)
+        self.assertRaises(TypeError, self.phi1._operate, '1', 'divide', False)
+        self.assertRaises(TypeError, self.phi1._operate,
+                          self.phi2._pdf, 'product', False)
+        self.assertRaises(TypeError, self.phi1._operate,
+                          self.phi2._pdf, 'divide', False)
+        self.assertRaises(TypeError, self.phi1._operate, [1], 'product', False)
+        self.assertRaises(TypeError, self.phi1._operate, [1], 'divide', False)
+
+        self.assertRaises(TypeError, self.phi4._operate, 1, 'product', False)
+        self.assertRaises(TypeError, self.phi4._operate, 1, 'divide', False)
+        self.assertRaises(TypeError, self.phi4._operate, '1', 'product', False)
+        self.assertRaises(TypeError, self.phi4._operate, '1', 'divide', False)
+        self.assertRaises(TypeError, self.phi4._operate,
+                          self.phi2._pdf, 'product', False)
+        self.assertRaises(TypeError, self.phi4._operate,
+                          self.phi2._pdf, 'divide', False)
+        self.assertRaises(TypeError, self.phi4._operate, [1], 'product', False)
+        self.assertRaises(TypeError, self.phi4._operate, [1], 'divide', False)
+
+        self.assertRaises(ValueError, self.phi1.__truediv__, self.phi2)
+        self.assertRaises(ValueError, self.phi1.__truediv__, self.phi3)
+        self.assertRaises(ValueError, self.phi1.__truediv__, self.phi4)
+        self.assertRaises(ValueError, self.phi2.__truediv__, self.phi3)
+        self.assertRaises(ValueError, self.phi2.__truediv__, self.phi4)
+
+    def test_copy(self):
+        copy1 = self.phi1.copy()
+        copy2 = self.phi3.copy()
+
+        copy4 = copy1.copy()
+        copy5 = copy2.copy()
+
+        self.assertEqual(copy1.variables, copy4.variables)
+        self.assertEqual(copy1._pdf, copy4._pdf)
+        self.assertEqual(copy2.variables, copy5.variables)
+        self.assertEqual(copy2._pdf, copy5._pdf)
+
+        copy1.variables = ['A', 'B']
+        self.assertEqual(copy4.variables, self.phi1.variables)
+
+        def pdf(a, b):
+            return (a + b) / (a * a + b * b)
+        copy1._pdf = pdf
+        copy1_pdf = pdf
+        self.assertEqual(copy4._pdf, self.phi1._pdf)
+        copy4.variables = ['X', 'Y']
+        self.assertEqual(copy1.variables, ['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]))
+
+        copy2.reduce([('x', 7.7)])
+
+        def reduced_pdf(y, z):
+            return z*(np.power(7.7, 1) * np.power(y, 2)) / beta(7.7, y)
+        self.assertEqual(copy5.variables, self.phi3.variables)
+        self.assertEqual(copy5._pdf, self.phi3._pdf)
+        copy5.reduce([('x', 11), ('z', 13)])
+        self.assertEqual(copy2.variables, ['y', 'z'])
+        for inp in np.random.rand(4, 2):
+            self.assertEqual(copy2._pdf(inp[0], inp[1]),
+                             reduced_pdf(inp[0], inp[1]))
+
+    def tearDown(self):
+        del self.phi1
+        del self.phi2
+        del self.phi3
\ No newline at end of file
diff --git a/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py b/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
index d899334d..62918756 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.continuous import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import JointGaussianDistribution as JGD
 
 
 class TestJGDInit(unittest.TestCase):
diff --git a/pgmpy/tests/test_sampling/test_base_continuous.py b/pgmpy/tests/test_sampling/test_base_continuous.py
index 29392d9e..b3220644 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.continuous import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import JointGaussianDistribution 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 bb115be2..de159ba5 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.continuous import JointGaussianDistribution as JGD
+from pgmpy.factors.distributions import JointGaussianDistribution as JGD
 from pgmpy.sampling import (HamiltonianMC as HMC, HamiltonianMCDA as HMCda, GradLogPDFGaussian, NoUTurnSampler as NUTS,
                             NoUTurnSamplerDA as NUTSda)