Repository: pgmpy/pgmpy
Pull request #788 (2016-12-17 08:42:28)

Markov Chain accepts matrix PR#748


Patch:
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/MarkovChain.py b/pgmpy/models/MarkovChain.py
index e276b456..846ff993 100644
--- a/pgmpy/models/MarkovChain.py
+++ b/pgmpy/models/MarkovChain.py
@@ -177,8 +177,11 @@ class MarkovChain(object):
         -----------
         variable: any hashable python object
             must be an existing variable of the model.
-        transition_model: dict
-            representing valid transition probabilities defined for every possible state of the variable.
+        transition_model: dict or 2d array
+            dict representing valid transition probabilities defined for every possible state of the variable.
+            array represent a square matrix where every row sums to 1,
+            array[i,j] indicates the transition probalities from State i to State j
+
 
         Examples:
         ---------
@@ -186,11 +189,27 @@ class MarkovChain(object):
         >>> model = MC()
         >>> model.add_variable('grade', 3)
         >>> grade_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}}
+        >>> grade_tm_matrix = np.array([[0.1, 0.5, 0.4], [0.2, 0.2, 0.6], [0.7, 0.15, 0.15]])
         >>> model.add_transition_model('grade', grade_tm)
+        >>> model.add_transition_model('grade', grade_tm_matrix)
         """
+        if isinstance(transition_model, list):
+            transition_model = np.array(transition_model)
+
         # check if the transition model is valid
         if not isinstance(transition_model, dict):
-            raise ValueError('Transition model must be a dict.')
+            if not isinstance(transition_model, np.ndarray):
+                raise ValueError('Transition model must be a dict or numpy array')
+            elif len(transition_model.shape) != 2:
+                raise ValueError('Transition model must be 2d array.given {t}'.format(t=transition_model.shape))
+            elif transition_model.shape[0] != transition_model.shape[1]:
+                raise ValueError('Dimension mismatch {d1}!={d2}'.format(d1=transition_model.shape[0],
+                                 d2=transition_model.shape[1]))
+            else:
+                # convert the matrix to dict
+                size = transition_model.shape[0]
+                transition_model = dict((i, dict((j, float(transition_model[i][j]))
+                                         for j in range(0, size))) for i in range(0, size))
 
         exp_states = set(range(self.cardinalities[variable]))
         tm_states = set(transition_model.keys())
@@ -284,8 +303,9 @@ class MarkovChain(object):
         Examples:
         ---------
         >>> 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: 1}}
+        >>> 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}}
         >>> 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)
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


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 64528f69..fb5e6da1 100644
--- a/pgmpy/tests/test_models/test_MarkovChain.py
+++ b/pgmpy/tests/test_models/test_MarkovChain.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 import sys
 import unittest
-
+import numpy as np
 from pandas import DataFrame
 from mock import patch, call
 
@@ -16,8 +16,11 @@ class TestMarkovChain(unittest.TestCase):
         self.card = [3, 2, 3]
         self.cardinalities = {'intel': 3, 'diff': 2, 'grade': 3}
         self.intel_tm = {0: {0: 0.1, 1: 0.25, 2: 0.65}, 1: {0: 0.5, 1: 0.3, 2: 0.2}, 2: {0: 0.3, 1: 0.3, 2: 0.4}}
+        self.intel_tm_matrix = np.array([[0.1, 0.25, 0.65], [0.5, 0.3, 0.2], [0.3, 0.3, 0.4]])
         self.diff_tm = {0: {0: 0.3, 1: 0.7}, 1: {0: 0.75, 1: 0.25}}
+        self.diff_tm_matrix = np.array([[0.3, 0.7], [0.75, 0.25]])
         self.grade_tm = {0: {0: 0.4, 1: 0.2, 2: 0.4}, 1: {0: 0.9, 1: 0.05, 2: 0.05}, 2: {0: 0.1, 1: 0.4, 2: 0.5}}
+        self.grade_tm_matrix = [[0.4, 0.2, 0.4], [0.9, 0.05, 0.05], [0.1, 0.4, 0.5]]
         self.start_state = [State('intel', 0), State('diff', 1), State('grade', 2)]
         self.model = MC()
 
@@ -149,8 +152,10 @@ class TestMarkovChain(unittest.TestCase):
 
     def test_add_transition_model_bad_type(self):
         model = MC()
-        # if transition_model is not a dict
+        grade_tm_matrix_bad = [[0.1, 0.5, 0.4], [0.2, 0.2, 0.6], 'abc']
+        # if transition_model is not a dict or np.array
         self.assertRaises(ValueError, model.add_transition_model, 'var', 123)
+        self.assertRaises(ValueError, model.add_transition_model, 'var', grade_tm_matrix_bad)
 
     def test_add_transition_model_bad_states(self):
         model = MC(['var'], [2])
@@ -181,6 +186,21 @@ class TestMarkovChain(unittest.TestCase):
         model.add_transition_model('var', transition_model)
         self.assertDictEqual(model.transition_models['var'], transition_model)
 
+    def test_transition_model_bad_matrix_dimension(self):
+        model = MC(['var'], [2])
+        transition_model = np.array([0.3, 0.7])
+        # check for square dimension of the matrix
+        self.assertRaises(ValueError, model.add_transition_model, 'var', transition_model)
+        transition_model = np.array([[0.3, 0.6, 0.1], [0.3, 0.3, 0.4]])
+        self.assertRaises(ValueError, model.add_transition_model, 'var', transition_model)
+
+    def test_transition_model_dict_to_matrix(self):
+        model = MC(['var'], [2])
+        transition_model = {0: {0: 0.3, 1: 0.7}, 1: {0: 0.5, 1: 0.5}}
+        transition_model_matrix = np.array([[0.3, 0.7], [0.5, 0.5]])
+        model.add_transition_model('var', transition_model_matrix)
+        self.assertDictEqual(model.transition_models['var'], transition_model)
+
     def test_sample(self):
         model = MC(['a', 'b'], [2, 2])
         model.transition_models['a'] = {0: {0: 0.1, 1: 0.9}, 1: {0: 0.2, 1: 0.8}}
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)