Repository: pgmpy/pgmpy
Pull request #711 (2016-07-20 13:51:11)

Accomodate changes from PR#693


Patch:
diff --git a/pgmpy/inference/continuous/base.py b/pgmpy/inference/continuous/base.py
index 074dea8b..ad7fc5f7 100644
--- a/pgmpy/inference/continuous/base.py
+++ b/pgmpy/inference/continuous/base.py
@@ -21,7 +21,7 @@ class BaseGradLogPDF(object):
 
     Examples
     --------
-    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.factors import JointGaussianDistribution
     >>> from pgmpy.inference.continuous import BaseGradLogPDF
     >>> import numpy as np
     >>> class GradLogGaussian(BaseGradLogPDF):
@@ -73,7 +73,7 @@ class BaseGradLogPDF(object):
         --------
         >>> # Using implementation of GradLogPDFGaussian
         >>> from pgmpy.inference.continuous import GradLogPDFGaussian
-        >>> from pgmpy.models import JointGaussianDistribution
+        >>> from pgmpy.factors import JointGaussianDistribution
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, -5], [-5, 2]])
@@ -104,7 +104,7 @@ class GradLogPDFGaussian(BaseGradLogPDF):
     Example
     -------
     >>> from pgmpy.inference.continuous import GradLogPDFGaussian
-    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.factors import JointGaussianDistribution
     >>> import numpy as np
     >>> mean = np.array([3, 4])
     >>> covariance = np.array([[5, 4], [4, 5]])
@@ -153,9 +153,8 @@ class BaseSimulateHamiltonianDynamics(object):
     stepsize: Float
         stepsize for the simulating dynamics
 
-    grad_log_pdf : A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
+    grad_log_pdf : A subclass of pgmpy.inference.continuous.BaseGradLogPDF
         A class for finding gradient log and log of distribution
-        If none, then model.get_gradient_log_pdf will be used
 
     grad_log_position: A 1d array like object, defaults to None
         Vector representing gradient log at given position
@@ -164,7 +163,7 @@ class BaseSimulateHamiltonianDynamics(object):
     Examples
     --------
     >>> from pgmpy.inference.continuous import BaseSimulateHamiltonianDynamics
-    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.factors import JointGaussianDistribution
     >>> from pgmpy.inference.continuous import GradLogPDFGaussian
     >>> import numpy as np
     >>> # Class should initalize self.new_position, self.new_momentum and self.new_grad_logp
@@ -177,7 +176,7 @@ class BaseSimulateHamiltonianDynamics(object):
     ...     def _get_proposed_values(self):
     ...         momentum_bar = self.momentum + self.stepsize * self.grad_log_position
     ...         position_bar = self.position + self.stepsize * momentum_bar
-    ...         grad_log_position, _ = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+    ...         grad_log_position, _ = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
     ...         return position_bar, momentum_bar, grad_log_position
     >>> pos = np.array([1, 2])
     >>> momentum = np.array([0, 0])
@@ -194,13 +193,13 @@ class BaseSimulateHamiltonianDynamics(object):
     array([-0.9375, -1.875])
     """
 
-    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, grad_log_position=None):
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf, grad_log_position=None):
 
         position = _check_1d_array_object(position, 'position')
 
         momentum = _check_1d_array_object(momentum, 'momentum')
 
-        if grad_log_pdf is not None and not issubclass(grad_log_pdf, BaseGradLogPDF):
+        if not issubclass(grad_log_pdf, BaseGradLogPDF):
             raise TypeError("grad_log_pdf must be an instance" +
                             " of pgmpy.inference.continuous.base.BaseGradLogPDF")
 
@@ -208,7 +207,7 @@ class BaseSimulateHamiltonianDynamics(object):
         _check_length_equal(position, model.variables, 'position', 'model.variables')
 
         if grad_log_position is None:
-            grad_log_position, _ = model.get_gradient_log_pdf(position, grad_log_pdf)
+            grad_log_position, _ = grad_log_pdf(position, model).get_gradient_log_pdf()
 
         else:
             grad_log_positon = _check_1d_array_object(grad_log_position, 'grad_log_position')
@@ -242,15 +241,15 @@ class BaseSimulateHamiltonianDynamics(object):
         Example
         -------
         >>> # Using implementation of ModifiedEuler
-        >>> from pgmpy.inference.continuous import ModifiedEuler
-        >>> from pgmpy.models import JointGaussianDistribution
+        >>> from pgmpy.inference.continuous import ModifiedEuler, GradLogPDFGaussian as GLPG
+        >>> from pgmpy.factors import JointGaussianDistribution
         >>> 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)
-        >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum, 0.70).get_proposed_values()
+        >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum, 0.70, GLPG).get_proposed_values()
         >>> new_pos
         array([ 3.04666667,  4.21      ])
         >>> new_momentum
@@ -290,15 +289,15 @@ class LeapFrog(BaseSimulateHamiltonianDynamics):
 
     Example
     --------
-    >>> from pgmpy.models import JointGaussianDistribution
-    >>> from pgmpy.inference.continuous import LeapFrog
+    >>> from pgmpy.factors import JointGaussianDistribution
+    >>> from pgmpy.inference.continuous 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)
-    >>> new_pos, new_momentum, new_grad = LeapFrog(model, pos, momentum, 4.0).get_proposed_values()
+    >>> new_pos, new_momentum, new_grad = LeapFrog(model, pos, momentum, 4.0, GLPG).get_proposed_values()
     >>> new_pos
     array([ 70., -19.])
     >>> new_momentum
@@ -307,7 +306,7 @@ class LeapFrog(BaseSimulateHamiltonianDynamics):
     array([ 41., -58.])
     """
 
-    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, grad_log_position=None):
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf, grad_log_position=None):
 
         BaseSimulateHamiltonianDynamics.__init__(self, model, position, momentum,
                                                  stepsize, grad_log_pdf, grad_log_position)
@@ -324,7 +323,7 @@ class LeapFrog(BaseSimulateHamiltonianDynamics):
         # Take full step in time for updating position position
         position_bar = self.position + self.stepsize * momentum_bar
 
-        grad_log, _ = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+        grad_log, _ = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         # Take remaining half step in time for updating momentum
         momentum_bar = momentum_bar + 0.5 * self.stepsize * grad_log
@@ -361,7 +360,7 @@ class ModifiedEuler(BaseSimulateHamiltonianDynamics):
 
     Example
     --------
-    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.factors import JointGaussianDistribution
     >>> from pgmpy.inference.continuous import GradLogPDFGaussian, ModifiedEuler
     >>> import numpy as np
     >>> pos = np.array([2, 1])
@@ -379,7 +378,7 @@ class ModifiedEuler(BaseSimulateHamiltonianDynamics):
     array([-2.125, -1.1875])
     """
 
-    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, grad_log_position=None):
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf, grad_log_position=None):
 
         BaseSimulateHamiltonianDynamics.__init__(self, model, position, momentum,
                                                  stepsize, grad_log_pdf, grad_log_position)
@@ -396,6 +395,6 @@ class ModifiedEuler(BaseSimulateHamiltonianDynamics):
         # Take full step in time and update position
         position_bar = self.position + self.stepsize * momentum_bar
 
-        grad_log, _ = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+        grad_log, _ = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         return position_bar, momentum_bar, grad_log
diff --git a/pgmpy/inference/continuous/sampling.py b/pgmpy/inference/continuous/sampling.py
index 788acf6d..28945402 100644
--- a/pgmpy/inference/continuous/sampling.py
+++ b/pgmpy/inference/continuous/sampling.py
@@ -41,13 +41,13 @@ class HamiltonianMC(object):
 
     Example:
     --------
-    >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, LeapFrog
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
+    >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, LeapFrog, GradLogPDFGaussian
+    >>> from pgmpy.factors import JointGaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([-3, 4])
     >>> covariance = np.array([[3, 0.7], [0.7, 5]])
     >>> model = JGD(['x', 'y'], mean, covariance)
-    >>> sampler = HMC(model=model, grad_log_pdf=None, simulate_dynamics=LeapFrog)
+    >>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
     >>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
     ...                          trajectory_length=2, stepsize=0.4)
     >>> samples
@@ -74,9 +74,9 @@ class HamiltonianMC(object):
     CRC Press, 2011.
     """
 
-    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog):
+    def __init__(self, model, grad_log_pdf, simulate_dynamics=LeapFrog):
 
-        if grad_log_pdf is not None and not issubclass(grad_log_pdf, BaseGradLogPDF):
+        if not issubclass(grad_log_pdf, BaseGradLogPDF):
             raise TypeError("grad_log_pdf must be an instance of " +
                             "pgmpy.inference.base_continuous.BaseGradLogPDF")
 
@@ -96,14 +96,15 @@ class HamiltonianMC(object):
         """
 
         # Parameters to help in evaluating Joint distribution P(position, momentum)
-        _, logp = self.model.get_gradient_log_pdf(position, self.grad_log_pdf)
-        _, logp_bar = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+        _, logp = self.grad_log_pdf(position, self.model).get_gradient_log_pdf()
+        _, logp_bar = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         # acceptance_prob = P(position_bar, momentum_bar)/ P(position, momentum)
         potential_change = logp_bar - logp  # Negative change
         kinetic_change = 0.5 * np.float(np.dot(momentum_bar.T, momentum_bar) - np.dot(momentum.T, momentum))
 
-        return np.exp(potential_change - kinetic_change)  # acceptance probability
+        # acceptance probability
+        return np.exp(potential_change - kinetic_change)
 
     def _find_reasonable_stepsize(self, position, stepsize_app=1):
         """
@@ -158,7 +159,7 @@ class HamiltonianMC(object):
         if lsteps is None:
             lsteps = int(max(1, round(trajectory_length / stepsize, 0)))
 
-        grad_bar, _ = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+        grad_bar, _ = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         for _ in range(lsteps):
             position_bar, momentum_bar, grad_bar =\
@@ -212,7 +213,7 @@ class HamiltonianMC(object):
         --------
         >>> # Example if pandas is installed in working environment
         >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, GradLogPDFGaussian, ModifiedEuler
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, -1])
         >>> covariance = np.array([[1, 0.2], [0.2, 1]])
@@ -298,7 +299,7 @@ class HamiltonianMC(object):
         Examples
         --------
         >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, GradLogPDFGaussian as GLPG
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([4, -1])
         >>> covariance = np.array([[3, 0.4], [0.4, 3]])
@@ -370,7 +371,7 @@ class HamiltonianMCda(HamiltonianMC):
     Example:
     --------
     >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, LeapFrog
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
+    >>> from pgmpy.factors import JointGaussianDistribution 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]])
@@ -458,7 +459,7 @@ class HamiltonianMCda(HamiltonianMC):
         Examples
         ---------
         >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, 0.7], [0.7, 3]])
@@ -485,7 +486,8 @@ class HamiltonianMCda(HamiltonianMC):
             return HamiltonianMC.sample(self, initial_pos, num_samples, trajectory_length, stepsize)
 
         # stepsize is epsilon
-        mu = np.log(10.0 * stepsize)  # freely chosen point, after each iteration xt(/position) is shrunk towards it
+        # freely chosen point, after each iteration xt(/position) is shrunk towards it
+        mu = np.log(10.0 * stepsize)
         # log(10 * stepsize) large values to save computation
         # stepsize_bar is epsilon_bar
         stepsize_bar = 1.0
@@ -549,7 +551,7 @@ class HamiltonianMCda(HamiltonianMC):
         Examples
         --------
         >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, 1])
         >>> covariance = np.array([[1, 0.7], [0.7, 3]])
@@ -619,13 +621,13 @@ class NoUTurnSampler(HamiltonianMCda):
 
     Example:
     --------
-    >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, LeapFrog
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
+    >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, LeapFrog, GradLogPDFGaussian
+    >>> from pgmpy.factors import JointGaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([1, 2, 3])
     >>> covariance = np.array([[4, 0.1, 0.2], [0.1, 1, 0.3], [0.2, 0.3, 8]])
     >>> model = JGD(['x', 'y', 'z'], mean, covariance)
-    >>> sampler = NUTS(model=model, grad_log_pdf=None, simulate_dynamics=LeapFrog)
+    >>> sampler = NUTS(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
     >>> samples = sampler.sample(initial_pos=np.array([0.1, 0.9, 0.3]), num_samples=20000,stepsize=0.4)
     >>> samples
     rec.array([(0.1, 0.9, 0.3),
@@ -644,7 +646,7 @@ class NoUTurnSampler(HamiltonianMCda):
     Algorithm 3 : Efficient No-U-Turn Sampler
     """
 
-    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog):
+    def __init__(self, model, grad_log_pdf, simulate_dynamics=LeapFrog):
 
         super(NoUTurnSampler, self).__init__(model=model, grad_log_pdf=grad_log_pdf,
                                              simulate_dynamics=simulate_dynamics)
@@ -657,7 +659,7 @@ class NoUTurnSampler(HamiltonianMCda):
         position_bar, momentum_bar, _ = self.simulate_dynamics(self.model, position, momentum, stepsize,
                                                                self.grad_log_pdf).get_proposed_values()
 
-        _, logp_bar = self.model.get_gradient_log_pdf(position_bar, self.grad_log_pdf)
+        _, logp_bar = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         hamiltonian = logp_bar - 0.5 * np.dot(momentum_bar, momentum_bar)
 
@@ -737,7 +739,7 @@ class NoUTurnSampler(HamiltonianMCda):
         position_backward, position_forward = position, position
         momentum_backward, momentum_forward = momentum, momentum
         candidate_set_size = accept_set_bool = 1
-        _, log_pdf = self.model.get_gradient_log_pdf(position, self.grad_log_pdf)
+        _, log_pdf = self.grad_log_pdf(position, self.model).get_gradient_log_pdf()
 
         # Resample slice variable `u`
         slice_var = np.random.uniform(0, np.exp(log_pdf - 0.5 * np.dot(momentum, momentum)))
@@ -794,13 +796,13 @@ class NoUTurnSampler(HamiltonianMCda):
         Examples
         ---------
         >>> # If environment has a installation of pandas
-        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, GradLogPDFGaussian, LeapFrog
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([0, 0, 0])
         >>> covariance = np.array([[6, 0.7, 0.2], [0.7, 3, 0.9], [0.2, 0.9, 1]])
         >>> model = JGD(['x', 'y', 'z'], mean, covariance)
-        >>> sampler = NUTS(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
+        >>> sampler = NUTS(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
         >>> samples = sampler.sample(initial_pos=np.array([1, 1, 1]), num_samples=10, stepsize=0.4)
         >>> samples
                   x         y         z
@@ -860,13 +862,13 @@ class NoUTurnSampler(HamiltonianMCda):
 
         Examples
         ---------
-        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.inference.continuous import NoUTurnSampler as NUTS, GradLogPDFGaussian
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([11, -6])
         >>> covariance = np.array([[0.7, 0.2], [0.2, 14]])
         >>> model = JGD(['x', 'y'], mean, covariance)
-        >>> sampler = NUTS(model=model)
+        >>> sampler = NUTS(model=model, grad_log_pdf=GradLogPDFGaussian)
         >>> samples = sampler.generate_sample(initial_pos=np.array([1, 1]), num_samples=10, stepsize=0.4)
         >>> samples = np.array([sample for sample in samples])
         >>> samples
@@ -924,13 +926,13 @@ class NoUTurnSamplerDA(NoUTurnSampler):
 
     Example:
     --------
-    >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda
-    >>> from pgmpy.models import JointGaussianDistribution as JGD
+    >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
+    >>> from pgmpy.factors import JointGaussianDistribution as JGD
     >>> import numpy as np
     >>> mean = np.array([-1, 12, -3])
     >>> covariance = np.array([[-2, 7, 2], [7, 14, 4], [2, 4, -1]])
     >>> model = JGD(['x', 'v', 't'], mean, covariance)
-    >>> sampler = NUTSda(model=model)
+    >>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian)
     >>> samples = sampler.sample(initial_pos=np.array([0, 0, 0]), num_adapt=10, num_samples=10, stepsize=0.25)
     >>> samples
     rec.array([(0.0, 0.0, 0.0),
@@ -953,7 +955,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
     Algorithm 6 : No-U-Turn Sampler with Dual Averaging
     """
 
-    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog, delta=0.65):
+    def __init__(self, model, grad_log_pdf, simulate_dynamics=LeapFrog, delta=0.65):
 
         if not isinstance(delta, float) or delta > 1.0 or delta < 0.0:
             raise ValueError(
@@ -1023,7 +1025,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
         momentum_backward, momentum_forward = momentum, momentum
         candidate_set_size = accept_set_bool = 1
         position_m_1 = position
-        _, log_pdf = self.model.get_gradient_log_pdf(position, self.grad_log_pdf)
+        _, log_pdf = self.grad_log_pdf(position, self.model).get_gradient_log_pdf()
 
         # Resample slice variable `u`
         slice_var = np.random.uniform(0, np.exp(log_pdf - 0.5 * np.dot(momentum, momentum)))
@@ -1085,13 +1087,13 @@ class NoUTurnSamplerDA(NoUTurnSampler):
         Examples
         ---------
         >>> # If environment has a installation of pandas
-        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian as GLPG, LeapFrog
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian, LeapFrog
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([10, -13])
         >>> covariance = np.array([[16, -3], [-3, 13]])
         >>> model = JGD(['x', 'y'], mean, covariance)
-        >>> sampler = NUTSda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
+        >>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
         >>> samples = sampler.sample(initial_pos=np.array([12, -4]), num_adapt=10, num_samples=10, stepsize=0.1)
         >>> samples
                    x          y
@@ -1131,7 +1133,8 @@ class NoUTurnSamplerDA(NoUTurnSampler):
             samples[i] = position_m
 
             if i <= num_adapt:
-                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha, n_alpha)
+                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu,
+                                                                   i, alpha, n_alpha)
             else:
                 stepsize = stepsize_bar
 
@@ -1166,13 +1169,13 @@ class NoUTurnSamplerDA(NoUTurnSampler):
 
         Examples
         --------
-        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda
-        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> from pgmpy.inference.continuous import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
+        >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([1, -100])
         >>> covariance = np.array([[-12, 45], [45, -10]])
         >>> model = JGD(['a', 'b'], mean, covariance)
-        >>> sampler = NUTSda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
+        >>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
         >>> samples = sampler.generate_sample(initial_pos=np.array([12, -4]), num_adapt=10,
         ...                                   num_samples=10, stepsize=0.1)
         >>> samples
@@ -1214,7 +1217,8 @@ class NoUTurnSamplerDA(NoUTurnSampler):
             position_m, alpha, n_alpha = self._sample(position_m, stepsize)
 
             if i <= num_adapt:
-                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha, n_alpha)
+                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu,
+                                                                   i, alpha, n_alpha)
             else:
                 stepsize = stepsize_bar
 
diff --git a/pgmpy/models/JointGaussian.py b/pgmpy/models/JointGaussian.py
deleted file mode 100644
index c0fb9e79..00000000
--- a/pgmpy/models/JointGaussian.py
+++ /dev/null
@@ -1,83 +0,0 @@
-import numpy as np
-
-from pgmpy.utils import _check_1d_array_object, _check_length_equal
-from pgmpy.inference.continuous import BaseGradLogPDF
-
-
-class JointGaussianDistribution(object):
-    """
-    A naive gaussian container class.
-    Will be removed when Multivariate Distributions will be implemented
-
-    Paramters
-    ---------
-    variables: iterable of any hashable python object
-        The variables for which the distribution is defined.
-
-    mean: 1d array type of structure
-        Represents the mean of the distribution
-
-    covariance: A 2d array type object of size len(mean) x len(mean),
-        Covariance of the distribution.
-
-    """
-
-    def __init__(self, variables, mean, covariance):
-
-        mean = _check_1d_array_object(mean, 'mean')
-        _check_length_equal(mean, variables, 'mean', 'variables')
-
-        if not isinstance(covariance, (np.matrix, np.ndarray, list)):
-            raise TypeError(
-                "covariance must be a 2d array type object")
-        covariance = np.array(covariance)
-        if covariance.shape[0] != covariance.shape[1]:
-            raise ValueError(
-                "covariance must be a square in shape")
-
-        if mean.shape[0] != covariance.shape[0]:
-            raise ValueError("shape of mean vector should be d X 1 and" +
-                             " shape of covariance matrix should be d X d")
-        self.variables = variables
-        self.mean = mean
-        self.covariance = covariance
-        self.precision_matrix = np.linalg.inv(covariance)
-
-    def get_gradient_log_pdf(self, variable_assignment, grad_log_pdf=None):
-        """
-        A method for finding gradient log and log of distribution for given assignment
-
-        Parameters
-        ----------
-        variable_assignment: A 1d array like object
-            Assignment of variables at which gradient is to be computed
-
-        grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
-            A coustom class for finding gradient log and log for given assignment
-            If None, the will be computed
-
-        Example
-        ---------
-
-        Returns
-        --------
-        A tuple of following types (in order)
-
-        numpy.array: A 1d numpy.array representing value of gradient log of JointGaussianDistribution
-
-        float: A floating value representin log of JointGaussianDistribution
-        """
-        variable_assignment = _check_1d_array_object(variable_assignment, 'variable_assignment')
-        _check_length_equal(variable_assignment, self.variables, 'variables_assignment', 'variables')
-
-        if grad_log_pdf is not None:
-            if not issubclass(grad_log_pdf, BaseGradLogPDF):
-                raise TypeError("grad_log_pdf must be an instance" +
-                                " of pgmpy.inference.continuous.base.BaseGradLogPDF")
-            return grad_log_pdf(variable_assignment, self).get_gradient_log_pdf()
-
-        sub_vec = variable_assignment - self.mean
-        grad = - np.dot(self.precision_matrix, sub_vec)
-        log_pdf = 0.5 * np.dot(sub_vec, grad)
-
-        return grad, log_pdf
diff --git a/pgmpy/models/__init__.py b/pgmpy/models/__init__.py
index c0b5a484..f0ab5135 100644
--- a/pgmpy/models/__init__.py
+++ b/pgmpy/models/__init__.py
@@ -7,7 +7,6 @@ from .JunctionTree import JunctionTree
 from .DynamicBayesianNetwork import DynamicBayesianNetwork
 from .MarkovChain import MarkovChain
 from .NaiveBayes import NaiveBayes
-from .JointGaussian import JointGaussianDistribution
 
 __all__ = ['BayesianModel',
            'NoisyOrModel',
@@ -17,5 +16,4 @@ __all__ = ['BayesianModel',
            'ClusterGraph',
            'DynamicBayesianNetwork',
            'MarkovChain',
-           'NaiveBayes',
-           'JointGaussianDistribution']
+           'NaiveBayes']


Test patch:
diff --git a/pgmpy/tests/test_inference/test_base_continuous.py b/pgmpy/tests/test_inference/test_base_continuous.py
index 3eb56aa7..013244ec 100644
--- a/pgmpy/tests/test_inference/test_base_continuous.py
+++ b/pgmpy/tests/test_inference/test_base_continuous.py
@@ -2,7 +2,7 @@ import unittest
 
 import numpy as np
 
-from pgmpy.models import JointGaussianDistribution as JGD
+from pgmpy.factors import JointGaussianDistribution as JGD
 from pgmpy.inference.continuous import LeapFrog, ModifiedEuler, GradLogPDFGaussian
 
 
@@ -36,26 +36,32 @@ class TestLeapFrog(unittest.TestCase):
         momentum = [-1, -1, -1]
         self.test_with_grad_log = LeapFrog(model=self.test_model, position=position, momentum=momentum,
                                            stepsize=0.3, grad_log_pdf=GradLogPDFGaussian, grad_log_position=None)
-        grad_log_position, _ = self.test_model.get_gradient_log_pdf(variable_assignment=position)
+        grad_log_position, _ = GradLogPDFGaussian(position, self.test_model).get_gradient_log_pdf()
         self.test_without_grad_log = LeapFrog(model=self.test_model, position=position, momentum=momentum,
-                                              stepsize=0.4, grad_log_pdf=None, grad_log_position=grad_log_position)
+                                              stepsize=0.4, grad_log_pdf=GradLogPDFGaussian,
+                                              grad_log_position=grad_log_position)
 
     def test_errors(self):
         with self.assertRaises(TypeError):
-            LeapFrog(model=self.test_model, position=1, momentum=[1, 1], stepsize=0.1)
+            LeapFrog(model=self.test_model, position=1, momentum=[1, 1], stepsize=0.1,
+                     grad_log_pdf=GradLogPDFGaussian)
         with self.assertRaises(TypeError):
-            LeapFrog(model=self.test_model, position=[1, 1], momentum=1, stepsize=0.1)
+            LeapFrog(model=self.test_model, position=[1, 1], momentum=1, stepsize=0.1,
+                     grad_log_pdf=GradLogPDFGaussian)
         with self.assertRaises(ValueError):
-            LeapFrog(model=self.test_model, position=[1, 1], momentum=[1], stepsize=0.1)
+            LeapFrog(model=self.test_model, position=[1, 1], momentum=[1], stepsize=0.1,
+                     grad_log_pdf=GradLogPDFGaussian)
         with self.assertRaises(TypeError):
             LeapFrog(model=self.test_model, position=[1], momentum=[1], stepsize=0.1, grad_log_pdf=1)
         with self.assertRaises(ValueError):
-            LeapFrog(model=self.test_model, position=[1, 1], momentum=[1, 1], stepsize=0.1)
+            LeapFrog(model=self.test_model, position=[1, 1], momentum=[1, 1], stepsize=0.1,
+                     grad_log_pdf=GradLogPDFGaussian)
         with self.assertRaises(TypeError):
-            LeapFrog(model=self.test_model, position=[1, 1, 1], momentum=[1, 1, 1], stepsize=0.1, grad_log_position=1)
+            LeapFrog(model=self.test_model, position=[1, 1, 1], momentum=[1, 1, 1], stepsize=0.1,
+                     grad_log_pdf=GradLogPDFGaussian, grad_log_position=1)
         with self.assertRaises(ValueError):
             LeapFrog(model=self.test_model, position=[1, 1, 1], momentum=[1, 1, 1], stepsize=0.1,
-                     grad_log_position=[1, 1])
+                     grad_log_pdf=GradLogPDFGaussian, grad_log_position=[1, 1])
 
     def test_leapfrog_methods(self):
         new_pos, new_momentum, new_grad = self.test_with_grad_log.get_proposed_values()
@@ -83,9 +89,10 @@ class TestModifiedEuler(unittest.TestCase):
         momentum = [-2, 1]
         self.test_with_grad_log = ModifiedEuler(model=self.test_model, position=position, momentum=momentum,
                                                 stepsize=0.5, grad_log_pdf=GradLogPDFGaussian, grad_log_position=None)
-        grad_log_position, _ = self.test_model.get_gradient_log_pdf(variable_assignment=position)
+        grad_log_position, _ = GradLogPDFGaussian(position, self.test_model).get_gradient_log_pdf()
         self.test_without_grad_log = ModifiedEuler(model=self.test_model, position=position, momentum=momentum,
-                                                   stepsize=0.3, grad_log_pdf=None, grad_log_position=grad_log_position)
+                                                   stepsize=0.3, grad_log_pdf=GradLogPDFGaussian,
+                                                   grad_log_position=grad_log_position)
 
     def test_modified_euler_methods(self):
         new_pos, new_momentum, new_grad = self.test_with_grad_log.get_proposed_values()
diff --git a/pgmpy/tests/test_inference/test_continuous_sampling.py b/pgmpy/tests/test_inference/test_continuous_sampling.py
index 2d900f46..25239f6e 100644
--- a/pgmpy/tests/test_inference/test_continuous_sampling.py
+++ b/pgmpy/tests/test_inference/test_continuous_sampling.py
@@ -3,9 +3,9 @@ import unittest
 import numpy as np
 
 from pgmpy.inference.continuous import (HamiltonianMC as HMC, HamiltonianMCda as HMCda, ModifiedEuler,
-                                        GradLogPDFGaussian as GradGaussian, NoUTurnSampler as NUTS,
+                                        GradLogPDFGaussian, NoUTurnSampler as NUTS,
                                         NoUTurnSamplerDA as NUTSda)
-from pgmpy.models import JointGaussianDistribution as JGD
+from pgmpy.factors import JointGaussianDistribution as JGD
 
 
 class TestHMCInference(unittest.TestCase):
@@ -14,13 +14,13 @@ class TestHMCInference(unittest.TestCase):
         mean = [-1, 1, -1]
         covariance = np.array([[3, 0.8, 0.2], [0.8, 2, 0.3], [0.2, 0.3, 1]])
         self.test_model = JGD(['x', 'y', 'z'], mean, covariance)
-        self.hmc_sampler = HMCda(model=self.test_model)
+        self.hmc_sampler = HMCda(model=self.test_model, grad_log_pdf=GradLogPDFGaussian)
 
     def test_errors(self):
         with self.assertRaises(TypeError):
             HMCda(model=self.test_model, grad_log_pdf=1)
         with self.assertRaises(TypeError):
-            HMCda(model=self.test_model, grad_log_pdf=GradGaussian, simulate_dynamics=1)
+            HMCda(model=self.test_model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=1)
         with self.assertRaises(ValueError):
             HMCda(model=self.test_model, delta=-1)
         with self.assertRaises(TypeError):
@@ -31,15 +31,6 @@ class TestHMCInference(unittest.TestCase):
             HMC(model=self.test_model).sample(initial_pos=1, num_samples=1, trajectory_length=1)
         with self.assertRaises(TypeError):
             HMC(model=self.test_model).generate_sample(1, 1, 1).send(None)
-        # TODO: Remove them after implementation of JointGaussianDistribution
-        with self.assertRaises(TypeError):
-            JGD(['x', 'y'], mean=1, covariance=1)
-        with self.assertRaises(TypeError):
-            JGD(['x', 'y'], mean=np.ones(2), covariance=1)
-        with self.assertRaises(ValueError):
-            JGD(['x', 'y'], mean=np.ones(3), covariance=np.eye(2))
-        with self.assertRaises(ValueError):
-            JGD(['x', 'y'], mean=np.ones(3), covariance=[[1, 2, 3], [1, 2, 3]])
 
     def test_acceptance_prob(self):
         acceptance_probability = self.hmc_sampler._acceptance_prob(np.array([1, 2, 3]), np.array([2, 3, 4]),
@@ -100,33 +91,36 @@ class TestNUTSInference(unittest.TestCase):
         mean = np.array([-1, 1, 0])
         covariance = np.array([[6, 0.7, 0.2], [0.7, 3, 0.9], [0.2, 0.9, 1]])
         self.test_model = JGD(['x', 'y', 'z'], mean, covariance)
-        self.nuts_sampler = NUTSda(model=self.test_model, grad_log_pdf=GradGaussian)
+        self.nuts_sampler = NUTSda(model=self.test_model, grad_log_pdf=GradLogPDFGaussian)
 
     def test_errors(self):
         with self.assertRaises(TypeError):
             NUTS(model=self.test_model, grad_log_pdf=JGD)
         with self.assertRaises(TypeError):
-            NUTS(model=self.test_model, grad_log_pdf=None, simulate_dynamics=GradGaussian)
+            NUTS(model=self.test_model, grad_log_pdf=None, simulate_dynamics=GradLogPDFGaussian)
         with self.assertRaises(ValueError):
-            NUTSda(model=self.test_model, delta=-0.2)
+            NUTSda(model=self.test_model, delta=-0.2, grad_log_pdf=None)
         with self.assertRaises(ValueError):
-            NUTSda(model=self.test_model, delta=1.1)
+            NUTSda(model=self.test_model, delta=1.1, grad_log_pdf=GradLogPDFGaussian)
         with self.assertRaises(TypeError):
-            NUTS(self.test_model).sample(initial_pos={1, 1, 1}, num_samples=1)
+            NUTS(self.test_model, GradLogPDFGaussian).sample(initial_pos={1, 1, 1}, num_samples=1)
         with self.assertRaises(ValueError):
-            NUTS(self.test_model).sample(initial_pos=[1, 1], num_samples=1)
+            NUTS(self.test_model, GradLogPDFGaussian).sample(initial_pos=[1, 1], num_samples=1)
         with self.assertRaises(TypeError):
-            NUTSda(self.test_model).sample(initial_pos=1, num_samples=1, num_adapt=1)
+            NUTSda(self.test_model, GradLogPDFGaussian).sample(initial_pos=1, num_samples=1, num_adapt=1)
         with self.assertRaises(ValueError):
-            NUTSda(self.test_model).sample(initial_pos=[1, 1, 1, 1], num_samples=1, num_adapt=1)
+            NUTSda(self.test_model, GradLogPDFGaussian).sample(initial_pos=[1, 1, 1, 1], num_samples=1, num_adapt=1)
         with self.assertRaises(TypeError):
-            NUTS(self.test_model).generate_sample(initial_pos=0.1, num_samples=1).send(None)
+            NUTS(self.test_model, GradLogPDFGaussian).generate_sample(initial_pos=0.1, num_samples=1).send(None)
         with self.assertRaises(ValueError):
-            NUTS(self.test_model).generate_sample(initial_pos=(0, 1, 1, 1), num_samples=1).send(None)
+            NUTS(self.test_model, GradLogPDFGaussian).generate_sample(initial_pos=(0, 1, 1, 1),
+                                                                      num_samples=1).send(None)
         with self.assertRaises(TypeError):
-            NUTSda(self.test_model).generate_sample(initial_pos=[[1, 2, 3]], num_samples=1, num_adapt=1).send(None)
+            NUTSda(self.test_model, GradLogPDFGaussian).generate_sample(initial_pos=[[1, 2, 3]], num_samples=1,
+                                                                        num_adapt=1).send(None)
         with self.assertRaises(ValueError):
-            NUTSda(self.test_model).generate_sample(initial_pos=[1], num_samples=1, num_adapt=1).send(None)
+            NUTSda(self.test_model, GradLogPDFGaussian).generate_sample(initial_pos=[1], num_samples=1,
+                                                                        num_adapt=1).send(None)
 
     def test_sampling(self):
         np.random.seed(1010101)