Repository: pgmpy/pgmpy
Pull request #796 (2016-12-22 20:10:39)

Fixes Sampling Markov Chain PR#791 #789
 This PR solves solves 2 Issues #791 Sampling problem and #789 Stationarity Test for Distribution. In Stationarity Test we check against steady state values obtained from the Eigen value and check the steady state distribution obtained from the samples that we generate these values have to be within a user defined parameter called tolerance to constitute a Passing of the stationarity test.


Patch:
diff --git a/.travis.yml b/.travis.yml
index f1972807..3dae1288 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,8 +2,8 @@ language: python
 virtualenv:
     system_site_packages: false
 env:
-    - DISTRIB="conda" PYTHON_VERSION="3.5" COVERAGE="true"
     - DISTRIB="conda" PYTHON_VERSION="3.4" COVERAGE="true"
+    - DISTRIB="conda" PYTHON_VERSION="3.3" COVERAGE="true"
     - DISTRIB="conda" PYTHON_VERSION="2.7" COVERAGE="true"
 
 before_install:
diff --git a/appveyor.yml b/appveyor.yml
index afe897ba..3e08acec 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -21,7 +21,7 @@ build_script:
   - SET PATH=C:\Py;C:\Py\Scripts;C:\Py\Library\bin;%PATH%
   - conda config --set always_yes yes
   - conda update conda --quiet
-  - conda install numpy scipy pandas pip nose mock wrapt --quiet
+  - conda install numpy scipy pandas pip nose mock --quiet
   - pip install pyparsing
   - python setup.py develop
   - set "GIT_DIR=%cd%"
diff --git a/pgmpy/base/OptionalDependency.py b/pgmpy/base/OptionalDependency.py
deleted file mode 100644
index 96729e91..00000000
--- a/pgmpy/base/OptionalDependency.py
+++ /dev/null
@@ -1,6 +0,0 @@
-#This module initializes flags for optional dependencies
-try:
-    import pandas
-    HAS_PANDAS = True
-except ImportError:
-    HAS_PANDAS = False
diff --git a/pgmpy/base/__init__.py b/pgmpy/base/__init__.py
index 29433287..d2246743 100644
--- a/pgmpy/base/__init__.py
+++ b/pgmpy/base/__init__.py
@@ -1,7 +1,5 @@
 from .UndirectedGraph import UndirectedGraph
 from .DirectedGraph import DirectedGraph
-from .OptionalDependency import HAS_PANDAS
 
 __all__ = ['UndirectedGraph',
-           'DirectedGraph',
-           'HAS_PANDAS']
+           'DirectedGraph']
diff --git a/pgmpy/factors/continuous/CanonicalFactor.py b/pgmpy/factors/continuous/CanonicalFactor.py
index 993696ec..2c64f4ae 100644
--- a/pgmpy/factors/continuous/CanonicalFactor.py
+++ b/pgmpy/factors/continuous/CanonicalFactor.py
@@ -260,7 +260,9 @@ class CanonicalFactor(ContinuousFactor):
         h_j = self.h[index_to_reduce]
 
         # The values for the reduced variables.
-        y = np.array([value for var, value in values]).reshape(len(index_to_reduce), 1)
+        y = np.array([value for index,
+                      value in sorted([(self.variables.index(var), value) for var,
+                                      value in values])]).reshape(len(index_to_reduce), 1)
 
         phi.variables = [self.variables[index] for index in index_to_keep]
         phi.K = K_i_i
diff --git a/pgmpy/factors/continuous/JointGaussianDistribution.py b/pgmpy/factors/continuous/JointGaussianDistribution.py
index 4aaabc7c..6a093742 100644
--- a/pgmpy/factors/continuous/JointGaussianDistribution.py
+++ b/pgmpy/factors/continuous/JointGaussianDistribution.py
@@ -230,8 +230,9 @@ class JointGaussianDistribution(ContinuousFactor):
 
         mu_j = self.mean[index_to_keep]
         mu_i = self.mean[index_to_reduce]
-        x_i = np.array([value for var, value in values]).reshape(len(index_to_reduce), 1)
-
+        x_i = np.array([value for index,
+                        value in sorted([(self.variables.index(var), value) for var,
+                                         value in values])]).reshape(len(index_to_reduce), 1)
         sig_i_j = self.covariance[np.ix_(index_to_reduce, index_to_keep)]
         sig_j_i = self.covariance[np.ix_(index_to_keep, index_to_reduce)]
         sig_i_i_inv = np.linalg.inv(self.covariance[np.ix_(index_to_reduce, index_to_reduce)])
@@ -330,7 +331,7 @@ class JointGaussianDistribution(ContinuousFactor):
         >>> phi.g
         -6.51533
         """
-        from pgmpy.factors.continuous import CanonicalFactor
+        from pgmpy.factors import CanonicalFactor
 
         mu = self.mean
         sigma = self.covariance
diff --git a/pgmpy/models/BayesianModel.py b/pgmpy/models/BayesianModel.py
index aa0f6501..6e6923ef 100644
--- a/pgmpy/models/BayesianModel.py
+++ b/pgmpy/models/BayesianModel.py
@@ -456,7 +456,7 @@ class BayesianModel(DirectedGraph):
 
     def local_independencies(self, variables):
         """
-        Returns an instance of Independencies containing the local independencies
+        Returns a independencies object containing the local independencies
         of each of the variables.
 
         Parameters
@@ -471,8 +471,12 @@ class BayesianModel(DirectedGraph):
         >>> student.add_edges_from([('diff', 'grade'), ('intel', 'grade'),
         >>>                         ('grade', 'letter'), ('intel', 'SAT')])
         >>> ind = student.local_independencies('grade')
-        >>> ind
-        (grade _|_ SAT | diff, intel)
+        >>> ind.event1
+        {'grade'}
+        >>> ind.event2
+        {'SAT'}
+        >>> ind.event3
+        {'diff', 'intel'}
         """
         def dfs(node):
             """
diff --git a/pgmpy/models/MarkovChain.py b/pgmpy/models/MarkovChain.py
index 846ff993..9392f699 100644
--- a/pgmpy/models/MarkovChain.py
+++ b/pgmpy/models/MarkovChain.py
@@ -4,6 +4,7 @@ from warnings import warn
 
 import numpy as np
 from pandas import DataFrame
+from scipy.linalg import eig
 
 from pgmpy.factors.discrete import State
 from pgmpy.utils import sample_discrete
@@ -282,11 +283,11 @@ class MarkovChain(object):
             for st in self.transition_models[var]:
                 var_states[var][st] = list(self.transition_models[var][st].keys())
                 var_values[var][st] = list(self.transition_models[var][st].values())
-                samples[var][st] = sample_discrete(var_states[var][st], var_values[var][st])[0]
+                samples[var][st] = sample_discrete(var_states[var][st], var_values[var][st], size=size)
 
         for i in range(size - 1):
             for j, (var, st) in enumerate(self.state):
-                next_st = samples[var][st]
+                next_st = samples[var][st][i]
                 self.state[j] = State(var, next_st)
             sampled.loc[i + 1] = [st for var, st in self.state]
 
@@ -305,7 +306,7 @@ class MarkovChain(object):
         >>> from pgmpy.models.MarkovChain import MarkovChain as MC
         >>> from pgmpy.factors.discrete import State
         >>> model = MC(['intel', 'diff'], [3, 2])
-        >>> intel_tm = {0: {0: 0.2, 1: 0.4, 2:0.4}, 1: {0: 0, 1: 0.5, 2: 0.5}, 2: {2: 0.1, 1:0.9}}
+        >>> intel_tm = {0: {0: 0.2, 1: 0.4, 2:0.4}, 1: {0: 0, 1: 0.5, 2: 0.5}, 2: {2: 0.5, 1:0.5}}
         >>> model.add_transition_model('intel', intel_tm)
         >>> diff_tm = {0: {0: 0.5, 1: 0.5}, 1: {0: 0.25, 1:0.75}}
         >>> model.add_transition_model('diff', diff_tm)
@@ -367,6 +368,58 @@ class MarkovChain(object):
                 self.state[j] = State(var, next_st)
             yield self.state[:]
 
+    def is_stationarity(self, tolerance=0.2, sample=None):
+        """
+        Checks if the given markov chain is stationary and checks the steady state
+        probablity values for the state are consistent.
+
+        Parameters:
+        -----------
+        tolerance: float
+            represents the diff between actual steady state value and the computed value
+        sample: [State(i,j)]
+            represents the list of state which the markov chain has sampled
+
+        Return Type:
+        ------------
+        Boolean
+        True, if the markov chain converges to steady state distribution within the tolerance
+        False, if the markov chain does not converge to steady state distribution within tolerance
+
+        Examples:
+        ---------
+        >>> from pgmpy.models.MarkovChain import MarkovChain
+        >>> from pgmpy.factors.discrete import State
+        >>> model = MarkovChain()
+        >>> model.add_variables_from(['intel', 'diff'], [3, 2])
+        >>> intel_tm = {0: {0: 0.2, 1: 0.4, 2:0.4}, 1: {0: 0, 1: 0.5, 2: 0.5}, 2: {0: 0.3, 1: 0.3, 2: 0.4}}
+        >>> model.add_transition_model('intel', intel_tm)
+        >>> diff_tm = {0: {0: 0.5, 1: 0.5}, 1: {0: 0.25, 1:0.75}}
+        >>> model.add_transition_model('diff', diff_tm)
+        >>> model.is_stationarity()
+        True
+        """
+        keys = self.transition_models.keys()
+        return_val = True
+        for k in keys:
+            # convert dict to numpy matrix
+            transition_mat = np.array([np.array(list(self.transition_models[k][i].values()))
+                                       for i in self.transition_models[k].keys()], dtype=np.float)
+            S, U = eig(transition_mat.T)
+            stationary = np.array(U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)
+            stationary = (stationary / np.sum(stationary)).real
+
+            probabilites = []
+            window_size = 10000 if sample is None else len(sample)
+            for i in range(0, transition_mat.shape[0]):
+                probabilites.extend(self.prob_from_sample([State(k, i)], window_size=window_size))
+            if any(np.abs(i) > tolerance for i in np.subtract(probabilites, stationary)):
+                return_val = return_val and False
+            else:
+                return_val = return_val and True
+
+        return return_val
+
     def random_state(self):
         """
         Generates a random state of the Markov Chain.
diff --git a/pgmpy/models/NaiveBayes.py b/pgmpy/models/NaiveBayes.py
index 8f7ab219..bcc67e18 100644
--- a/pgmpy/models/NaiveBayes.py
+++ b/pgmpy/models/NaiveBayes.py
@@ -141,8 +141,9 @@ class NaiveBayes(BayesianModel):
 
     def local_independencies(self, variables):
         """
-        Returns an instance of Independencies containing the local independencies
-        of each of the variables.
+        Returns a list of independencies objects containing the local independencies
+        of each of the variables. If local independencies does not exist for a variable
+        it gives a None for that variable.
 
 
         Parameters
@@ -157,12 +158,15 @@ class NaiveBayes(BayesianModel):
         >>> model.add_edges_from([('a', 'b'), ('a', 'c'), ('a', 'd')])
         >>> ind = model.local_independencies('b')
         >>> ind
-        (b _|_ d, c | a)
+        [(b _|_ d, c | a)]
         """
-        independencies = Independencies()
+        independencies = []
         for variable in [variables] if isinstance(variables, str) else variables:
             if variable != self.parent_node:
-                independencies.add_assertions([variable, list(set(self.children_nodes) - set(variable)), self.parent_node])
+                independencies.append(Independencies(
+                    [variable, list(set(self.children_nodes) - set(variable)), self.parent_node]))
+            else:
+                independencies.append(None)
         return independencies
 
     def fit(self, data, parent_node=None, estimator_type=None):
diff --git a/pgmpy/sampling/HMC.py b/pgmpy/sampling/HMC.py
index 0b48c87f..c0f6493b 100644
--- a/pgmpy/sampling/HMC.py
+++ b/pgmpy/sampling/HMC.py
@@ -6,9 +6,15 @@ from __future__ import division
 from math import sqrt
 
 import numpy as np
+try:
+    import pandas as pd
+    HAS_PANDAS = True
+except ImportError:
+    HAS_PANDAS = False
 
 from pgmpy.utils import _check_1d_array_object, _check_length_equal
-from pgmpy.sampling import LeapFrog, BaseSimulateHamiltonianDynamics, BaseGradLogPDF, _return_samples
+from pgmpy.sampling import LeapFrog, BaseSimulateHamiltonianDynamics, BaseGradLogPDF
+
 
 
 class HamiltonianMC(object):
@@ -44,7 +50,7 @@ class HamiltonianMC(object):
     >>> model = JGD(['x', 'y'], mean, covariance)
     >>> 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, return_type='recarray')
+    ...                          trajectory_length=2, stepsize=0.4)
     >>> samples
     rec.array([(1.0, 1.0), (-3.1861687131079086, 3.7940994520145654),
      (-1.6920542547310844, 6.347410703806017), ...,
@@ -173,7 +179,7 @@ class HamiltonianMC(object):
 
         return position, alpha
 
-    def sample(self, initial_pos, num_samples, trajectory_length, stepsize=None, return_type='dataframe'):
+    def sample(self, initial_pos, num_samples, trajectory_length, stepsize=None):
         """
         Method to return samples using Hamiltonian Monte Carlo
 
@@ -195,16 +201,18 @@ class HamiltonianMC(object):
             The stepsize for proposing new values of position and momentum in simulate_dynamics
             If None, then will be choosen suitably
 
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
 
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        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.sampling import HamiltonianMC as HMC, GradLogPDFGaussian, ModifiedEuler
         >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
         >>> import numpy as np
@@ -213,7 +221,7 @@ class HamiltonianMC(object):
         >>> 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, return_type='dataframe')
+        ...                          trajectory_length=6, stepsize=0.25)
         >>> samples
                        x              y
         0   1.000000e+00   1.000000e+00
@@ -226,7 +234,7 @@ class HamiltonianMC(object):
         >>> 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, return_type='dataframe')
+        ...                          trajectory_length=6, stepsize=0.25)
         >>> np.cov(samples.values.T)
         array([[ 1.00795398,  0.71384233,  0.79802097],
                [ 0.71384233,  1.00633524,  0.21313767],
@@ -257,7 +265,10 @@ class HamiltonianMC(object):
 
         self.acceptance_rate = self.accepted_proposals / num_samples
 
-        return _return_samples(return_type, samples)
+        if HAS_PANDAS is True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
 
     def generate_sample(self, initial_pos, num_samples, trajectory_length, stepsize=None):
         """
@@ -288,7 +299,7 @@ class HamiltonianMC(object):
 
         Examples
         --------
-        >>> from pgmpy.sampling import HamiltonianMC as HMC, GradLogPDFGaussian as GLPG
+        >>> from pgmpy.inference.continuous import HamiltonianMC as HMC, GradLogPDFGaussian as GLPG
         >>> from pgmpy.factors import JointGaussianDistribution as JGD
         >>> import numpy as np
         >>> mean = np.array([4, -1])
@@ -360,15 +371,14 @@ class HamiltonianMCDA(HamiltonianMC):
 
     Example:
     --------
-    >>> from pgmpy.sampling import HamiltonianMCDA as HMCda, LeapFrog, GradLogPDFGaussian as GLPG
+    >>> from pgmpy.sampling import HamiltonianMCDA as HMCda, LeapFrog
     >>> from pgmpy.factors.continuous 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, grad_log_pdf=GLPG)
-    >>> samples = sampler.sample(np.array([0, 0, 0]), num_adapt=10000, num_samples = 10000, trajectory_length=7,
-    ...                          return_type='recarray')
+    >>> 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],
@@ -414,7 +424,7 @@ class HamiltonianMCDA(HamiltonianMC):
 
         return stepsize, stepsize_bar, h_bar
 
-    def sample(self, initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None, return_type='dataframe'):
+    def sample(self, initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None):
         """
         Method to return samples using Hamiltonian Monte Carlo
 
@@ -439,13 +449,13 @@ class HamiltonianMCDA(HamiltonianMC):
             The stepsize for proposing new values of position and momentum in simulate_dynamics
             If None, then will be choosen suitably
 
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
-
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        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
         ---------
@@ -456,8 +466,8 @@ class HamiltonianMCDA(HamiltonianMC):
         >>> 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, return_type='recarray')
+        >>> 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],
@@ -504,7 +514,10 @@ class HamiltonianMCDA(HamiltonianMC):
 
         self.acceptance_rate = self.accepted_proposals / num_samples
 
-        return _return_samples(return_type, samples)
+        if HAS_PANDAS is True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
 
     def generate_sample(self, initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None):
         """
diff --git a/pgmpy/sampling/NUTS.py b/pgmpy/sampling/NUTS.py
index c941c944..2efcb22e 100644
--- a/pgmpy/sampling/NUTS.py
+++ b/pgmpy/sampling/NUTS.py
@@ -3,8 +3,14 @@
 from __future__ import division
 
 import numpy as np
-
-from pgmpy.sampling import HamiltonianMCDA, LeapFrog, _return_samples
+try:
+    import pandas as pd
+    HAS_PANDAS = True
+except ImportError:
+    HAS_PANDAS = False
+
+from pgmpy.sampling import HamiltonianMCDA
+from pgmpy.sampling import LeapFrog
 from pgmpy.utils import _check_1d_array_object, _check_length_equal
 
 
@@ -39,8 +45,7 @@ class NoUTurnSampler(HamiltonianMCDA):
     >>> 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=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
-    >>> samples = sampler.sample(initial_pos=np.array([0.1, 0.9, 0.3]), num_samples=20000,
-    ...                          stepsize=0.4, return_type='recarray')
+    >>> samples = sampler.sample(initial_pos=np.array([0.1, 0.9, 0.3]), num_samples=20000,stepsize=0.4)
     >>> samples
     rec.array([(0.1, 0.9, 0.3),
      (-0.27303886844752756, 0.5028580705249155, 0.2895768065049909),
@@ -180,7 +185,7 @@ class NoUTurnSampler(HamiltonianMCDA):
 
         return position
 
-    def sample(self, initial_pos, num_samples, stepsize=None, return_type='dataframe'):
+    def sample(self, initial_pos, num_samples, stepsize=None):
         """
         Method to return samples using No U Turn Sampler
 
@@ -197,26 +202,25 @@ class NoUTurnSampler(HamiltonianMCDA):
             The stepsize for proposing new values of position and momentum in simulate_dynamics
             If None, then will be choosen suitably
 
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
-
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        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.sampling import NoUTurnSampler as NUTS, GradLogPDFGaussian, LeapFrog
-        >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
+        >>> # If environment has a installation of pandas
+        >>> 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=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
-        >>> samples = sampler.sample(initial_pos=np.array([1, 1, 1]), num_samples=10,
-        ...                          stepsize=0.4, return_type='dataframe')
+        >>> samples = sampler.sample(initial_pos=np.array([1, 1, 1]), num_samples=10, stepsize=0.4)
         >>> samples
                   x         y         z
         0  1.000000  1.000000  1.000000
@@ -247,7 +251,10 @@ class NoUTurnSampler(HamiltonianMCDA):
             position_m = self._sample(position_m, stepsize)
             samples[i] = position_m
 
-        return _return_samples(return_type, samples)
+        if HAS_PANDAS is True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
 
     def generate_sample(self, initial_pos, num_samples, stepsize=None):
         """
@@ -343,8 +350,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
     >>> covariance = np.array([[-2, 7, 2], [7, 14, 4], [2, 4, -1]])
     >>> model = JGD(['x', 'v', 't'], mean, covariance)
     >>> 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, return_type='recarray')
+    >>> samples = sampler.sample(initial_pos=np.array([0, 0, 0]), num_adapt=10, num_samples=10, stepsize=0.25)
     >>> samples
     rec.array([(0.0, 0.0, 0.0),
      (0.06100992691638076, -0.17118088764170125, 0.14048470935160887),
@@ -467,7 +473,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
 
         return position, alpha, n_alpha
 
-    def sample(self, initial_pos, num_adapt, num_samples, stepsize=None, return_type='dataframe'):
+    def sample(self, initial_pos, num_adapt, num_samples, stepsize=None):
         """
         Returns samples using No U Turn Sampler with dual averaging
 
@@ -487,17 +493,17 @@ class NoUTurnSamplerDA(NoUTurnSampler):
             The stepsize for proposing new values of position and momentum in simulate_dynamics
             If None, then will be choosen suitably
 
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
-
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        Returns two different types (based on installations)
 
+        pandas.DataFrame: Returns samples as pandas.DataFrame if environment has a installation of pandas
+
+        numpy.recarray: Returns samples in form of numpy recorded arrays (numpy.recarray)
 
         Examples
         ---------
+        >>> # If environment has a installation of pandas
         >>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian, LeapFrog
         >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
         >>> import numpy as np
@@ -505,8 +511,7 @@ class NoUTurnSamplerDA(NoUTurnSampler):
         >>> covariance = np.array([[16, -3], [-3, 13]])
         >>> model = JGD(['x', 'y'], mean, covariance)
         >>> 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, return_type='dataframe')
+        >>> samples = sampler.sample(initial_pos=np.array([12, -4]), num_adapt=10, num_samples=10, stepsize=0.1)
         >>> samples
                    x          y
         0  12.000000  -4.000000
@@ -550,7 +555,10 @@ class NoUTurnSamplerDA(NoUTurnSampler):
             else:
                 stepsize = stepsize_bar
 
-        return _return_samples(return_type, samples)
+        if HAS_PANDAS is True:
+            return pd.DataFrame.from_records(samples)
+
+        return samples
 
     def generate_sample(self, initial_pos, num_adapt, num_samples, stepsize=None):
         """
diff --git a/pgmpy/sampling/Sampling.py b/pgmpy/sampling/Sampling.py
index 1a31e566..7d2294ea 100644
--- a/pgmpy/sampling/Sampling.py
+++ b/pgmpy/sampling/Sampling.py
@@ -1,16 +1,15 @@
-from collections import namedtuple
 import itertools
+from collections import namedtuple
 
 import networkx as nx
 import numpy as np
+from pandas import DataFrame
+from pgmpy.extern.six.moves import map, range
 
 from pgmpy.factors.discrete import factor_product
 from pgmpy.inference import Inference
 from pgmpy.models import BayesianModel, MarkovChain, MarkovModel
 from pgmpy.utils.mathext import sample_discrete
-from pgmpy.extern.six.moves import map, range
-from pgmpy.sampling import _return_samples
-
 
 State = namedtuple('State', ['var', 'state'])
 
@@ -36,7 +35,7 @@ class BayesianModelSampling(Inference):
         self.topological_order = list(nx.topological_sort(model))
         super(BayesianModelSampling, self).__init__(model)
 
-    def forward_sample(self, size=1, return_type='dataframe'):
+    def forward_sample(self, size=1):
         """
         Generates sample(s) from joint distribution of the bayesian network.
 
@@ -45,16 +44,11 @@ class BayesianModelSampling(Inference):
         size: int
             size of sample to be generated
 
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
-
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        sampled: pandas.DataFrame
             the generated samples
 
-
         Examples
         --------
         >>> from pgmpy.models.BayesianModel import BayesianModel
@@ -68,26 +62,24 @@ class BayesianModelSampling(Inference):
         ...                ['intel', 'diff'], [2, 2])
         >>> student.add_cpds(cpd_d, cpd_i, cpd_g)
         >>> inference = BayesianModelSampling(student)
-        >>> inference.forward_sample(size=2, return_type='recarray')
-        rec.array([(0, 0, 1), (1, 0, 2)], 
-          dtype=[('diff', '<i8'), ('intel', '<i8'), ('grade', '<i8')])
+        >>> inference.forward_sample(2)
+                diff       intel       grade
+        0        1           0          1
+        1        1           0          2
         """
-        types = [(var_name, 'int') for var_name in self.topological_order]
-        sampled = np.zeros(size, dtype=types).view(np.recarray)
-
+        sampled = DataFrame(index=range(size), columns=self.topological_order)
         for node in self.topological_order:
             cpd = self.model.get_cpds(node)
             states = range(self.cardinality[node])
             evidence = cpd.variables[:0:-1]
             if evidence:
                 cached_values = self.pre_compute_reduce(variable=node)
-                evidence = np.vstack([sampled[i] for i in evidence])
-                weights = list(map(lambda t: cached_values[tuple(t)], evidence.T))
+                evidence = sampled.ix[:, evidence].values
+                weights = list(map(lambda t: cached_values[tuple(t)], evidence))
             else:
                 weights = cpd.values
             sampled[node] = sample_discrete(states, weights, size)
-
-        return _return_samples(return_type, sampled)
+        return sampled
 
     def pre_compute_reduce(self, variable):
         variable_cpd = self.model.get_cpds(variable)
@@ -100,7 +92,7 @@ class BayesianModelSampling(Inference):
 
         return cached_values
 
-    def rejection_sample(self, evidence=None, size=1, return_type="dataframe"):
+    def rejection_sample(self, evidence=None, size=1):
         """
         Generates sample(s) from joint distribution of the bayesian network,
         given the evidence.
@@ -111,13 +103,10 @@ class BayesianModelSampling(Inference):
             None if no evidence
         size: int
             size of sample to be generated
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
 
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        sampled: pandas.DataFrame
             the generated samples
 
         Examples
@@ -135,32 +124,26 @@ class BayesianModelSampling(Inference):
         >>> student.add_cpds(cpd_d, cpd_i, cpd_g)
         >>> inference = BayesianModelSampling(student)
         >>> evidence = [State(var='diff', state=0)]
-        >>> inference.rejection_sample(evidence=evidence, size=2, return_type='dataframe')
+        >>> inference.rejection_sample(evidence, 2)
                 intel       diff       grade
         0         0          0          1
         1         0          0          1
         """
         if evidence is None:
             return self.forward_sample(size)
-        types = [(var_name, 'int') for var_name in self.topological_order]
-        sampled = np.zeros(0, dtype=types).view(np.recarray)
+        sampled = DataFrame(columns=self.topological_order)
         prob = 1
-        i = 0
-        while i < size:
-            _size = int(((size - i) / prob) * 1.5)
-            _sampled = self.forward_sample(_size, 'recarray')
-
+        while len(sampled) < size:
+            _size = int(((size - len(sampled)) / prob) * 1.5)
+            _sampled = self.forward_sample(_size)
             for evid in evidence:
-                _sampled = _sampled[_sampled[evid[0]] == evid[1]]
-
+                _sampled = _sampled[_sampled.ix[:, evid[0]] == evid[1]]
             prob = max(len(_sampled) / _size, 0.01)
-            sampled = np.append(sampled, _sampled)[:size]
-
-            i += len(_sampled)
-
-        return _return_samples(return_type, sampled)
+            sampled = sampled.append(_sampled)
+        sampled.reset_index(inplace=True, drop=True)
+        return sampled[:size]
 
-    def likelihood_weighted_sample(self, evidence=None, size=1, return_type="dataframe"):
+    def likelihood_weighted_sample(self, evidence=None, size=1):
         """
         Generates weighted sample(s) from joint distribution of the bayesian
         network, that comply with the given evidence.
@@ -173,13 +156,10 @@ class BayesianModelSampling(Inference):
             None if no evidence
         size: int
             size of sample to be generated
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
 
         Returns
         -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
+        sampled: pandas.DataFrame
             the generated samples with corresponding weights
 
         Examples
@@ -197,40 +177,36 @@ class BayesianModelSampling(Inference):
         >>> student.add_cpds(cpd_d, cpd_i, cpd_g)
         >>> inference = BayesianModelSampling(student)
         >>> evidence = [State('diff', 0)]
-        >>> inference.likelihood_weighted_sample(evidence=evidence, size=2, return_type='recarray')
-        rec.array([(0, 0, 1, 0.6), (0, 0, 2, 0.6)], 
-          dtype=[('diff', '<i8'), ('intel', '<i8'), ('grade', '<i8'), ('_weight', '<f8')])
+        >>> inference.likelihood_weighted_sample(evidence, 2)
+                intel       diff       grade  _weight
+        0         0          0          1        0.6
+        1         1          0          1        0.6
         """
-        types = [(var_name, 'int') for var_name in self.topological_order]
-        types.append(('_weight', 'float'))
-        sampled = np.zeros(size, dtype=types).view(np.recarray)
+        sampled = DataFrame(index=range(size), columns=self.topological_order)
         sampled['_weight'] = np.ones(size)
         evidence_dict = {var: st for var, st in evidence}
-
         for node in self.topological_order:
             cpd = self.model.get_cpds(node)
             states = range(self.cardinality[node])
             evidence = cpd.get_evidence()
-
             if evidence:
-                evidence_values = np.vstack([sampled[i] for i in evidence])
+                evidence_values = sampled.ix[:, evidence].values
                 cached_values = self.pre_compute_reduce(node)
-                weights = list(map(lambda t: cached_values[tuple(t)], evidence_values.T))
+                weights = list(map(lambda t: cached_values[tuple(t)], evidence_values))
                 if node in evidence_dict:
                     sampled[node] = evidence_dict[node]
                     for i in range(size):
-                        sampled['_weight'][i] *= weights[i][evidence_dict[node]]
+                        sampled.loc[i, '_weight'] *= weights[i][evidence_dict[node]]
                 else:
                     sampled[node] = sample_discrete(states, weights)
             else:
                 if node in evidence_dict:
                     sampled[node] = evidence_dict[node]
                     for i in range(size):
-                        sampled['_weight'][i] *= cpd.values[evidence_dict[node]]
+                        sampled.loc[i, '_weight'] *= cpd.values[evidence_dict[node]]
                 else:
                     sampled[node] = sample_discrete(states, cpd.values, size)
-
-        return _return_samples(return_type, sampled)
+        return sampled
 
 
 class GibbsSampling(MarkovChain):
@@ -259,8 +235,10 @@ class GibbsSampling(MarkovChain):
     >>> student.add_nodes_from(['intel', 'sat'])
     >>> student.add_edge('intel', 'sat')
     >>> student.add_cpds(intel_cpd, sat_cpd)
+
     >>> from pgmpy.inference import GibbsSampling
     >>> gibbs_chain = GibbsSampling(student)
+
     Sample from it:
     >>> gibbs_chain.sample(size=3)
        intel  sat
@@ -336,7 +314,7 @@ class GibbsSampling(MarkovChain):
                 kernel[tup] = reduced_factor.values / sum(reduced_factor.values)
             self.transition_models[var] = kernel
 
-    def sample(self, start_state=None, size=1, return_type="dataframe"):
+    def sample(self, start_state=None, size=1):
         """
         Sample from the Markov Chain.
 
@@ -346,14 +324,10 @@ class GibbsSampling(MarkovChain):
             Representing the starting states of the variables. If None is passed, a random start_state is chosen.
         size: int
             Number of samples to be generated.
-        return_type: string (dataframe | recarray)
-            Return type for samples, either of 'dataframe' or 'recarray'.
-            Defaults to 'dataframe'
 
-        Returns
-        -------
-        sampled: A pandas.DataFrame or a numpy.recarray object depending upon return_type argument
-            the generated samples
+        Return Type:
+        ------------
+        pandas.DataFrame
 
         Examples:
         ---------
@@ -365,7 +339,7 @@ class GibbsSampling(MarkovChain):
         >>> factor_cb = DiscreteFactor(['C', 'B'], [2, 2], [5, 6, 7, 8])
         >>> model.add_factors(factor_ab, factor_cb)
         >>> gibbs = GibbsSampling(model)
-        >>> gibbs.sample(size=4, return_tupe='dataframe')
+        >>> gibbs.sample(size=4)
            A  B  C
         0  0  1  1
         1  1  0  0
@@ -377,18 +351,16 @@ class GibbsSampling(MarkovChain):
         elif start_state is not None:
             self.set_start_state(start_state)
 
-        types = [(var_name, 'int') for var_name in self.variables]
-        sampled = np.zeros(size, dtype=types).view(np.recarray)
-        sampled[0] = np.array([st for var, st in self.state])
+        sampled = DataFrame(index=range(size), columns=self.variables)
+        sampled.loc[0] = [st for var, st in self.state]
         for i in range(size - 1):
             for j, (var, st) in enumerate(self.state):
                 other_st = tuple(st for v, st in self.state if var != v)
                 next_st = sample_discrete(list(range(self.cardinalities[var])),
                                           self.transition_models[var][other_st])[0]
                 self.state[j] = State(var, next_st)
-            sampled[i + 1] = np.array([st for var, st in self.state])
-
-        return _return_samples(return_type, sampled)
+            sampled.loc[i + 1] = [st for var, st in self.state]
+        return sampled
 
     def generate_sample(self, start_state=None, size=1):
         """
@@ -413,10 +385,9 @@ class GibbsSampling(MarkovChain):
         [[State(var='C', state=1), State(var='B', state=1), State(var='A', state=0)],
          [State(var='C', state=0), State(var='B', state=1), State(var='A', state=1)]]
         """
-
         if start_state is None and self.state is None:
             self.state = self.random_state()
-        elif start_state is not None:
+        else:
             self.set_start_state(start_state)
 
         for i in range(size):
diff --git a/pgmpy/sampling/__init__.py b/pgmpy/sampling/__init__.py
index 318d2a18..9ca9e2e8 100644
--- a/pgmpy/sampling/__init__.py
+++ b/pgmpy/sampling/__init__.py
@@ -1,5 +1,5 @@
 from .base import (BaseGradLogPDF, GradLogPDFGaussian, LeapFrog,
-                                 ModifiedEuler, BaseSimulateHamiltonianDynamics, _return_samples)
+                                 ModifiedEuler, BaseSimulateHamiltonianDynamics)
 from .HMC import HamiltonianMC, HamiltonianMCDA
 from .NUTS import NoUTurnSampler, NoUTurnSamplerDA
 from .Sampling import GibbsSampling, BayesianModelSampling
@@ -9,7 +9,6 @@ __all__ = ['LeapFrog',
            'BaseSimulateHamiltonianDynamics',
            'BaseGradLogPDF',
            'GradLogPDFGaussian',
-           '_return_samples',
            'HamiltonianMC',
            'HamiltonianMCDA',
            'NoUTurnSampler',
diff --git a/pgmpy/sampling/base.py b/pgmpy/sampling/base.py
index 48df455e..ae92d470 100644
--- a/pgmpy/sampling/base.py
+++ b/pgmpy/sampling/base.py
@@ -1,12 +1,6 @@
-from warnings import warn
-
 import numpy as np
 
 from pgmpy.utils import _check_1d_array_object, _check_length_equal
-from pgmpy.base import HAS_PANDAS
-
-if HAS_PANDAS:
-    import pandas
 
 
 class BaseGradLogPDF(object):
@@ -404,17 +398,3 @@ class ModifiedEuler(BaseSimulateHamiltonianDynamics):
         grad_log, _ = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()
 
         return position_bar, momentum_bar, grad_log
-
-
-def _return_samples(return_type, samples):
-    """
-        A utility function to return samples according to type
-    """
-    if return_type.lower() == "dataframe":
-        if HAS_PANDAS:
-            return pandas.DataFrame.from_records(samples)
-        else:
-            warn("Pandas installation not found. Returning numpy.recarray object")
-            return samples
-    else:
-        return samples
diff --git a/pgmpy/utils/state_name.py b/pgmpy/utils/state_name.py
index 97dd8ffa..01c28c8a 100644
--- a/pgmpy/utils/state_name.py
+++ b/pgmpy/utils/state_name.py
@@ -1,4 +1,3 @@
-from wrapt import decorator
 """
 The file contains decorators to manage the user defined
 variable state names. It maps the internal representaion
@@ -29,16 +28,17 @@ class StateNameInit():
     >>> print(phi.state_names)
     {'speed': ['low', 'medium', 'high'], 'switch': ['on', 'off'], 'time': ['day', 'night']}
     """
-    @decorator
-    def __call__(self, f, instance, args, kwargs):
-        # Case, when no state names dict is provided.
-        if 'state_names' not in kwargs:
-            # instance represents the self parameter of the __init__ method
-            instance.state_names = None
-        else:
-            instance.state_names = kwargs['state_names']
-            del kwargs['state_names']
-        f(*args, **kwargs)
+    def __call__(self, f):
+        def wrapper(*args, **kwargs):
+            # Case, when no state names dict is provided.
+            if 'state_names' not in kwargs:
+                # args[0] represents the self parameter of the __init__ method
+                args[0].state_names = None
+            else:
+                args[0].state_names = kwargs['state_names']
+                del kwargs['state_names']
+            f(*args, **kwargs)
+        return wrapper
 
 
 class StateNameDecorator():
@@ -158,42 +158,44 @@ class StateNameDecorator():
                     arg_val[var] = self.state_names[var].index(arg_val[var])
             return arg_val
 
-    @decorator
-    def __call__(self, f, instance, args, kwargs):
-        # instance represents the self parameter of the __init__ method
-        method_self = instance
-
-        if not method_self.state_names:
-            return f(*args, **kwargs)
-        else:
-            self.state_names = method_self.state_names
-
-        if self.arg and not self.return_val:
-            # If input parameters are in kwargs format
-            if self.arg in kwargs:
-                arg_val = kwargs[self.arg]
-                kwargs[self.arg] = self.get_mapped_value(arg_val)
+    def __call__(self, f):
+        def wrapper(*args, **kwargs):
+            # args[0] represents the self parameter of the __init__ method
+            method_self = args[0]
+
+            if not method_self.state_names:
                 return f(*args, **kwargs)
-            # If input parameters are in args format
             else:
-                for arg_val in args[:]:
-                    mapped_arg_val = self.get_mapped_value(arg_val)
-                    if mapped_arg_val:
-                        mapped_args = list(args)
-                        mapped_args[args.index(arg_val)] = mapped_arg_val
-                return f(*mapped_args, **kwargs)
-
-        elif not self.arg and self.return_val:
-            return_val = f(*args, **kwargs)
-
-            mapped_return_val = self.get_mapped_value(return_val)
-            # If the function returns only one output
-            if mapped_return_val:
-                return mapped_return_val
-
-            # If the function returns more than one output.
-            for ret_val in list(return_val):
-                mapped_ret_val = self.get_mapped_value(ret_val)
-                if ret_val:
-                    return_val[return_val.index(ret_val)] = mapped_ret_val
-            return return_val
+                self.state_names = method_self.state_names
+
+            if self.arg and not self.return_val:
+                # If input parameters are in kwargs format
+                if self.arg in kwargs:
+                    arg_val = kwargs[self.arg]
+                    kwargs[self.arg] = self.get_mapped_value(arg_val)
+                    return f(*args, **kwargs)
+                # If input parameters are in args format
+                else:
+                    for arg_val in args[1:]:
+                        mapped_arg_val = self.get_mapped_value(arg_val)
+                        if mapped_arg_val:
+                            mapped_args = list(args)
+                            mapped_args[args.index(arg_val)] = mapped_arg_val
+                    return f(*mapped_args, **kwargs)
+
+            elif not self.arg and self.return_val:
+                return_val = f(*args, **kwargs)
+
+                mapped_return_val = self.get_mapped_value(return_val)
+                # If the function returns only one output
+                if mapped_return_val:
+                    return mapped_return_val
+
+                # If the function returns more than one output.
+                for ret_val in list(return_val):
+                    mapped_ret_val = self.get_mapped_value(ret_val)
+                    if ret_val:
+                        return_val[return_val.index(ret_val)] = mapped_ret_val
+                return return_val
+
+        return wrapper
diff --git a/requirements-dev.txt b/requirements-dev.txt
index ab15a26a..206656f9 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -5,4 +5,3 @@ pandas
 pyparsing
 mock
 nose
-wrapt
diff --git a/requirements.txt b/requirements.txt
index fb4c749e..6eed30de 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,5 +3,3 @@ numpy==1.11.1
 scipy==0.18.0
 pandas==0.18.1
 pyparsing==2.1.8
-wrapt==1.10.8
-
diff --git a/setup.py b/setup.py
index cace146c..4c3ef836 100644
--- a/setup.py
+++ b/setup.py
@@ -10,9 +10,7 @@ setup(
     author=open("AUTHORS.rst").read(),
     url="https://github.com/pgmpy/pgmpy",
     classifiers=[
-        "Programming Language :: Python :: 3.5",
-        "Programming Language :: Python :: 3.4",
-        "Programming Language :: Python :: 2.7",
+        "Programming Language :: Python :: 3",
         "Intended Audience :: Developers",
         "Operating System :: Unix",
         "Operating System :: POSIX",


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_BaseEstimator.py b/pgmpy/tests/test_estimators/test_BaseEstimator.py
index b14d4b05..37ec9b9a 100644
--- a/pgmpy/tests/test_estimators/test_BaseEstimator.py
+++ b/pgmpy/tests/test_estimators/test_BaseEstimator.py
@@ -40,12 +40,12 @@ class TestBaseEstimator(unittest.TestCase):
     def test_test_conditional_independence_titanic(self):
         est = BaseEstimator(self.titanic_data)
 
-        np.testing.assert_almost_equal(est.test_conditional_independence('Embarked', 'Sex'),
-                                       (13.355630515001746, 0.020264556044311655, True))
-        np.testing.assert_almost_equal(est.test_conditional_independence('Pclass', 'Survived', ['Embarked']),
-                                       (96.403283942888635, 4.1082315854166553e-13, True))
-        np.testing.assert_almost_equal(est.test_conditional_independence('Embarked', 'Survived', ["Sex", "Pclass"]),
-                                       (21.537481934494085, 0.96380273702382602, True))
+        self.assertEqual(est.test_conditional_independence('Embarked', 'Sex'),
+                         (13.355630515001746, 0.020264556044311655, True))
+        self.assertEqual(est.test_conditional_independence('Pclass', 'Survived', ['Embarked']),
+                         (96.403283942888635, 4.1082315854166553e-13, True))
+        self.assertEqual(est.test_conditional_independence('Embarked', 'Survived', ["Sex", "Pclass"]),
+                         (21.537481934494085, 0.96380273702382602, True))
         # insufficient data test commented out, because generates warning
         # self.assertEqual(est.test_conditional_independence('Sex', 'Survived', ["Age", "Embarked"]),
         #                 (235.51133052530713, 0.99999999683394869, False))
diff --git a/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
index b1a5cef0..99a9bbda 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_Canonical_Factor.py
@@ -96,19 +96,14 @@ class TestJGDMethods(unittest.TestCase):
         np_test.assert_almost_equal(self.phi1.h, np.array([[11.7], [-1.0]]))
         np_test.assert_almost_equal(self.phi1.g, -21.95)
 
-        self.phi1 = self.phi3.copy()
+        self.phi1 = self.phi3
         self.phi1.reduce([('x1', 4), ('x2', 1.23)])
         self.assertEqual(self.phi1.variables, ['x3'])
         np_test.assert_almost_equal(self.phi1.K, np.array([[4.0]]))
         np_test.assert_almost_equal(self.phi1.h, np.array([[1.46]]))
         np_test.assert_almost_equal(self.phi1.g, 0.8752)
 
-        self.phi1 = self.phi3.copy()
-        self.phi1.reduce([('x2', 1.23), ('x1', 4)])
-        self.assertEqual(self.phi1.variables, ['x3'])
-        np_test.assert_almost_equal(self.phi1.K, np.array([[4.0]]))
-        np_test.assert_almost_equal(self.phi1.h, np.array([[1.46]]))
-        np_test.assert_almost_equal(self.phi1.g, 0.8752)
+        self.phi1 = self.phi3
 
     def test_marginalize(self):
         phi = self.phi1.marginalize(['x1'], inplace=False)
diff --git a/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py b/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
index d899334d..5fef3df8 100644
--- a/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
+++ b/pgmpy/tests/test_factors/test_continuous/test_JointGaussianDistribution.py
@@ -183,7 +183,7 @@ class TestJGDMethods(unittest.TestCase):
                                    np.asarray([[4, -4], [-4, 7]], dtype=float))
         self.assertEqual(self.phi1._precision_matrix, None)
 
-        self.phi1 = self.phi3.copy()
+        self.phi1 = self.phi3
         self.phi1.reduce([('x1', 3), ('x2', 1)])
         self.assertEqual(self.phi1.variables, ['x3'])
         np_test.assert_array_equal(self.phi1.mean, np.asarray([[0]], dtype=float))
@@ -191,13 +191,7 @@ class TestJGDMethods(unittest.TestCase):
                                    np.asarray([[3]], dtype=float))
         self.assertEqual(self.phi1._precision_matrix, None)
 
-        self.phi1 = self.phi3.copy()
-        self.phi1.reduce([('x2', 1), ('x1', 3)])
-        self.assertEqual(self.phi1.variables, ['x3'])
-        np_test.assert_array_equal(self.phi1.mean, np.asarray([[0]], dtype=float))
-        np_test.assert_array_equal(self.phi1.covariance,
-                                   np.asarray([[3]], dtype=float))
-        self.assertEqual(self.phi1._precision_matrix, None)
+        self.phi1 = self.phi3
 
     def test_normalize(self):
         phi = self.phi1.copy()
diff --git a/pgmpy/tests/test_models/test_MarkovChain.py b/pgmpy/tests/test_models/test_MarkovChain.py
index fb5e6da1..87e57850 100644
--- a/pgmpy/tests/test_models/test_MarkovChain.py
+++ b/pgmpy/tests/test_models/test_MarkovChain.py
@@ -5,6 +5,7 @@ import numpy as np
 from pandas import DataFrame
 from mock import patch, call
 
+
 from pgmpy.factors.discrete import State
 from pgmpy.models import MarkovChain as MC
 from pgmpy.extern.six.moves import range, zip
@@ -229,6 +230,24 @@ class TestMarkovChain(unittest.TestCase):
         probabilites = model.prob_from_sample([State('a', 1), State('b', 0)])
         self.assertEqual(list(probabilites), [1] * 50 + [0] * 50)
 
+    def test_is_stationarity_success(self):
+        model = MC(['intel', 'diff'], [2, 3])
+        model.set_start_state([State('intel', 0), State('diff', 2)])
+        intel_tm = {0: {0: 0.25, 1: 0.75}, 1: {0: 0.5, 1: 0.5}}
+        model.add_transition_model('intel', intel_tm)
+        diff_tm = {0: {0: 0.1, 1: 0.5, 2: 0.4}, 1: {0: 0.2, 1: 0.2, 2: 0.6}, 2: {0: 0.7, 1: 0.15, 2: 0.15}}
+        model.add_transition_model('diff', diff_tm)
+        self.assertTrue(model.is_stationarity)
+
+    def test_is_stationarity_failure(self):
+        model = MC(['intel', 'diff'], [2, 3])
+        model.set_start_state([State('intel', 0), State('diff', 2)])
+        intel_tm = {0: {0: 0.25, 1: 0.75}, 1: {0: 0.5, 1: 0.5}}
+        model.add_transition_model('intel', intel_tm)
+        diff_tm = {0: {0: 0.1, 1: 0.5, 2: 0.4}, 1: {0: 0.2, 1: 0.2, 2: 0.6}, 2: {0: 0.7, 1: 0.15, 2: 0.15}}
+        model.add_transition_model('diff', diff_tm)
+        self.assertFalse(model.is_stationarity(0.002, None))
+
     @patch.object(sys.modules["pgmpy.models.MarkovChain"], "sample_discrete")
     def test_generate_sample(self, sample_discrete):
         model = MC(['a', 'b'], [2, 2])
diff --git a/pgmpy/tests/test_models/test_NaiveBayes.py b/pgmpy/tests/test_models/test_NaiveBayes.py
index 3a74d382..3f1b064c 100644
--- a/pgmpy/tests/test_models/test_NaiveBayes.py
+++ b/pgmpy/tests/test_models/test_NaiveBayes.py
@@ -129,10 +129,13 @@ class TestNaiveBayesMethods(unittest.TestCase):
                               ('d', 's')])
 
     def test_local_independencies(self):
-        self.assertEqual(self.G1.local_independencies('a'), Independencies())
-        self.assertEqual(self.G1.local_independencies('b'), Independencies(['b', ['e', 'c', 'd'], 'a']))
-        self.assertEqual(self.G1.local_independencies('c'), Independencies(['c', ['e', 'b', 'd'], 'a']))
-        self.assertEqual(self.G1.local_independencies('d'), Independencies(['d', ['b', 'c', 'e'], 'a']))
+        self.assertListEqual(self.G1.local_independencies('a'), [None])
+        self.assertListEqual(self.G1.local_independencies('b'),
+                             [Independencies(['b', ['e', 'c', 'd'], 'a'])])
+        self.assertListEqual(self.G1.local_independencies('c'),
+                             [Independencies(['c', ['e', 'b', 'd'], 'a'])])
+        self.assertListEqual(self.G1.local_independencies('d'),
+                             [Independencies(['d', ['b', 'c', 'e'], 'a'])])
 
     def test_active_trail_nodes(self):
         self.assertListEqual(sorted(self.G2.active_trail_nodes('d')), ['d', 'g', 'l', 's'])
@@ -146,12 +149,6 @@ class TestNaiveBayesMethods(unittest.TestCase):
         self.assertListEqual(sorted(self.G2.active_trail_nodes('s', observed=['g', 'l'])), ['d', 's'])
         self.assertListEqual(sorted(self.G2.active_trail_nodes('s', observed=['d', 'l'])), ['s'])
 
-    def test_get_ancestors_of(self):
-        self.assertListEqual(sorted(self.G1._get_ancestors_of('b')), ['a', 'b'])
-        self.assertListEqual(sorted(self.G1._get_ancestors_of('e')), ['a', 'e'])
-        self.assertListEqual(sorted(self.G1._get_ancestors_of('a')), ['a'])
-        self.assertListEqual(sorted(self.G1._get_ancestors_of(['b', 'e'])), ['a', 'b', 'e'])
-
     def tearDown(self):
         del self.G1
         del self.G2
diff --git a/pgmpy/tests/test_sampling/test_continuous_sampling.py b/pgmpy/tests/test_sampling/test_continuous_sampling.py
index bb115be2..4dcb8c16 100644
--- a/pgmpy/tests/test_sampling/test_continuous_sampling.py
+++ b/pgmpy/tests/test_sampling/test_continuous_sampling.py
@@ -123,10 +123,8 @@ class TestNUTSInference(unittest.TestCase):
 
     def test_sampling(self):
         np.random.seed(1010101)
-        samples = self.nuts_sampler.sample(initial_pos=[-0.4, 1, 3.6], num_adapt=0, num_samples=10000,
-                                            return_type='recarray')
-        sample_array = np.array([samples[var_name] for var_name in self.test_model.variables])
-        sample_covariance = np.cov(sample_array)
+        samples = self.nuts_sampler.sample(initial_pos=[-0.4, 1, 3.6], num_adapt=0, num_samples=10000)
+        sample_covariance = np.cov(samples.values.T)
         self.assertTrue(np.linalg.norm(sample_covariance - self.test_model.covariance) < 3)
 
         np.random.seed(1210161)