Repository: pgmpy/pgmpy
Pull request #830 (2017-01-21 22:03:18)

Add tests for _get_ancestors
Fixes PR #522.
Adds tests for _get_ancestors method in Bayesian Model.
Also adds a way to check if the argument given to _get_ancestors method is present in the graph.
ping @yashu-seth @ankurankan


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/estimators/HillClimbSearch.py b/pgmpy/estimators/HillClimbSearch.py
index 09e57d1c..fe47c310 100644
--- a/pgmpy/estimators/HillClimbSearch.py
+++ b/pgmpy/estimators/HillClimbSearch.py
@@ -136,7 +136,7 @@ class HillClimbSearch(StructureEstimator):
         >>> est.estimate(max_indegree=1).edges()
         [('J', 'A'), ('B', 'J')]
         """
-        epsilon = 1e-8
+
         nodes = self.state_names.keys()
         if start is None:
             start = BayesianModel()
@@ -156,7 +156,7 @@ class HillClimbSearch(StructureEstimator):
                     best_operation = operation
                     best_score_delta = score_delta
 
-            if best_operation is None or best_score_delta < epsilon:
+            if best_operation is None:
                 break
             elif best_operation[0] == '+':
                 current_model.add_edge(*best_operation[1])
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/factors/discrete/CPD.py b/pgmpy/factors/discrete/CPD.py
index 49ad297d..4b678be6 100644
--- a/pgmpy/factors/discrete/CPD.py
+++ b/pgmpy/factors/discrete/CPD.py
@@ -94,7 +94,7 @@ class TabularCPD(DiscreteFactor):
 
     Public Methods
     --------------
-    get_values()
+    get_cpd()
     marginalize([variables_list])
     normalize()
     reduce([values_list])
@@ -146,18 +146,15 @@ class TabularCPD(DiscreteFactor):
 
         return var_str + evidence_str + ') at {address}>'.format(address=hex(id(self)))
 
-    def get_values(self):
+    def get_cpd(self):
         """
         Returns the cpd
-        
-        Examples
-        --------
-        >>> from pgmpy.factors.discrete import TabularCPD
+        >>> from pgmpy.factors import TabularCPD
         >>> cpd = TabularCPD('grade', 3, [[0.1, 0.1],
         ...                               [0.1, 0.1],
         ...                               [0.8, 0.8]],
         ...                  evidence='evi1', evidence_card=2)
-        >>> cpd.get_values()
+        >>> cpd.get_cpd()
         array([[ 0.1,  0.1],
                [ 0.1,  0.1],
                [ 0.8,  0.8]])
@@ -205,7 +202,7 @@ class TabularCPD(DiscreteFactor):
         else:
             variable_array = [['{s}_{d}'.format(s=self.variable, d=i) for i in range(self.variable_card)]]
         # Stack with data
-        labeled_rows = np.hstack((np.array(variable_array).T, self.get_values())).tolist()
+        labeled_rows = np.hstack((np.array(variable_array).T, self.get_cpd())).tolist()
         # No support for multi-headers in tabulate
         cdf_str = tabulate(headers_list + labeled_rows, tablefmt=tablefmt)
         return cdf_str
@@ -236,7 +233,7 @@ class TabularCPD(DiscreteFactor):
         """
         evidence = self.variables[1:] if len(self.variables) > 1 else None
         evidence_card = self.cardinality[1:] if len(self.variables) > 1 else None
-        return TabularCPD(self.variable, self.variable_card, self.get_values(),
+        return TabularCPD(self.variable, self.variable_card, self.get_cpd(),
                           evidence, evidence_card)
 
     def normalize(self, inplace=True):
@@ -256,12 +253,12 @@ class TabularCPD(DiscreteFactor):
         ...                        [[0.7, 0.2, 0.6, 0.2],[0.4, 0.4, 0.4, 0.8]],
         ...                        ['intel', 'diff'], [2, 2])
         >>> cpd_table.normalize()
-        >>> cpd_table.get_values()
+        >>> cpd_table.get_cpd()
         array([[ 0.63636364,  0.33333333,  0.6       ,  0.2       ],
                [ 0.36363636,  0.66666667,  0.4       ,  0.8       ]])
         """
         tabular_cpd = self if inplace else self.copy()
-        cpd = tabular_cpd.get_values()
+        cpd = tabular_cpd.get_cpd()
         tabular_cpd.values = (cpd / cpd.sum(axis=0)).reshape(tabular_cpd.cardinality)
         if not inplace:
             return tabular_cpd
@@ -286,7 +283,7 @@ class TabularCPD(DiscreteFactor):
         ...                        [[0.7, 0.6, 0.6, 0.2],[0.3, 0.4, 0.4, 0.8]],
         ...                        ['intel', 'diff'], [2, 2])
         >>> cpd_table.marginalize(['diff'])
-        >>> cpd_table.get_values()
+        >>> cpd_table.get_cpd()
         array([[ 0.65,  0.4 ],
                 [ 0.35,  0.6 ]])
         """
@@ -322,7 +319,7 @@ class TabularCPD(DiscreteFactor):
         ...                        [[0.7, 0.6, 0.6, 0.2],[0.3, 0.4, 0.4, 0.8]],
         ...                        ['intel', 'diff'], [2, 2])
         >>> cpd_table.reduce([('diff', 0)])
-        >>> cpd_table.get_values()
+        >>> cpd_table.get_cpd()
         array([[ 0.7,  0.6],
                [ 0.3,  0.4]])
         """
@@ -459,12 +456,12 @@ class TabularCPD(DiscreteFactor):
                     variables = [self.variables[0]] + new_order
                     cardinality = [self.variable_card] + [card_map[var] for var in new_order]
                     super(TabularCPD, self).__init__(variables, cardinality, new_values.flatten('C'))
-                    return self.get_values()
+                    return self.get_cpd()
                 else:
                     return new_values.reshape(self.cardinality[0], np.prod([card_map[var] for var in new_order]))
             else:
                 warn("Same ordering provided as current")
-                return self.get_values()
+                return self.get_cpd()
 
     def get_evidence(self):
         return self.variables[:0:-1]
diff --git a/pgmpy/models/BayesianModel.py b/pgmpy/models/BayesianModel.py
index 413abe6f..6a82e988 100644
--- a/pgmpy/models/BayesianModel.py
+++ b/pgmpy/models/BayesianModel.py
@@ -361,7 +361,7 @@ class BayesianModel(DirectedGraph):
 
     def _get_ancestors_of(self, obs_nodes_list):
         """
-        Returns a list of all ancestors of all the observed nodes including the
+        Returns a dictionary of all ancestors of all the observed nodes including the
         node itself.
 
         Parameters
@@ -382,6 +382,10 @@ class BayesianModel(DirectedGraph):
         if not isinstance(obs_nodes_list, (list, tuple)):
             obs_nodes_list = [obs_nodes_list]
 
+        for node in obs_nodes_list:
+            if node not in self.nodes():
+                raise ValueError('Node {s} not in not in graph'.format(s=node))
+
         ancestors_list = set()
         nodes_list = set(obs_nodes_list)
         while nodes_list:
@@ -456,7 +460,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 +475,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):
             """
@@ -710,20 +718,11 @@ class BayesianModel(DirectedGraph):
         >>> predict_data.drop('E', axis=1, inplace=True)
         >>> y_pred = model.predict(predict_data)
         >>> y_pred
-            E
-        800 0
-        801 1
-        802 1
-        803 1
-        804 0
-        ... ...
-        993 0
-        994 0
-        995 1
-        996 1
-        997 0
-        998 0
-        999 0
+        array([0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
+               1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1,
+               1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1,
+               1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1,
+               1, 1, 1, 0, 0, 0, 1, 0])
         """
         from pgmpy.inference import VariableElimination
 
@@ -731,85 +730,18 @@ class BayesianModel(DirectedGraph):
             raise ValueError("No variable missing in data. Nothing to predict")
 
         elif set(data.columns) - set(self.nodes()):
-            raise ValueError("Data has variables which are not in the model")
+            raise ValueError("data has variables which are not in the model")
 
         missing_variables = set(self.nodes()) - set(data.columns)
         pred_values = defaultdict(list)
 
-        # Send state_names dict from one of the estimated CPDs to the inference class.
-        model_inference = VariableElimination(self, state_names=self.get_cpds()[0].state_names)
+        model_inference = VariableElimination(self)
         for index, data_point in data.iterrows():
             states_dict = model_inference.map_query(variables=missing_variables, evidence=data_point.to_dict())
             for k, v in states_dict.items():
                 pred_values[k].append(v)
         return pd.DataFrame(pred_values, index=data.index)
 
-    def predict_probability(self, data):
-        """
-        Predicts probabilities of all states of the missing variables.
-
-        Parameters
-        ----------
-        data : pandas DataFrame object
-            A DataFrame object with column names same as the variables in the model.
-
-        Examples
-        --------
-        >>> import numpy as np
-        >>> import pandas as pd
-        >>> from pgmpy.models import BayesianModel
-        >>> values = pd.DataFrame(np.random.randint(low=0, high=2, size=(100, 5)),
-        ...                       columns=['A', 'B', 'C', 'D', 'E'])
-        >>> train_data = values[:80]
-        >>> predict_data = values[80:]
-        >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D'), ('B', 'E')])
-        >>> model.fit(values)
-        >>> predict_data = predict_data.copy()
-        >>> predict_data.drop('B', axis=1, inplace=True)
-        >>> y_prob = model.predict_probability(predict_data)
-        >>> y_prob 
-            B_0         B_1
-        80  0.439178    0.560822
-        81  0.581970    0.418030
-        82  0.488275    0.511725
-        83  0.581970    0.418030
-        84  0.510794    0.489206
-        85  0.439178    0.560822
-        86  0.439178    0.560822
-        87  0.417124    0.582876
-        88  0.407978    0.592022
-        89  0.429905    0.570095
-        90  0.581970    0.418030
-        91  0.407978    0.592022
-        92  0.429905    0.570095
-        93  0.429905    0.570095
-        94  0.439178    0.560822
-        95  0.407978    0.592022
-        96  0.559904    0.440096
-        97  0.417124    0.582876
-        98  0.488275    0.511725
-        99  0.407978    0.592022
-        """
-        from pgmpy.inference import VariableElimination
-
-        if set(data.columns) == set(self.nodes()):
-            raise ValueError("No variable missing in data. Nothing to predict")
-
-        elif set(data.columns) - set(self.nodes()):
-            raise ValueError("Data has variables which are not in the model")
-
-        missing_variables = set(self.nodes()) - set(data.columns)
-        pred_values = defaultdict(list)
-
-        model_inference = VariableElimination(self)
-        for index, data_point in data.iterrows():
-            states_dict = model_inference.query(variables=missing_variables, evidence=data_point.to_dict())
-            for k, v in states_dict.items():
-                for l in range(len(v.values)):
-                    state = self.get_cpds(k).state_names[k][l]
-                    pred_values[k + '_' + str(state)].append(v.values[l])
-        return pd.DataFrame(pred_values, index=data.index)
-
     def get_factorized_product(self, latex=False):
         # TODO: refer to IMap class for explanation why this is not implemented.
         pass
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/readwrite/BIF.py b/pgmpy/readwrite/BIF.py
index fac86ab9..63ae4bf2 100644
--- a/pgmpy/readwrite/BIF.py
+++ b/pgmpy/readwrite/BIF.py
@@ -61,7 +61,7 @@ class BIFReader(object):
         self.variable_states = self.get_states()
         self.variable_properties = self.get_property()
         self.variable_parents = self.get_parents()
-        self.variable_cpds = self.get_values()
+        self.variable_cpds = self.get_cpd()
         self.variable_edges = self.get_edges()
 
     def get_variable_grammar(self):
@@ -214,7 +214,7 @@ class BIFReader(object):
             variable_parents[names[0]] = names[1:]
         return variable_parents
 
-    def get_values(self):
+    def get_cpd(self):
         """
         Returns the CPD of the variables present in the network
 
@@ -222,7 +222,7 @@ class BIFReader(object):
         --------
         >>> from pgmpy.readwrite import BIFReader
         >>> reader = BIFReader("bif_test.bif")
-        >>> reader.get_values()
+        >>> reader.get_cpd()
         {'bowel-problem': np.array([[0.01],
                                     [0.99]]),
         'dog-out': np.array([[0.99, 0.97, 0.9, 0.3],
diff --git a/pgmpy/readwrite/XMLBIF.py b/pgmpy/readwrite/XMLBIF.py
index 3dc5bb2d..94c7dc92 100644
--- a/pgmpy/readwrite/XMLBIF.py
+++ b/pgmpy/readwrite/XMLBIF.py
@@ -60,7 +60,7 @@ class XMLBIFReader(object):
         self.variable_parents = self.get_parents()
         self.edge_list = self.get_edges()
         self.variable_states = self.get_states()
-        self.variable_CPD = self.get_values()
+        self.variable_CPD = self.get_cpd()
         self.variable_property = self.get_property()
 
     def get_variables(self):
@@ -129,14 +129,14 @@ class XMLBIFReader(object):
                             for definition in self.network.findall('DEFINITION')}
         return variable_parents
 
-    def get_values(self):
+    def get_cpd(self):
         """
         Returns the CPD of the variables present in the network
 
         Examples
         --------
         >>> reader = XMLBIF.XMLBIFReader("xmlbif_test.xml")
-        >>> reader.get_values()
+        >>> reader.get_cpd()
         {'bowel-problem': array([[ 0.01],
                                  [ 0.99]]),
          'dog-out': array([[ 0.99,  0.01,  0.97,  0.03],
@@ -239,7 +239,7 @@ class XMLBIFWriter(object):
         self.states = self.get_states()
         self.properties = self.get_properties()
         self.definition = self.get_definition()
-        self.tables = self.get_values()
+        self.tables = self.get_cpd()
 
     def __str__(self):
         """
@@ -396,7 +396,7 @@ class XMLBIFWriter(object):
 
         return definition_tag
 
-    def get_values(self):
+    def get_cpd(self):
         """
         Add Table to XMLBIF.
 
@@ -407,7 +407,7 @@ class XMLBIFWriter(object):
         Examples
         -------
         >>> writer = XMLBIFWriter(model)
-        >>> writer.get_values()
+        >>> writer.get_cpd()
         {'dog-out': <Element TABLE at 0x7f240726f3c8>,
          'light-on': <Element TABLE at 0x7f240726f488>,
          'bowel-problem': <Element TABLE at 0x7f240726f388>,
@@ -420,7 +420,7 @@ class XMLBIFWriter(object):
         for cpd in cpds:
             table_tag[cpd.variable] = etree.SubElement(definition_tag[cpd.variable], "TABLE")
             table_tag[cpd.variable].text = ''
-            for val in cpd.get_values().ravel(order="F"):
+            for val in cpd.get_cpd().ravel(order="F"):
                 table_tag[cpd.variable].text += str(val) + ' '
 
         return table_tag
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_factors/test_discrete/test_Factor.py b/pgmpy/tests/test_factors/test_discrete/test_Factor.py
index deda235b..b8b17192 100644
--- a/pgmpy/tests/test_factors/test_discrete/test_Factor.py
+++ b/pgmpy/tests/test_factors/test_discrete/test_Factor.py
@@ -616,30 +616,30 @@ class TestTabularCPDMethods(unittest.TestCase):
 
     def test_copy(self):
         copy_cpd = self.cpd.copy()
-        np_test.assert_array_equal(self.cpd.get_values(), copy_cpd.get_values())
+        np_test.assert_array_equal(self.cpd.get_cpd(), copy_cpd.get_cpd())
 
     def test_copy_original_safe(self):
         copy_cpd = self.cpd.copy()
         copy_cpd.reorder_parents(['diff', 'intel'])
-        np_test.assert_array_equal(self.cpd.get_values(), np.array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
+        np_test.assert_array_equal(self.cpd.get_cpd(), np.array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
                                                                  [0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
                                                                  [0.8, 0.8, 0.8, 0.8, 0.8, 0.8]]))
 
     def test_reduce_1(self):
         self.cpd.reduce([('diff', 0)])
-        np_test.assert_array_equal(self.cpd.get_values(), np.array([[0.1, 0.1, 0.1],
+        np_test.assert_array_equal(self.cpd.get_cpd(), np.array([[0.1, 0.1, 0.1],
                                                                  [0.1, 0.1, 0.1],
                                                                  [0.8, 0.8, 0.8]]))
 
     def test_reduce_2(self):
         self.cpd.reduce([('intel', 0)])
-        np_test.assert_array_equal(self.cpd.get_values(), np.array([[0.1, 0.1],
+        np_test.assert_array_equal(self.cpd.get_cpd(), np.array([[0.1, 0.1],
                                                                  [0.1, 0.1],
                                                                  [0.8, 0.8]]))
 
     def test_reduce_3(self):
         self.cpd.reduce([('intel', 0), ('diff', 0)])
-        np_test.assert_array_equal(self.cpd.get_values(), np.array([[0.1],
+        np_test.assert_array_equal(self.cpd.get_cpd(), np.array([[0.1],
                                                                  [0.1],
                                                                  [0.8]]))
 
@@ -653,8 +653,8 @@ class TestTabularCPDMethods(unittest.TestCase):
         self.cpd.reduce([('diff', 1)])
         np_test.assert_array_almost_equal(self.cpd.values, copy_cpd.values)
 
-    def test_get_values(self):
-        np_test.assert_array_equal(self.cpd.get_values(), np.array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
+    def test_get_cpd(self):
+        np_test.assert_array_equal(self.cpd.get_cpd(), np.array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
                                                                  [0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
                                                                  [0.8, 0.8, 0.8, 0.8, 0.8, 0.8]]))
 
@@ -662,7 +662,7 @@ class TestTabularCPDMethods(unittest.TestCase):
         new_vals = self.cpd2.reorder_parents(['B', 'A', 'C'])
         np_test.assert_array_equal(new_vals, np.array([[0.9, 0.3, 0.8, 0.8, 0.9, 0.3, 0.4, 0.4],
                                                        [0.1, 0.7, 0.2, 0.2, 0.1, 0.7, 0.6, 0.6]]))
-        np_test.assert_array_equal(self.cpd2.get_values(), np.array([[0.9, 0.3, 0.8, 0.8, 0.9, 0.3, 0.4, 0.4],
+        np_test.assert_array_equal(self.cpd2.get_cpd(), np.array([[0.9, 0.3, 0.8, 0.8, 0.9, 0.3, 0.4, 0.4],
                                                                   [0.1, 0.7, 0.2, 0.2, 0.1, 0.7, 0.6, 0.6]]))
 
     def test_reorder_parents(self):
@@ -672,7 +672,7 @@ class TestTabularCPDMethods(unittest.TestCase):
 
     def test_reorder_parents_no_effect(self):
         self.cpd2.reorder_parents(['C', 'A', 'B'], inplace=False)
-        np_test.assert_array_equal(self.cpd2.get_values(), np.array([[0.9, 0.3, 0.9, 0.3, 0.8, 0.8, 0.4, 0.4],
+        np_test.assert_array_equal(self.cpd2.get_cpd(), np.array([[0.9, 0.3, 0.9, 0.3, 0.8, 0.8, 0.4, 0.4],
                                                                   [0.1, 0.7, 0.1, 0.7, 0.2, 0.2, 0.6, 0.6]]))
 
     def test_reorder_parents_warning(self):
@@ -680,7 +680,7 @@ class TestTabularCPDMethods(unittest.TestCase):
             warnings.simplefilter("always")
             self.cpd2.reorder_parents(['A', 'B', 'C'], inplace=False)
             assert("Same ordering provided as current" in str(w[-1].message))
-            np_test.assert_array_equal(self.cpd2.get_values(), np.array([[0.9, 0.3, 0.9, 0.3, 0.8, 0.8, 0.4, 0.4],
+            np_test.assert_array_equal(self.cpd2.get_cpd(), np.array([[0.9, 0.3, 0.9, 0.3, 0.8, 0.8, 0.4, 0.4],
                                                                       [0.1, 0.7, 0.1, 0.7, 0.2, 0.2, 0.6, 0.6]]))
 
     def tearDown(self):
diff --git a/pgmpy/tests/test_models/test_BayesianModel.py b/pgmpy/tests/test_models/test_BayesianModel.py
index cc66420c..43ac8ab2 100644
--- a/pgmpy/tests/test_models/test_BayesianModel.py
+++ b/pgmpy/tests/test_models/test_BayesianModel.py
@@ -116,6 +116,7 @@ class TestBayesianModelMethods(unittest.TestCase):
                                                    [0.8, 0.8, 0.8, 0.8, 0.8, 0.8]],
                                evidence=['diff', 'intel'], evidence_card=[2, 3])
         self.G1.add_cpds(diff_cpd, intel_cpd, grade_cpd)
+        self.G2 = BayesianModel([('d', 'g'), ('g', 'l'), ('i', 'g'), ('i', 'l')])
 
     def test_moral_graph(self):
         moral_graph = self.G.moralize()
@@ -132,6 +133,17 @@ class TestBayesianModelMethods(unittest.TestCase):
             self.assertTrue(edge in [('a', 'b'), ('c', 'b'), ('d', 'a'), ('d', 'b'), ('d', 'e')] or
                             (edge[1], edge[0]) in [('a', 'b'), ('c', 'b'), ('d', 'a'), ('d', 'b'), ('d', 'e')])
 
+    def test_get_ancestors_of_success(self):
+        ancenstors1 = self.G2._get_ancestors_of('g')
+        ancenstors2 = self.G2._get_ancestors_of('d')
+        ancenstors3 = self.G2._get_ancestors_of(['i', 'l'])
+        self.assertEqual(ancenstors1, {'d', 'i', 'g'})
+        self.assertEqual(ancenstors2, {'d'})
+        self.assertEqual(ancenstors3, {'g', 'i', 'l', 'd'})
+
+    def test_get_ancestors_of_failure(self):
+        self.assertRaises(ValueError, self.G2._get_ancestors_of, 'h')
+
     def test_local_independencies(self):
         self.assertEqual(self.G.local_independencies('a'), Independencies(['a', ['b', 'c']]))
         self.assertEqual(self.G.local_independencies('c'), Independencies(['c', ['a', 'd', 'e'], 'b']))
@@ -398,10 +410,6 @@ class TestBayesianModelFitPredict(unittest.TestCase):
                                         'C': [1, 1, np.NaN],
                                         'D': [np.NaN, 'Y', np.NaN]})
 
-        # data_link - "https://www.kaggle.com/c/titanic/download/train.csv"
-        self.titanic_data = pd.read_csv('pgmpy/tests/test_estimators/testdata/titanic_train.csv', dtype=str)
-        self.titanic_data2 = self.titanic_data[["Survived", "Sex", "Pclass"]]
-
     def test_bayesian_fit(self):
         print(isinstance(BayesianEstimator, BaseEstimator))
         print(isinstance(MaximumLikelihoodEstimator, BaseEstimator))
@@ -430,35 +438,9 @@ class TestBayesianModelFitPredict(unittest.TestCase):
             value = value.reindex(sorted(value.index)).values
             np_test.assert_array_equal(cpd.values, value)
 
-    def test_predict(self):
-        titanic = BayesianModel()
-        titanic.add_edges_from([("Sex", "Survived"), ("Pclass", "Survived")])
-        titanic.fit(self.titanic_data2[500:])
-
-        p1 = titanic.predict(self.titanic_data2[["Sex", "Pclass"]][:30])
-        p2 = titanic.predict(self.titanic_data2[["Survived", "Pclass"]][:30])
-        p3 = titanic.predict(self.titanic_data2[["Survived", "Sex"]][:30])
-
-        p1_res =  np.array(['0', '1', '0', '1', '0', '0', '0', '0', '0', '1', '0', '1', '0',
-                            '0', '0', '1', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0',
-                            '0', '0', '0', '0'])
-        p2_res = np.array(['male', 'female', 'female', 'female', 'male', 'male', 'male',
-                           'male', 'female', 'female', 'female', 'female', 'male', 'male',
-                           'male', 'female', 'male', 'female', 'male', 'female', 'male',
-                           'female', 'female', 'female', 'male', 'female', 'male', 'male',
-                           'female', 'male'])
-        p3_res = np.array(['3', '1', '1', '1', '3', '3', '3', '3', '1', '1', '1', '1', '3',
-                           '3', '3', '1', '3', '1', '3', '1', '3', '1', '1', '1', '3', '1',
-                           '3', '3', '1', '3'])
-
-        np_test.assert_array_equal(p1.values.ravel(), p1_res)
-        np_test.assert_array_equal(p2.values.ravel(), p2_res)
-        np_test.assert_array_equal(p3.values.ravel(), p3_res)
-
     def test_connected_predict(self):
         np.random.seed(42)
-        values = pd.DataFrame(np.array(np.random.randint(low=0, high=2, size=(1000, 5)),
-                                       dtype=str),
+        values = pd.DataFrame(np.random.randint(low=0, high=2, size=(1000, 5)),
                               columns=['A', 'B', 'C', 'D', 'E'])
         fit_data = values[:800]
         predict_data = values[800:].copy()
@@ -481,40 +463,7 @@ class TestBayesianModelFitPredict(unittest.TestCase):
                                              0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0,
                                              1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
                                              1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1,
-                                             1, 1, 1, 0], dtype=str))
-
-    def test_connected_predict_probability(self):
-        np.random.seed(42)
-        values = pd.DataFrame(np.random.randint(low=0, high=2, size=(100, 5)),
-                              columns=['A', 'B', 'C', 'D', 'E'])
-        fit_data = values[:80]
-        predict_data = values[80:].copy()
-        self.model_connected.fit(fit_data)
-        predict_data.drop('E', axis=1, inplace=True)
-        e_prob = self.model_connected.predict_probability(predict_data)
-        np_test.assert_allclose(e_prob.values.ravel(),
-                                    np.array([0.57894737,  0.42105263,  0.57894737,  0.42105263,  0.57894737,
-                                             0.42105263,  0.5       ,  0.5       ,  0.57894737,  0.42105263,
-                                             0.5       ,  0.5       ,  0.57894737,  0.42105263,  0.57894737,
-                                             0.42105263,  0.57894737,  0.42105263,  0.5       ,  0.5       ,
-                                             0.57894737,  0.42105263,  0.57894737,  0.42105263,  0.5       ,
-                                             0.5       ,  0.57894737,  0.42105263,  0.57894737,  0.42105263,
-                                             0.5       ,  0.5       ,  0.57894737,  0.42105263,  0.5       ,
-                                             0.5       ,  0.5       ,  0.5       ,  0.5       ,  0.5       ]), atol = 0)
-        predict_data = pd.DataFrame(np.random.randint(low=0, high=2, size=(1, 5)),
-                              columns=['A', 'B', 'C', 'F', 'E'])[:]
-
-    def test_predict_probability_errors(self):
-        np.random.seed(42)
-        values = pd.DataFrame(np.random.randint(low=0, high=2, size=(2, 5)),
-                              columns=['A', 'B', 'C', 'D', 'E'])
-        fit_data = values[:1]
-        predict_data = values[1:].copy()
-        self.model_connected.fit(fit_data)
-        self.assertRaises(ValueError, self.model_connected.predict_probability, predict_data)
-        predict_data = pd.DataFrame(np.random.randint(low=0, high=2, size=(1, 5)),
-                              columns=['A', 'B', 'C', 'F', 'E'])[:]
-        self.assertRaises(ValueError, self.model_connected.predict_probability, predict_data)
+                                             1, 1, 1, 0]))
 
     def tearDown(self):
         del self.model_connected
@@ -582,7 +531,7 @@ class TestDirectedGraphCPDOperations(unittest.TestCase):
         self.graph.remove_cpds('diff', 'grade')
         self.assertListEqual(self.graph.get_cpds(), [cpd2])
 
-    def test_get_values_for_node(self):
+    def test_get_cpd_for_node(self):
         cpd1 = TabularCPD('diff', 2, values=np.random.rand(2, 1))
         cpd2 = TabularCPD('intel', 2, values=np.random.rand(2, 1))
         cpd3 = TabularCPD('grade', 2, values=np.random.rand(2, 4),
@@ -593,7 +542,7 @@ class TestDirectedGraphCPDOperations(unittest.TestCase):
         self.assertEqual(self.graph.get_cpds('intel'), cpd2)
         self.assertEqual(self.graph.get_cpds('grade'), cpd3)
 
-    def test_get_values_raises_error(self):
+    def test_get_cpd_raises_error(self):
         cpd1 = TabularCPD('diff', 2, values=np.random.rand(2, 1))
         cpd2 = TabularCPD('intel', 2, values=np.random.rand(2, 1))
         cpd3 = TabularCPD('grade', 2, values=np.random.rand(2, 4),
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_readwrite/test_BIF.py b/pgmpy/tests/test_readwrite/test_BIF.py
index a7e1f5ff..39560ff2 100644
--- a/pgmpy/tests/test_readwrite/test_BIF.py
+++ b/pgmpy/tests/test_readwrite/test_BIF.py
@@ -93,7 +93,7 @@ probability (  "family-out" ) { //1 variable(s) and 2 values
             self.assertListEqual(property_expected[variable],
                                  prop[variable])
 
-    def test_get_values(self):
+    def test_get_cpd(self):
 
         cpd_expected = {'bowel-problem': np.array([[0.01],
                                                    [0.99]]),
@@ -166,7 +166,7 @@ probability (  "family-out" ) { //1 variable(s) and 2 values
                                    [0.4, 0.95]])]
         model = self.reader.get_model()
         for cpd_index in range(0, len(cpds_expected)):
-            np_test.assert_array_equal(model.get_cpds()[cpd_index].get_values(),
+            np_test.assert_array_equal(model.get_cpds()[cpd_index].get_cpd(),
                                        cpds_expected[cpd_index])
         self.assertDictEqual(model.node, node_expected)
         self.assertDictEqual(model.edge, edge_expected)
diff --git a/pgmpy/tests/test_readwrite/test_ProbModelXML.py b/pgmpy/tests/test_readwrite/test_ProbModelXML.py
index 39a8fb40..60400afa 100644
--- a/pgmpy/tests/test_readwrite/test_ProbModelXML.py
+++ b/pgmpy/tests/test_readwrite/test_ProbModelXML.py
@@ -611,7 +611,7 @@ class TestProbModelXMLReaderString(unittest.TestCase):
                          np.array([[0.99, 0.01], [0.95, 0.05]]),
                          np.array([[1, 0, 0, 1], [0, 1, 0, 1]])]
         for cpd_index in range(0, len(cpds_expected)):
-            np_test.assert_array_equal(model.get_cpds()[cpd_index].get_values(),
+            np_test.assert_array_equal(model.get_cpds()[cpd_index].get_cpd(),
                                        cpds_expected[cpd_index])
         self.assertDictEqual(model.node, node_expected)
         self.assertDictEqual(model.edge, edge_expected)
diff --git a/pgmpy/tests/test_readwrite/test_XMLBIF.py b/pgmpy/tests/test_readwrite/test_XMLBIF.py
index b0007445..425472eb 100644
--- a/pgmpy/tests/test_readwrite/test_XMLBIF.py
+++ b/pgmpy/tests/test_readwrite/test_XMLBIF.py
@@ -164,7 +164,7 @@ class TestXMLBIFReaderMethods(unittest.TestCase):
         self.assertListEqual(sorted(self.reader.edge_list),
                              sorted(edges_expected))
 
-    def test_get_values(self):
+    def test_get_cpd(self):
         cpd_expected = {'bowel_problem': np.array([[0.01],
                                                    [0.99]]),
                         'dog_out': np.array([[0.99, 0.97, 0.9, 0.3],
@@ -240,7 +240,7 @@ class TestXMLBIFReaderMethodsFile(unittest.TestCase):
         self.assertListEqual(sorted(self.reader.edge_list),
                              sorted(edges_expected))
 
-    def test_get_values(self):
+    def test_get_cpd(self):
         cpd_expected = {'bowel_problem': np.array([[0.01],
                                                    [0.99]]),
                         'dog_out': np.array([[0.99, 0.97, 0.9, 0.3],
diff --git a/pgmpy/tests/test_readwrite/test_XMLBeliefNetwork.py b/pgmpy/tests/test_readwrite/test_XMLBeliefNetwork.py
index 7141359b..7aeedae5 100644
--- a/pgmpy/tests/test_readwrite/test_XMLBeliefNetwork.py
+++ b/pgmpy/tests/test_readwrite/test_XMLBeliefNetwork.py
@@ -229,7 +229,7 @@ class TestXBNReader(unittest.TestCase):
                          'd': np.array([[0.8, 0.2, 0.9, 0.1],
                                         [0.7, 0.3, 0.05, 0.95]])}
         for cpd in model.get_cpds():
-            np_test.assert_array_equal(cpd.get_values(), cpds_expected[cpd.variable])
+            np_test.assert_array_equal(cpd.get_cpd(), cpds_expected[cpd.variable])
         self.assertListEqual(sorted(model.edges()), sorted([('b', 'd'), ('a', 'b'), ('a', 'c'),
                                                             ('c', 'd'), ('c', 'e')]))
         self.assertDictEqual(model.node, node_expected)
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)
diff --git a/pgmpy/tests/test_utils/test_state_name.py b/pgmpy/tests/test_utils/test_state_name.py
index 871e917a..40fced54 100644
--- a/pgmpy/tests/test_utils/test_state_name.py
+++ b/pgmpy/tests/test_utils/test_state_name.py
@@ -144,7 +144,7 @@ class StateNameDecorator(unittest.TestCase):
         cpd.reduce([('diff', 'high')])
         self.assertEqual(cpd.variable, 'grade')
         self.assertEqual(cpd.variables, ['grade', 'intel'])
-        np_test.assert_array_equal(cpd.get_values(), np.array([[0.1, 0.1, 0.1],
+        np_test.assert_array_equal(cpd.get_cpd(), np.array([[0.1, 0.1, 0.1],
                                                             [0.1, 0.1, 0.1],
                                                             [0.8, 0.8, 0.8]]))
 
@@ -156,7 +156,7 @@ class StateNameDecorator(unittest.TestCase):
         cpd.reduce([('diff', 0)])
         self.assertEqual(cpd.variable, 'grade')
         self.assertEqual(cpd.variables, ['grade', 'intel'])
-        np_test.assert_array_equal(cpd.get_values(), np.array([[0.1, 0.1, 0.1],
+        np_test.assert_array_equal(cpd.get_cpd(), np.array([[0.1, 0.1, 0.1],
                                                             [0.1, 0.1, 0.1],
                                                             [0.8, 0.8, 0.8]]))
 
@@ -168,7 +168,7 @@ class StateNameDecorator(unittest.TestCase):
         cpd = cpd.reduce([('diff', 'high')], inplace=False)
         self.assertEqual(cpd.variable, 'grade')
         self.assertEqual(cpd.variables, ['grade', 'intel'])
-        np_test.assert_array_equal(cpd.get_values(), np.array([[0.1, 0.1, 0.1],
+        np_test.assert_array_equal(cpd.get_cpd(), np.array([[0.1, 0.1, 0.1],
                                                             [0.1, 0.1, 0.1],
                                                             [0.8, 0.8, 0.8]]))
 
@@ -180,7 +180,7 @@ class StateNameDecorator(unittest.TestCase):
         cpd = cpd.reduce([('diff', 0)], inplace=False)
         self.assertEqual(cpd.variable, 'grade')
         self.assertEqual(cpd.variables, ['grade', 'intel'])
-        np_test.assert_array_equal(cpd.get_values(), np.array([[0.1, 0.1, 0.1],
+        np_test.assert_array_equal(cpd.get_cpd(), np.array([[0.1, 0.1, 0.1],
                                                             [0.1, 0.1, 0.1],
                                                             [0.8, 0.8, 0.8]]))