Repository: pgmpy/pgmpy
Pull request #702 (2016-06-17 18:25:09)

Hamiltonian Monte Carlo & Hamiltonian Monte Carlo with dual averaging
@ankurankan I have send this PR to aid us in discussion.
 I was experimenting things with how to handle gradients (removing the gradient argument). I tested with two ways:
- First I tried handled `grad_log_pdf` argument on places itself depending upon how user passed the argument, if None was passed then I created a lambda function to call `model.get_gradient_log_pdf` otherwise I created a lambda function to use the custom class. But with this things were messy as I have to handle this parameter at two places, one in `sampling class` and other in `BaseSimulateHamiltonianDynamics` class.

- Second ( this PR implements it). Handle everything in `model.get_gradient_log_pdf`. This code is less messy, because every call is made to model.get_gradient_log_pdf and the method internally handles the rest so need of making suitable changes at different places.

How do you suggest I should handle the gradients ? You can look at the last commit to specifically see the changes I made https://github.com/pgmpy/pgmpy/pull/702/commits/748eb1fe13488bb8f0cf27a7064a67384ec3315e

After the discussion I'll close one of the PR. 



Patch:
diff --git a/pgmpy/factors/Factor.py b/pgmpy/factors/Factor.py
index 9fd5b696..6ef99206 100644
--- a/pgmpy/factors/Factor.py
+++ b/pgmpy/factors/Factor.py
@@ -96,7 +96,10 @@ class Factor(object):
         if isinstance(variables, six.string_types):
             raise TypeError("Variables: Expected type list or array like, got string")
 
-        values = np.array(values, dtype=float)
+        values = np.array(values)
+
+        if values.dtype != int and values.dtype != float:
+            raise TypeError("Values: Expected type int or type float, got ", values.dtype)
 
         if len(cardinality) != len(variables):
             raise ValueError("Number of elements in cardinality must be equal to number of variables")
diff --git a/pgmpy/inference/__init__.py b/pgmpy/inference/__init__.py
index f819359e..b9bc8ec8 100644
--- a/pgmpy/inference/__init__.py
+++ b/pgmpy/inference/__init__.py
@@ -4,6 +4,7 @@ from .ExactInference import BeliefPropagation
 from .Sampling import BayesianModelSampling, GibbsSampling
 from .dbn_inference import DBNInference
 from .mplp import Mplp
+from . import continuous
 
 __all__ = ['Inference',
            'VariableElimination',
@@ -11,4 +12,5 @@ __all__ = ['Inference',
            'BeliefPropagation',
            'BayesianModelSampling',
            'GibbsSampling',
-           'Mplp']
+           'Mplp',
+           'continuous']
diff --git a/pgmpy/inference/continuous/__init__.py b/pgmpy/inference/continuous/__init__.py
new file mode 100644
index 00000000..b5114e61
--- /dev/null
+++ b/pgmpy/inference/continuous/__init__.py
@@ -0,0 +1,12 @@
+from .base import (LeapFrog, ModifiedEuler, BaseSimulateHamiltonianDynamics,
+                   BaseGradLogPDF, GradLogPDFGaussian)
+from .sampling import HamiltonianMCda, HamiltonianMC
+
+
+__all__ = ['LeapFrog',
+           'ModifiedEuler',
+           'BaseSimulateHamiltonianDynamics',
+           'BaseGradLogPDF',
+           'GradLogPDFGaussian',
+           'HamiltonianMC',
+           'HamiltonianMCda']
diff --git a/pgmpy/inference/continuous/base.py b/pgmpy/inference/continuous/base.py
new file mode 100644
index 00000000..074dea8b
--- /dev/null
+++ b/pgmpy/inference/continuous/base.py
@@ -0,0 +1,401 @@
+import numpy as np
+
+from pgmpy.utils import _check_1d_array_object, _check_length_equal
+
+
+class BaseGradLogPDF(object):
+    """
+    Base class for evaluating gradient log of probability density function/ distribution
+
+    Classes inheriting this base class can be passed as an argument for
+    finding gradient log of probability distribution in inference algorithms
+
+    The class should initialize  self.grad_log and self.log_pdf
+
+    Parameters
+    ----------
+    variable_assignments : A 1d array like object (numpy.ndarray or list)
+        Vector representing values(assignments) of variables at which we want to find gradient and log
+
+    model : An instance of pgmpy.models
+
+    Examples
+    --------
+    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.inference.continuous import BaseGradLogPDF
+    >>> import numpy as np
+    >>> class GradLogGaussian(BaseGradLogPDF):
+    ...     def __init__(self, position, model):
+    ...         BaseGradLogPDF.__init__(self, position, model)
+    ...         self.grad_log, self.log_pdf = self._get_gradient_log_pdf()
+    ...     def _get_gradient_log_pdf(self):
+    ...         sub_vec = self.position - self.model.mean.flatten()
+    ...         grad = - np.dot(self.model.precision_matrix, sub_vec)
+    ...         log_pdf = 0.5 * np.float(np.dot(sub_vec, grad))
+    ...         return grad, log_pdf
+    >>> mean = np.array([1, 1])
+    >>> covariance = np.array([[1, 0.2], [0.2, 7]])
+    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> dist_param = np.array([0.1, 0.9])
+    >>> grad_logp, logp = GradLogGaussian(dist_param, model).get_gradient_log_pdf()
+    >>> logp
+    -0.4054597701149426
+    >>> grad_logp
+    array([ 0.90229885, -0.01149425])
+    """
+
+    def __init__(self, variable_assignments, model):
+
+        self.variable_assignments = _check_1d_array_object(variable_assignments, 'variable_assignments')
+        _check_length_equal(variable_assignments, model.variables, 'variable_assignments', 'model.variables')
+
+        self.model = model
+
+        # The gradient log of probability distribution at position
+        self.grad_log = None
+
+        # The gradient log of probability distribution at position
+        self.log_pdf = None
+
+    def get_gradient_log_pdf(self):
+        """
+        Returns the gradient log and log of model at given position
+
+        Returns
+        -------
+        Returns a tuple of following types
+
+        numpy.array: Representing gradient log of model at given position
+
+        float: Representing log of model at given position
+
+        Example
+        --------
+        >>> # Using implementation of GradLogPDFGaussian
+        >>> from pgmpy.inference.continuous import GradLogPDFGaussian
+        >>> from pgmpy.models import JointGaussianDistribution
+        >>> import numpy as np
+        >>> mean = np.array([1, 1])
+        >>> covariance = np.array([[1, -5], [-5, 2]])
+        >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+        >>> dist_param = np.array([0.6, 0.8])
+        >>> grad_logp, logp = GradLogPDFGaussian(dist_param, model).get_gradient_log_pdf()
+        >>> logp
+        0.025217391304347823
+        >>> grad_logp
+        array([-0.07826087, -0.09565217])
+        """
+        return self.grad_log, self.log_pdf
+
+
+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
+
+    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
+
+    Example
+    -------
+    >>> from pgmpy.inference.continuous import GradLogPDFGaussian
+    >>> from pgmpy.models import JointGaussianDistribution
+    >>> import numpy as np
+    >>> mean = np.array([3, 4])
+    >>> covariance = np.array([[5, 4], [4, 5]])
+    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> dist_param = np.array([12, 21])
+    >>> grad_logp, logp = GradLogPDFGaussian(dist_param, model).get_gradient_log_pdf()
+    >>> logp
+    -34.777777777777771
+    >>> grad_logp
+    array([ 2.55555556, -5.44444444])
+    """
+
+    def __init__(self, variable_assignments, model):
+        BaseGradLogPDF.__init__(self, variable_assignments, model)
+        self.grad_log, self.log_pdf = self._get_gradient_log_pdf()
+
+    def _get_gradient_log_pdf(self):
+        """
+        Method that finds gradient and its log at position
+        """
+        sub_vec = self.variable_assignments - self.model.mean.flatten()
+        grad = - np.dot(self.model.precision_matrix, sub_vec)
+        log_pdf = 0.5 * np.dot(sub_vec, grad)
+
+        return grad, log_pdf
+
+
+class BaseSimulateHamiltonianDynamics(object):
+    """
+    Base class for proposing new values of position and momentum by simulating Hamiltonian Dynamics.
+
+    Classes inheriting this base class can be passed as an argument for
+    simulate_dynamics in inference algorithms.
+
+    Parameters
+    ----------
+    model : An instance of pgmpy.models
+        Model for which DiscretizeTime object is initialized
+
+    position : A 1d array like object (numpy.ndarray or list)
+        Vector representing values of parameter position( or X)
+
+    momentum: A 1d array like object (numpy.ndarray or list)
+        Vector representing the proposed value for momentum (velocity)
+
+    stepsize: Float
+        stepsize for the simulating dynamics
+
+    grad_log_pdf : A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
+        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
+        If None, then will be calculated
+
+    Examples
+    --------
+    >>> from pgmpy.inference.continuous import BaseSimulateHamiltonianDynamics
+    >>> from pgmpy.models 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
+    >>> # self.new_grad_logp represents gradient log at new proposed value of position
+    >>> class ModifiedEuler(BaseSimulateHamiltonianDynamics):
+    ...     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)
+    ...         self.new_position, self.new_momentum, self.new_grad_logp = self._get_proposed_values()
+    ...     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)
+    ...         return position_bar, momentum_bar, grad_log_position
+    >>> pos = np.array([1, 2])
+    >>> momentum = np.array([0, 0])
+    >>> mean = np.array([0, 0])
+    >>> covariance = np.eye(2)
+    >>> model = JointGaussianDistribution(['x', 'y'], mean, covariance)
+    >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum,
+    ...                                                 0.25, GradLogPDFGaussian).get_proposed_values()
+    >>> new_pos
+    array([0.9375, 1.875])
+    >>> new_momentum
+    array([-0.25, -0.5])
+    >>> new_grad
+    array([-0.9375, -1.875])
+    """
+
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, 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):
+            raise TypeError("grad_log_pdf must be an instance" +
+                            " of pgmpy.inference.continuous.base.BaseGradLogPDF")
+
+        _check_length_equal(position, momentum, 'position', 'momentum')
+        _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)
+
+        else:
+            grad_log_positon = _check_1d_array_object(grad_log_position, 'grad_log_position')
+            _check_length_equal(grad_log_position, position, 'grad_log_position', 'position')
+
+        self.position = position
+        self.momentum = momentum
+        self.stepsize = stepsize
+        self.model = model
+        self.grad_log_pdf = grad_log_pdf
+        self.grad_log_position = grad_log_position
+
+        # new_position is the new proposed position, new_momentum is the new proposed momentum, new_grad_lop
+        # is the value of grad log at new_position
+        self.new_position = self.new_momentum = self.new_grad_logp = None
+
+    def get_proposed_values(self):
+        """
+        Returns the new proposed values of position and momentum
+
+        Returns
+        -------
+        Returns a tuple of following type (in order)
+
+        numpy.array: A 1d numpy.array representing new proposed value of position
+
+        numpy.array: A 1d numpy.array representing new proposed value of momentum
+
+        numpy.array: A 1d numpy.array representing gradient of log distribution at new proposed value of position
+
+        Example
+        -------
+        >>> # Using implementation of ModifiedEuler
+        >>> from pgmpy.inference.continuous import ModifiedEuler
+        >>> from pgmpy.models 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
+        array([ 3.04666667,  4.21      ])
+        >>> new_momentum
+        array([ 0.06666667,  0.3       ])
+        >>> new_grad
+        array([-1.34888889, -1.07      ])
+        """
+        return self.new_position, self.new_momentum, self.new_grad_logp
+
+
+class LeapFrog(BaseSimulateHamiltonianDynamics):
+    """
+    Class for simulating hamiltonian dynamics using leapfrog method
+    Inherits pgmpy.inference.base_continuous.BaseSimulateHamiltonianDynamics
+
+    Parameters
+    ----------
+    model : An instance of pgmpy.models
+        Model for which DiscretizeTime object is initialized
+
+    position : A 1d array like object (numpy.ndarray or list)
+        Vector representing values of parameter position( or X)
+
+    momentum: A 1d array like object (numpy.ndarray or list)
+        Vector representing the proposed value for momentum (velocity)
+
+    stepsize: Float
+        stepsize for the simulating dynamics
+
+    grad_log_pdf : A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
+        A class for evaluating gradient log and log of distribution for a given assignment of variables
+        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
+        If None, then will be calculated
+
+    Example
+    --------
+    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.inference.continuous import LeapFrog
+    >>> 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
+    array([ 70., -19.])
+    >>> new_momentum
+    array([  99., -121.])
+    >>> new_grad
+    array([ 41., -58.])
+    """
+
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, grad_log_position=None):
+
+        BaseSimulateHamiltonianDynamics.__init__(self, model, position, momentum,
+                                                 stepsize, grad_log_pdf, grad_log_position)
+
+        self.new_position, self.new_momentum, self.new_grad_logp = self._get_proposed_values()
+
+    def _get_proposed_values(self):
+        """
+        Method to perform time splitting using leapfrog
+        """
+        # Take half step in time for updating momentum
+        momentum_bar = self.momentum + 0.5 * self.stepsize * self.grad_log_position
+
+        # 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)
+
+        # Take remaining half step in time for updating momentum
+        momentum_bar = momentum_bar + 0.5 * self.stepsize * grad_log
+
+        return position_bar, momentum_bar, grad_log
+
+
+class ModifiedEuler(BaseSimulateHamiltonianDynamics):
+    """
+    Class for simulating Hamiltonian Dynamics using Modified euler method
+    Inherits pgmpy.inference.base_continuous.BaseSimulateHamiltonianDynamics
+
+    Parameters
+    ----------
+    model: An instance of pgmpy.models
+        Model for which DiscretizeTime object is initialized
+
+    position: A 1d array like object (numpy.ndarray or list)
+        Vector representing values of parameter position( or X)
+
+    momentum: A 1d array like object (numpy.ndarray or list)
+        Vector representing the proposed value for momentum (velocity)
+
+    stepsize: Float
+        stepsize for the simulating dynamics
+
+    grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
+        A class for finding gradient log and log of distribution
+        If None, then will use model.get_gradient_log_pdf
+
+    grad_log_position: A 1d array like object, defaults to None
+        Vector representing gradient log at given position
+        If None, then will be calculated
+
+    Example
+    --------
+    >>> from pgmpy.models import JointGaussianDistribution
+    >>> from pgmpy.inference.continuous 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)
+    >>> new_pos, new_momentum, new_grad = ModifiedEuler(model, pos, momentum,
+    ...                                                 0.25, GradLogPDFGaussian).get_proposed_values()
+    >>> new_pos
+    array([2.125, 1.1875])
+    >>> new_momentum
+    array([0.5, 0.75])
+    >>> new_grad
+    array([-2.125, -1.1875])
+    """
+
+    def __init__(self, model, position, momentum, stepsize, grad_log_pdf=None, grad_log_position=None):
+
+        BaseSimulateHamiltonianDynamics.__init__(self, model, position, momentum,
+                                                 stepsize, grad_log_pdf, grad_log_position)
+
+        self.new_position, self.new_momentum, self.new_grad_logp = self._get_proposed_values()
+
+    def _get_proposed_values(self):
+        """
+        Method to perform time splitting using Modified euler method
+        """
+        # Take full step in time and update momentum
+        momentum_bar = self.momentum + self.stepsize * self.grad_log_position
+
+        # 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)
+
+        return position_bar, momentum_bar, grad_log
diff --git a/pgmpy/inference/continuous/sampling.py b/pgmpy/inference/continuous/sampling.py
new file mode 100644
index 00000000..9dcbcc37
--- /dev/null
+++ b/pgmpy/inference/continuous/sampling.py
@@ -0,0 +1,593 @@
+"""
+    A collection of methods for sampling from continuous models in pgmpy
+"""
+from math import sqrt
+
+import numpy as np
+try:
+    import pandas as pd
+    HAS_PANDAS = True
+except ImportError:
+    HAS_PANDAS = False
+
+from pgmpy.inference.continuous import LeapFrog, BaseGradLogPDF, BaseSimulateHamiltonianDynamics
+from pgmpy.utils import _check_1d_array_object, _check_length_equal
+
+
+class HamiltonianMC(object):
+    """
+    Class for performing sampling using simple
+    Hamiltonian Monte Carlo
+
+    Parameters:
+    -----------
+    model: An instance pgmpy.models
+        Model from which sampling has to be done
+
+    grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
+        A class to find log and gradient of log distribution for a given assignment
+        If None, then will use model.get_gradient_log_pdf
+
+    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
+        A class to propose future values of momentum and position in time by simulating
+        Hamiltonian Dynamics
+
+    Public Methods:
+    ---------------
+    sample()
+    generate_sample()
+
+    Example:
+    --------
+    >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, LeapFrog
+    >>> from pgmpy.models 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)
+    >>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
+    ...                          trajectory_length=2, stepsize=0.4)
+    >>> samples
+    array([(5e-324, 5e-324), (-2.2348941964735225, 4.43066330647519),
+           (-2.316454719617516, 7.430291195678112), ...,
+           (-1.1443831048872348, 3.573135519428842),
+           (-0.2325915892988598, 4.155961788010201),
+           (-0.7582492446601238, 3.5416519297297056)], 
+          dtype=[('x', '<f8'), ('y', '<f8')])
+
+    >>> samples = np.array([samples[var_name] for var_name in model.variables])
+    >>> np.cov(samples)
+    array([[ 3.0352818 ,  0.71379304],
+           [ 0.71379304,  4.91776713]])
+    >>> sampler.accepted_proposals
+    9932.0
+    >>> sampler.acceptance_rate
+    0.9932
+
+    References
+    ----------
+    R.Neal. Handbook of Markov Chain Monte Carlo,
+    chapter 5: MCMC Using Hamiltonian Dynamics.
+    CRC Press, 2011.
+    """
+
+    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog):
+
+        if grad_log_pdf is not None and not issubclass(grad_log_pdf, BaseGradLogPDF):
+            raise TypeError("grad_log_pdf must be an instance of " +
+                            "pgmpy.inference.base_continuous.BaseGradLogPDF")
+
+        if not issubclass(simulate_dynamics, BaseSimulateHamiltonianDynamics):
+            raise TypeError("split_time must be an instance of " +
+                            "pgmpy.inference.base_continuous.BaseSimulateHamiltonianDynamics")
+
+        self.model = model
+        self.grad_log_pdf = grad_log_pdf
+        self.simulate_dynamics = simulate_dynamics
+        self.accepted_proposals = 0.0
+        self.acceptance_rate = 0
+
+    def _acceptance_prob(self, position, position_bar, momentum, momentum_bar):
+        """
+        Returns the acceptance probability for given new position(position) and momentum
+        """
+
+        # 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)
+
+        # 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
+
+    def _find_reasonable_stepsize(self, position, stepsize_app=1):
+        """
+        Method for choosing initial value of stepsize
+
+        References
+        -----------
+        Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively
+        Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
+        Machine Learning Research 15 (2014) 1351-1381
+        Algorithm 4 : Heuristic for choosing an initial value of epsilon
+        """
+        # momentum = N(0, I)
+        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)
+
+        # Take a single step in time
+        position_bar, momentum_bar, _ =\
+            self.simulate_dynamics(self.model, position, momentum,
+                                   stepsize_app, self.grad_log_pdf).get_proposed_values()
+
+        acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)
+
+        # a = 2I[acceptance_prob] -1
+        a = 2 * (acceptance_prob > 0.5) - 1
+
+        condition = (acceptance_prob ** a) > (2 ** (-a))
+
+        while condition:
+            stepsize_app = (2 ** a) * stepsize_app
+
+            position_bar, momentum_bar, _ =\
+                self.simulate_dynamics(self.model, position, momentum,
+                                       stepsize_app, self.grad_log_pdf).get_proposed_values()
+
+            acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)
+
+            condition = (acceptance_prob ** a) > (2 ** (-a))
+
+        return stepsize_app
+
+    def _sample(self, position, trajectory_length, stepsize, lsteps=None):
+        """
+        Runs a single sampling iteration to return a sample
+        """
+        # Resampling momentum
+        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)
+
+        # position_m here will be the previous sampled value of position
+        position_bar, momentum_bar = position.copy(), momentum
+
+        # Number of steps L to simulate dynamics
+        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)
+
+        for _ in range(lsteps):
+            position_bar, momentum_bar, grad_bar =\
+                self.simulate_dynamics(self.model, position_bar, momentum_bar,
+                                       stepsize, self.grad_log_pdf, grad_bar).get_proposed_values()
+
+        acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)
+
+        # Metropolis acceptance probability
+        alpha = min(1, acceptance_prob)
+
+        # Accept or reject the new proposed value of position, i.e position_bar
+        if np.random.rand() < alpha:
+            position = position_bar.copy()
+            self.accepted_proposals += 1.0
+
+        return position, alpha
+
+    def sample(self, initial_pos, num_samples, trajectory_length, stepsize=None):
+        """
+        Method to return samples using Hamiltonian Monte Carlo
+
+        Parameters
+        ----------
+        initial_pos: A 1d array like object
+            Vector representing values of parameter position, the starting
+            state in markov chain.
+
+        num_samples: int
+            Number of samples to be generated
+
+        trajectory_length: int or float
+            Target trajectory length, stepsize * number of steps(L),
+            where L is the number of steps taken per HMC iteration,
+            and stepsize is step size for splitting time method.
+
+        stepsize: float , defaults to None
+            The stepsize for proposing new values of position and momentum in simulate_dynamics
+            If None, then will be choosen suitably
+
+
+        Returns
+        -------
+        Returns two different types (based on installations)
+
+        pandas.DataFrame: Returns samples as pandas.DataFrame if environment has a installation of pandas
+
+        numpy.recarray: Returns samples in form of numpy recorded arrays (numpy.recarray)
+
+        Examples
+        --------
+        >>> # 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
+        >>> import numpy as np
+        >>> mean = np.array([1, -1])
+        >>> covariance = np.array([[1, 0.2], [0.2, 1]])
+        >>> model = JGD(['x', 'y'], mean, covariance)
+        >>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=ModifiedEuler)
+        >>> samples = sampler.sample(np.array([1, 1]), num_samples = 5,
+        ...                          trajectory_length=6, stepsize=0.25)
+        >>> samples
+                       x              y
+        0   1.000000e+00   1.000000e+00
+        1   1.592133e+00   1.152911e+00
+        2   1.608700e+00   1.315349e+00
+        3   1.608700e+00   1.315349e+00
+        4   6.843856e-01   6.237043e-01
+        >>> mean = np.array([4, 1, -1])
+        >>> covariance = np.array([[1, 0.7 , 0.8], [0.7, 1, 0.2], [0.8, 0.2, 1]])
+        >>> model = JGD(['x', 'y', 'z'], mean, covariance)
+        >>> sampler = HMC(model=model, grad_log_pdf=GLPG)
+        >>> samples = sampler.sample(np.array([1, 1]), num_samples = 10000,
+        ...                          trajectory_length=6, stepsize=0.25)
+        >>> np.cov(samples.values.T)
+        array([[ 1.00795398,  0.71384233,  0.79802097],
+               [ 0.71384233,  1.00633524,  0.21313767],
+               [ 0.79802097,  0.21313767,  0.98519017]])
+        """
+
+        self.accepted_proposals = 1.0
+        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
+        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
+
+        if stepsize is None:
+            stepsize = self._find_reasonable_stepsize(initial_pos)
+
+        types = [(var_name, 'float') for var_name in self.model.variables]
+        samples = np.zeros(num_samples, dtype=types).view(np.recarray)
+
+        # Assigning after converting into tuple because value was being changed after assignment
+        # Reason for this is unknown
+        samples[0] = tuple(initial_pos)
+        position_m = initial_pos
+
+        lsteps = int(max(1, round(trajectory_length / stepsize, 0)))
+        for i in range(1, num_samples):
+
+            # Genrating sample
+            position_m, _ = self._sample(position_m, trajectory_length, stepsize, lsteps)
+            samples[i] = position_m
+
+        self.acceptance_rate = self.accepted_proposals / num_samples
+
+        if HAS_PANDAS == True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
+
+    def generate_sample(self, initial_pos, num_samples, trajectory_length, stepsize=None):
+        """
+        Method returns a generator type object whose each iteration yields a sample
+        using Hamiltonian Monte Carlo
+
+        Parameters
+        ----------
+        initial_pos: A 1d array like object
+            Vector representing values of parameter position, the starting
+            state in markov chain.
+
+        num_samples: int
+            Number of samples to be generated
+
+        trajectory_length: int or float
+            Target trajectory length, stepsize * number of steps(L),
+            where L is the number of steps taken per HMC iteration,
+            and stepsize is step size for splitting time method.
+
+        stepsize: float , defaults to None
+            The stepsize for proposing new values of position and momentum in simulate_dynamics
+            If None, then will be choosen suitably
+
+        Returns
+        -------
+        genrator: yielding a 1d numpy.array type object for a sample
+
+        Examples
+        --------
+        >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, GradLogPDFGaussian as GLPG
+        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> import numpy as np
+        >>> mean = np.array([4, -1])
+        >>> covariance = np.array([[3, 0.4], [0.4, 3]])
+        >>> model = JGD(['x', 'y'], mean, covariance)
+        >>> sampler = HMC(model=model, grad_log_pdf=GLPG)
+        >>> gen_samples = sampler.generate_sample(np.array([-1, 1]), num_samples = 10000,
+        ...                                       trajectory_length=2, stepsize=0.25)
+        >>> samples_array = np.array([sample for sample in gen_samples])
+        >>> samples_array
+        array([[ 0.1467264 ,  0.27143857],
+               [ 4.0371448 ,  0.15871274],
+               [ 3.24656208, -1.03742621],
+               ..., 
+               [ 6.45975905,  1.97941306],
+               [ 4.89007171,  0.15413156],
+               [ 5.9528083 ,  1.92983158]])
+        >>> np.cov(samples_array.T)
+        array([[ 2.95692642,  0.4379419 ],
+               [ 0.4379419 ,  3.00939434]])
+        >>> sampler.acceptance_rate
+        0.9969
+        """
+
+        self.accepted_proposals = 0
+        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
+        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
+
+        if stepsize is None:
+            stepsize = self._find_reasonable_stepsize(initial_pos)
+
+        lsteps = int(max(1, round(trajectory_length / stepsize, 0)))
+        position_m = initial_pos.copy()
+
+        for i in range(0, num_samples):
+
+            position_m, _ = self._sample(position_m, trajectory_length, stepsize, lsteps)
+
+            yield position_m
+
+        self.acceptance_rate = self.accepted_proposals / num_samples
+
+
+class HamiltonianMCda(HamiltonianMC):
+    """
+    Class for performing sampling in Continuous model
+    using Hamiltonian Monte Carlo with dual averaging for
+    adaptaion of parameter stepsize.
+
+    Parameters:
+    -----------
+    model: An instance pgmpy.models
+        Model from which sampling has to be done
+
+    grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
+        Class to compute the log and gradient log of distribution
+
+    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
+        Class to propose future states of position and momentum in time by simulating
+        HamiltonianDynamics
+
+    delta: float (in between 0 and 1), defaults to 0.65
+        The target HMC acceptance probability
+
+    Public Methods:
+    ---------------
+    sample()
+    generate_sample()
+
+    Example:
+    --------
+    >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, LeapFrog
+    >>> from pgmpy.models 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]])
+    >>> model = JGD(['x', 'y', 'z'], mean, covariance)
+    >>> sampler = HMCda(model=model)
+    >>> samples = sampler.sample(np.array([0, 0, 0]), num_adapt=10000, num_samples = 10000, trajectory_length=7)
+    >>> samples_array = np.array([samples[var_name] for var_name in model.variables])
+    >>> np.cov(samples_array)
+    array([[ 1.83023816,  0.40449162,  0.51200707],
+           [ 0.40449162,  2.85863596,  0.76747343],
+           [ 0.51200707,  0.76747343,  3.87020982]])
+    >>> sampler.acceptance_rate
+    0.9929
+
+    References
+    -----------
+    Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively
+    Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
+    Machine Learning Research 15 (2014) 1351-1381
+    Algorithm 5 : Hamiltonian Monte Carlo with dual averaging
+    """
+
+    def __init__(self, model, grad_log_pdf=None, simulate_dynamics=LeapFrog, delta=0.65):
+
+        if not isinstance(delta, float) or delta > 1.0 or delta < 0.0:
+            raise AttributeError(
+                "delta should be a floating value in between 0 and 1")
+
+        self.delta = delta
+
+        super(HamiltonianMCda, self).__init__(model=model, grad_log_pdf=grad_log_pdf,
+                                              simulate_dynamics=simulate_dynamics)
+
+    def _adapt_params(self, stepsize, stepsize_bar, h_bar, mu, index_i, alpha):
+        """
+        Run tha adaptation for stepsize for better proposals of position
+        """
+        gamma = 0.05  # free parameter that controls the amount of shrinkage towards mu
+        t0 = 10.0  # free parameter that stabilizes the initial iterations of the algorithm
+        kappa = 0.75
+        # See equation (6) section 3.2.1 for details
+
+        estimate = 1.0 / (index_i + t0)
+        h_bar = (1 - estimate) * h_bar + estimate * (self.delta - alpha)
+
+        stepsize = np.exp(mu - sqrt(index_i) / gamma * h_bar)
+        i_kappa = index_i ** (-kappa)
+        stepsize_bar = np.exp(i_kappa * np.log(stepsize) + (1 - i_kappa) * np.log(stepsize_bar))
+
+        return stepsize, stepsize_bar, h_bar
+
+    def sample(self, initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None):
+        """
+        Method to return samples using Hamiltonian Monte Carlo
+
+        Parameters
+        ----------
+        initial_pos: A 1d array like object
+            Vector representing values of parameter position, the starting
+            state in markov chain.
+
+        num_adapt: int
+            The number of interations to run the adaptation of stepsize
+
+        num_samples: int
+            Number of samples to be generated
+
+        trajectory_length: int or float
+            Target trajectory length, stepsize * number of steps(L),
+            where L is the number of steps taken per HMC iteration,
+            and stepsize is step size for splitting time method.
+
+        stepsize: float , defaults to None
+            The stepsize for proposing new values of position and momentum in simulate_dynamics
+            If None, then will be choosen suitably
+
+        Returns
+        -------
+        Returns two different types (based on installations)
+
+        pandas.DataFrame: Returns samples as pandas.DataFrame if environment has a installation of pandas
+
+        numpy.recarray: Returns samples in form of numpy recorded arrays (numpy.recarray)
+
+        Examples
+        ---------
+        >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
+        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> import numpy as np
+        >>> mean = np.array([1, 1])
+        >>> covariance = np.array([[1, 0.7], [0.7, 3]])
+        >>> model = JGD(['x', 'y'], mean, covariance)
+        >>> sampler = HMCda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
+        >>> samples = sampler.sample(np.array([1, 1]), num_adapt=10000,
+        ...                          num_samples = 10000, trajectory_length=2, stepsize=None)
+        >>> samples_array = np.array([samples[var_name] for var_name in model.variables])
+        >>> np.cov(samples_array)
+        array([[ 0.98432155,  0.66517394],
+               [ 0.66517394,  2.95449533]])
+
+        """
+
+        self.accepted_proposals = 1.0
+
+        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
+        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
+
+        if stepsize is None:
+            stepsize = self._find_reasonable_stepsize(initial_pos)
+
+        if num_adapt <= 1:  # Return samples genrated using Simple HMC algorithm
+            return HamiltonianMC.sample(self, initial_pos, num_samples, stepsize)
+
+        # stepsize is epsilon
+        mu = np.log(10.0 * stepsize)  # freely chosen point, after each iteration xt(/position) is shrunk towards it
+        # log(10 * stepsize) large values to save computation
+        # stepsize_bar is epsilon_bar
+        stepsize_bar = 1.0
+        h_bar = 0.0
+        # See equation (6) section 3.2.1 for details
+
+        types = [(var_name, 'float') for var_name in self.model.variables]
+        samples = np.zeros(num_samples, dtype=types).view(np.recarray)
+        samples[0] = tuple(initial_pos)
+        position_m = initial_pos
+
+        for i in range(1, num_samples):
+
+            # Genrating sample
+            position_m, alpha = self._sample(position_m, trajectory_length, stepsize)
+            samples[i] = position_m
+
+            # Adaptation of stepsize till num_adapt iterations
+            if i <= num_adapt:
+                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha)
+            else:
+                stepsize = stepsize_bar
+
+        self.acceptance_rate = self.accepted_proposals / num_samples
+
+        if HAS_PANDAS == True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
+
+    def generate_sample(self, initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None):
+        """
+        Method returns a generator type object whose each iteration yields a sample
+        using Hamiltonian Monte Carlo
+
+        Parameters
+        ----------
+        initial_pos: A 1d array like object
+            Vector representing values of parameter position, the starting
+            state in markov chain.
+
+        num_adapt: int
+            The number of interations to run the adaptation of stepsize
+
+        num_samples: int
+            Number of samples to be generated
+
+        trajectory_length: int or float
+            Target trajectory length, stepsize * number of steps(L),
+            where L is the number of steps taken to propose new values of position and momentum
+            per HMC iteration and stepsize is step size.
+
+        stepsize: float , defaults to None
+            The stepsize for proposing new values of position and momentum in simulate_dynamics
+            If None, then will be choosen suitably
+
+        Returns
+        -------
+        genrator: yielding a numpy.array type object for a sample
+
+        Examples
+        --------
+        >>> from pgmpy.inference.continuous import HamiltonianMCda as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
+        >>> from pgmpy.models import JointGaussianDistribution as JGD
+        >>> import numpy as np
+        >>> mean = np.array([1, 1])
+        >>> covariance = np.array([[1, 0.7], [0.7, 3]])
+        >>> model = JGD(['x', 'y'], mean, covariance)
+        >>> sampler = HMCda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
+        >>> gen_samples = sampler.generate_sample(np.array([1, 1]), num_adapt=10000,
+        ...                                       num_samples = 10000, trajectory_length=2, stepsize=None)
+        >>> samples_array = np.array([sample for sample in gen_samples])
+        >>> np.cov(samples_array.T)
+        array([[ 0.98432155,  0.69517394],
+               [ 0.69517394,  2.95449533]])
+        """
+        self.accepted_proposals = 0
+        initial_pos = _check_1d_array_object(initial_pos, 'initial_pos')
+        _check_length_equal(initial_pos, self.model.variables, 'initial_pos', 'model.variables')
+
+        if stepsize is None:
+            stepsize = self._find_reasonable_stepsize(initial_pos)
+
+        if num_adapt <= 1:  # return sample generated using Simple HMC algorithm
+            for sample in HamiltonianMC.generate_sample(self, initial_pos, num_samples, stepsize, trajectory_length):
+                yield sample
+            return
+        mu = np.log(10.0 * stepsize)
+
+        stepsize_bar = 1.0
+        h_bar = 0.0
+
+        position_m = initial_pos.copy()
+        num_adapt += 1
+
+        for i in range(1, num_samples + 1):
+
+            position_m, alpha = self._sample(position_m, trajectory_length, stepsize)
+
+            if i <= num_adapt:
+                stepsize, stepsize_bar, h_bar = self._adapt_params(stepsize, stepsize_bar, h_bar, mu, i, alpha)
+            else:
+                stepsize = stepsize_bar
+
+            yield position_m
+
+        self.acceptance_rate = self.accepted_proposals / num_samples
diff --git a/pgmpy/models/JointGaussian.py b/pgmpy/models/JointGaussian.py
new file mode 100644
index 00000000..c0fb9e79
--- /dev/null
+++ b/pgmpy/models/JointGaussian.py
@@ -0,0 +1,83 @@
+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 f0ab5135..c0b5a484 100644
--- a/pgmpy/models/__init__.py
+++ b/pgmpy/models/__init__.py
@@ -7,6 +7,7 @@ from .JunctionTree import JunctionTree
 from .DynamicBayesianNetwork import DynamicBayesianNetwork
 from .MarkovChain import MarkovChain
 from .NaiveBayes import NaiveBayes
+from .JointGaussian import JointGaussianDistribution
 
 __all__ = ['BayesianModel',
            'NoisyOrModel',
@@ -16,4 +17,5 @@ __all__ = ['BayesianModel',
            'ClusterGraph',
            'DynamicBayesianNetwork',
            'MarkovChain',
-           'NaiveBayes']
+           'NaiveBayes',
+           'JointGaussianDistribution']
diff --git a/pgmpy/utils/__init__.py b/pgmpy/utils/__init__.py
index ac5faf30..d96f2d4b 100644
--- a/pgmpy/utils/__init__.py
+++ b/pgmpy/utils/__init__.py
@@ -1,7 +1,10 @@
 from .mathext import cartesian, sample_discrete
 from .state_name import StateNameInit, StateNameDecorator
+from .check_functions import _check_1d_array_object, _check_length_equal
 
 __all__ = ['cartesian',
            'sample_discrete',
            'StateNameInit',
-           'StateNameDecorator']
+           'StateNameDecorator',
+           '_check_1d_array_object',
+           '_check_length_equal']
diff --git a/pgmpy/utils/check_functions.py b/pgmpy/utils/check_functions.py
new file mode 100644
index 00000000..78e4bfca
--- /dev/null
+++ b/pgmpy/utils/check_functions.py
@@ -0,0 +1,24 @@
+"""
+Contains simple check functions
+"""
+import numpy as np
+
+def _check_1d_array_object(parameter, name_param):
+    """
+    Checks whether given parameter is a 1d array like object, and returns a numpy array object
+    """
+    if isinstance(parameter, (np.ndarray, list, tuple, np.matrix)):
+        parameter = np.array(parameter)
+        if parameter.ndim != 1:
+            raise TypeError("{} should be a 1d array type object".format(name_param))
+    else:
+        raise TypeError("{} should be a 1d array type object".format(name_param))
+
+    return parameter
+
+def _check_length_equal(param_1, param_2, name_param_1, name_param_2):
+    """
+    Raises an error when the length of given two arguments is not equal
+    """
+    if len(param_1) != len(param_2):
+        raise ValueError("Length of {} must be same as Length of {}".format(name_param_1, name_param_2))


Test patch:
diff --git a/pgmpy/tests/test_factors/test_Factor.py b/pgmpy/tests/test_factors/test_Factor.py
index fe13cfb2..3f6a66f8 100644
--- a/pgmpy/tests/test_factors/test_Factor.py
+++ b/pgmpy/tests/test_factors/test_Factor.py
@@ -33,6 +33,11 @@ class TestFactorInit(unittest.TestCase):
         self.assertRaises(ValueError, Factor, ['x1', 'x2', 'x3'], [2, 2, 2], np.ones(9))
 
     def test_class_init_typeerror(self):
+        self.assertRaises(TypeError, Factor, ['x1', 'x2', 'x3'], [2, 1, 1], ['val1', 'val2'])
+        self.assertRaises(TypeError, Factor, ['x1', 'x2', 'x3'], [2, 1, 1], [1, 'val1'])
+        self.assertRaises(TypeError, Factor, ['x1', 'x2', 'x3'], [2, 1, 1], ['val1', 1])
+        self.assertRaises(TypeError, Factor, ['x1', 'x2', 'x3'], [2, 1, 1], [0.1, 'val1'])
+        self.assertRaises(TypeError, Factor, ['x1', 'x2', 'x3'], [2, 1, 1], ['val1', 0.1])
         self.assertRaises(TypeError, Factor, 'x1', [3], [1, 2, 3])
         self.assertRaises(ValueError, Factor, ['x1', 'x1', 'x3'], [2, 3, 2], range(12))
 
diff --git a/pgmpy/tests/test_inference/test_base_continuous.py b/pgmpy/tests/test_inference/test_base_continuous.py
new file mode 100644
index 00000000..3eb56aa7
--- /dev/null
+++ b/pgmpy/tests/test_inference/test_base_continuous.py
@@ -0,0 +1,103 @@
+import unittest
+
+import numpy as np
+
+from pgmpy.models import JointGaussianDistribution as JGD
+from pgmpy.inference.continuous import LeapFrog, ModifiedEuler, GradLogPDFGaussian
+
+
+class TestGradLogPDFGaussian(unittest.TestCase):
+
+    def setUp(self):
+        mean = np.array([1, 2, 3, 4])
+        covariance = np.array([[1, 0.2, 0.4, 0.7], [0.2, 2, 0.5, 0.8], [0.4, 0.5, 3, 0.6], [0.7, 0.8, 0.6, 4]])
+        self.test_model = JGD(['x', 'y', 'z', 't'], mean, covariance)
+        self.test_gradient = GradLogPDFGaussian([0, 0, 0, 0], self.test_model)
+
+    def test_error(self):
+        with self.assertRaises(TypeError):
+            GradLogPDFGaussian(1, self.test_model)
+        with self.assertRaises(ValueError):
+            GradLogPDFGaussian([1, 1], self.test_model)
+
+    def test_gradient(self):
+        grad, log = self.test_gradient.get_gradient_log_pdf()
+        np.testing.assert_almost_equal(grad, np.array([0.05436475, 0.49454937, 0.75465073, 0.77837868]))
+        np.testing.assert_almost_equal(log, -3.21046521505)
+
+
+class TestLeapFrog(unittest.TestCase):
+
+    def setUp(self):
+        mean = np.array([-1, 1, -1])
+        covariance = np.array([[1, 0.6, 0.5], [0.6, 2, 0.3], [0.5, 0.3, 1]])
+        self.test_model = JGD(['x', 'y', 'z'], mean, covariance)
+        position = [0, 0, 0]
+        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)
+        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)
+
+    def test_errors(self):
+        with self.assertRaises(TypeError):
+            LeapFrog(model=self.test_model, position=1, momentum=[1, 1], stepsize=0.1)
+        with self.assertRaises(TypeError):
+            LeapFrog(model=self.test_model, position=[1, 1], momentum=1, stepsize=0.1)
+        with self.assertRaises(ValueError):
+            LeapFrog(model=self.test_model, position=[1, 1], momentum=[1], stepsize=0.1)
+        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)
+        with self.assertRaises(TypeError):
+            LeapFrog(model=self.test_model, position=[1, 1, 1], momentum=[1, 1, 1], stepsize=0.1, 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])
+
+    def test_leapfrog_methods(self):
+        new_pos, new_momentum, new_grad = self.test_with_grad_log.get_proposed_values()
+        np.testing.assert_almost_equal(new_pos, np.array([-0.35634146, -0.25609756, -0.33]))
+        np.testing.assert_almost_equal(new_momentum, np.array([-1.3396624, -0.70344884, -1.16963415]))
+        np.testing.assert_almost_equal(new_grad, np.array([-1.0123835, 1.00139798, -0.46422764]))
+        new_pos, new_momentum, new_grad = self.test_without_grad_log.get_proposed_values()
+        np.testing.assert_almost_equal(new_pos, np.array([-0.5001626, -0.32195122, -0.45333333]))
+        np.testing.assert_almost_equal(new_momentum, np.array([-1.42947981, -0.60709102, -1.21246612]))
+        np.testing.assert_almost_equal(new_grad, np.array([-0.89536651, 0.98893516, -0.39566396]))
+
+    def tearDown(self):
+        del self.test_model
+        del self.test_with_grad_log
+        del self.test_without_grad_log
+
+
+class TestModifiedEuler(unittest.TestCase):
+
+    def setUp(self):
+        mean = np.array([0, 0])
+        covariance = np.array([[-1, 0.8], [0.8, 3]])
+        self.test_model = JGD(['x', 'y'], mean, covariance)
+        position = [0, 0]
+        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)
+        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)
+
+    def test_modified_euler_methods(self):
+        new_pos, new_momentum, new_grad = self.test_with_grad_log.get_proposed_values()
+        np.testing.assert_almost_equal(new_pos, np.array([-1.0, 0.5]))
+        np.testing.assert_almost_equal(new_momentum, np.array([-2.0, 1.0]))
+        np.testing.assert_almost_equal(new_grad, np.array([-0.93406593, 0.08241758]))
+        new_pos, new_momentum, new_grad = self.test_without_grad_log.get_proposed_values()
+        np.testing.assert_almost_equal(new_pos, np.array([-0.6, 0.3]))
+        np.testing.assert_almost_equal(new_momentum, np.array([-2.0, 1.0]))
+        np.testing.assert_almost_equal(new_grad, np.array([-0.56043956, 0.04945055]))
+
+    def tearDown(self):
+        del self.test_model
+        del self.test_with_grad_log
+        del self.test_without_grad_log
diff --git a/pgmpy/tests/test_inference/test_continuous_sampling.py b/pgmpy/tests/test_inference/test_continuous_sampling.py
new file mode 100644
index 00000000..8951f5a5
--- /dev/null
+++ b/pgmpy/tests/test_inference/test_continuous_sampling.py
@@ -0,0 +1,91 @@
+import unittest
+
+import numpy as np
+
+from pgmpy.inference.continuous import (HamiltonianMC as HMC, HamiltonianMCda as HMCda, ModifiedEuler,
+                                        GradLogPDFGaussian as GradGaussian)
+from pgmpy.models import JointGaussianDistribution as JGD
+
+
+class TestHMCInference(unittest.TestCase):
+
+    def setUp(self):
+        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)
+
+    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)
+        with self.assertRaises(TypeError):
+            self.hmc_sampler.sample(initial_pos=1, num_adapt=1, num_samples=1, trajectory_length=1)
+        with self.assertRaises(TypeError):
+            self.hmc_sampler.generate_sample(1, 1, 1, 1).send(None)
+        with self.assertRaises(TypeError):
+            HMC(model=self.test_model).sample(position0=1, num_samples=1, trajectory_length=1)
+        with self.assertRaises(TypeError):
+            HMC(model=self.test_model).generate_sample(1, 1, 1).send(None)
+        # TODO: Remove them after implementation of JointGaussianDistribution
+        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]),
+                                                                   np.array([1, -1, 1]), np.array([0, 0, 0]))
+        np.testing.assert_almost_equal(acceptance_probability, 0.0347363)
+
+    def test_find_resonable_stepsize(self):
+        np.random.seed(987654321)
+        stepsize = self.hmc_sampler._find_reasonable_stepsize(np.array([-1, 1, -1]))
+        np.testing.assert_almost_equal(stepsize, 2.0)
+
+    def test_adapt_params(self):
+        stepsize, stepsize_bar, h_bar = self.hmc_sampler._adapt_params(0.0025, 1, 1, np.log(0.025), 2, 1)
+        np.testing.assert_almost_equal(stepsize, 3.13439452e-13)
+        np.testing.assert_almost_equal(stepsize_bar, 3.6742481e-08)
+        np.testing.assert_almost_equal(h_bar, 0.8875)
+
+    def test_sample(self):
+        # Seeding is done for _find_reasonable_stepsize method
+        # Testing sample method simple HMC
+        np.random.seed(3124141)
+        samples = self.hmc_sampler.sample(initial_pos=[0.3, 0.4, 0.2], num_adapt=0,
+                                          num_samples=10000, trajectory_length=4)
+        covariance = np.cov(samples.values.T)
+        self.assertTrue(np.linalg.norm(covariance - self.test_model.covariance) < 3)
+
+        # Testing sample of method of HMCda
+        np.random.seed(3124141)
+        samples = self.hmc_sampler.sample(initial_pos=[0.6, 0.2, 0.8], num_adapt=10000,
+                                          num_samples=10000, trajectory_length=4)
+        covariance = np.cov(samples.values.T)
+        self.assertTrue(np.linalg.norm(covariance - self.test_model.covariance) < 0.3)
+
+        # Testing generate_sample method of simple HMC
+        np.random.seed(3124141)
+        gen_samples = self.hmc_sampler.generate_sample(initial_pos=[0.3, 0.4, 0.2], num_adapt=0,
+                                                       num_samples=10000, trajectory_length=4)
+        samples = np.array([sample for sample in gen_samples])
+        covariance = np.cov(samples.T)
+        self.assertTrue(np.linalg.norm(covariance - self.test_model.covariance) < 3)
+
+        # Testing sample of method of HMCda
+        np.random.seed(3124141)
+        gen_samples = self.hmc_sampler.generate_sample(initial_pos=[0.6, 0.2, 0.8], num_adapt=10000,
+                                                       num_samples=10000, trajectory_length=4)
+        samples = np.array([sample for sample in gen_samples])
+        covariance = np.cov(samples.T)
+        self.assertTrue(np.linalg.norm(covariance - self.test_model.covariance) < 0.3)
+
+    def tearDown(self):
+        del self.hmc_sampler
+        del self.test_model