Sampling Methods

Gibbs Sampler

class pgmpy.sampling.Sampling.GibbsSampling(model=None)[source]

Class for performing Gibbs sampling.

model: BayesianModel or MarkovModel
Model from which variables are inherited and transition probabilites computed.

set_start_state(state) sample(start_state, size) generate_sample(start_state, size)

Initialization from a BayesianModel object: >>> from pgmpy.factors.discrete import TabularCPD >>> from pgmpy.models import BayesianModel >>> intel_cpd = TabularCPD(‘intel’, 2, [[0.7], [0.3]]) >>> sat_cpd = TabularCPD(‘sat’, 2, [[0.95, 0.2], [0.05, 0.8]], evidence=[‘intel’], evidence_card=[2]) >>> student = BayesianModel() >>> 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

0 0 0 1 0 0 2 1 1

generate_sample(start_state=None, size=1)[source]

Generator version of self.sample

List of State namedtuples, representing the assignment to all variables of the model.

>>> from pgmpy.factors.discrete import DiscreteFactor
>>> from pgmpy.sampling import GibbsSampling
>>> from pgmpy.models import MarkovModel
>>> model = MarkovModel([('A', 'B'), ('C', 'B')])
>>> factor_ab = DiscreteFactor(['A', 'B'], [2, 2], [1, 2, 3, 4])
>>> factor_cb = DiscreteFactor(['C', 'B'], [2, 2], [5, 6, 7, 8])
>>> model.add_factors(factor_ab, factor_cb)
>>> gibbs = GibbsSampling(model)
>>> gen = gibbs.generate_sample(size=2)
>>> [sample for sample in gen]
[[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)]]
sample(start_state=None, size=1)[source]

Sample from the Markov Chain.

start_state: dict or array-like iterable
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.

pandas.DataFrame

>>> from pgmpy.factors import DiscreteFactor
>>> from pgmpy.inference import GibbsSampling
>>> from pgmpy.models import MarkovModel
>>> model = MarkovModel([('A', 'B'), ('C', 'B')])
>>> factor_ab = DiscreteFactor(['A', 'B'], [2, 2], [1, 2, 3, 4])
>>> 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)
   A  B  C
0  0  1  1
1  1  0  0
2  1  1  0
3  1  1  1

Bayesian Model Samplers

class pgmpy.sampling.Sampling.BayesianModelSampling(model)[source]

Class for sampling methods specific to Bayesian Models

model: instance of BayesianModel
model on which inference queries will be computed

forward_sample(size)

forward_sample(size=1)[source]

Generates sample(s) from joint distribution of the bayesian network.

size: int
size of sample to be generated
sampled: pandas.DataFrame
the generated samples
>>> from pgmpy.models.BayesianModel import BayesianModel
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianModel([('diff', 'grade'), ('intel', 'grade')])
>>> cpd_d = TabularCPD('diff', 2, [[0.6], [0.4]])
>>> cpd_i = TabularCPD('intel', 2, [[0.7], [0.3]])
>>> cpd_g = TabularCPD('grade', 3, [[0.3, 0.05, 0.9, 0.5], [0.4, 0.25,
...                0.08, 0.3], [0.3, 0.7, 0.02, 0.2]],
...                ['intel', 'diff'], [2, 2])
>>> student.add_cpds(cpd_d, cpd_i, cpd_g)
>>> inference = BayesianModelSampling(student)
>>> inference.forward_sample(2)
        diff       intel       grade
0        1           0          1
1        1           0          2
likelihood_weighted_sample(evidence=None, size=1)[source]

Generates weighted sample(s) from joint distribution of the bayesian network, that comply with the given evidence. ‘Probabilistic Graphical Model Principles and Techniques’, Koller and Friedman, Algorithm 12.2 pp 493.

evidence: list of pgmpy.factor.State namedtuples
None if no evidence
size: int
size of sample to be generated
sampled: pandas.DataFrame
the generated samples with corresponding weights
>>> from pgmpy.factors.discrete import State
>>> from pgmpy.models.BayesianModel import BayesianModel
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianModel([('diff', 'grade'), ('intel', 'grade')])
>>> cpd_d = TabularCPD('diff', 2, [[0.6], [0.4]])
>>> cpd_i = TabularCPD('intel', 2, [[0.7], [0.3]])
>>> cpd_g = TabularCPD('grade', 3, [[0.3, 0.05, 0.9, 0.5], [0.4, 0.25,
...         0.08, 0.3], [0.3, 0.7, 0.02, 0.2]],
...         ['intel', 'diff'], [2, 2])
>>> student.add_cpds(cpd_d, cpd_i, cpd_g)
>>> inference = BayesianModelSampling(student)
>>> evidence = [State('diff', 0)]
>>> inference.likelihood_weighted_sample(evidence, 2)
        intel       diff       grade  _weight
0         0          0          1        0.6
1         1          0          1        0.6
rejection_sample(evidence=None, size=1)[source]

Generates sample(s) from joint distribution of the bayesian network, given the evidence.

evidence: list of pgmpy.factor.State namedtuples
None if no evidence
size: int
size of sample to be generated
sampled: pandas.DataFrame
the generated samples
>>> from pgmpy.models.BayesianModel import BayesianModel
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.factors.discrete import State
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianModel([('diff', 'grade'), ('intel', 'grade')])
>>> cpd_d = TabularCPD('diff', 2, [[0.6], [0.4]])
>>> cpd_i = TabularCPD('intel', 2, [[0.7], [0.3]])
>>> cpd_g = TabularCPD('grade', 3, [[0.3, 0.05, 0.9, 0.5], [0.4, 0.25,
...                0.08, 0.3], [0.3, 0.7, 0.02, 0.2]],
...                ['intel', 'diff'], [2, 2])
>>> student.add_cpds(cpd_d, cpd_i, cpd_g)
>>> inference = BayesianModelSampling(student)
>>> evidence = [State(var='diff', state=0)]
>>> inference.rejection_sample(evidence, 2)
        intel       diff       grade
0         0          0          1
1         0          0          1

Hamiltonian Monte Carlo

A collection of methods for sampling from continuous models in pgmpy

class pgmpy.sampling.HMC.HamiltonianMC(model, grad_log_pdf, simulate_dynamics=<class 'pgmpy.sampling.base.LeapFrog'>)[source]

Class for performing sampling using simple Hamiltonian Monte Carlo

model: An instance pgmpy.models
Model from which sampling has to be done
grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
A class to find log and gradient of log distribution for a given assignment If None, then will use model.get_gradient_log_pdf
simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
A class to propose future values of momentum and position in time by simulating Hamiltonian Dynamics

sample() generate_sample()

>>> from pgmpy.sampling import HamiltonianMC as HMC, LeapFrog, GradLogPDFGaussian
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([-3, 4])
>>> covariance = np.array([[3, 0.7], [0.7, 5]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
>>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
...                          trajectory_length=2, stepsize=0.4)
>>> samples
rec.array([(1.0, 1.0), (-3.1861687131079086, 3.7940994520145654),
 (-1.6920542547310844, 6.347410703806017), ...,
 (-1.8093621120575312, 5.940860883943261),
 (0.3933248026088032, 6.3853098838119235),
 (-0.8654072934719572, 6.023803629334816)],
          dtype=[('x', '<f8'), ('y', '<f8')])
>>> samples = np.array([samples[var_name] for var_name in model.variables])
>>> np.cov(samples)
array([[ 3.0352818 ,  0.71379304],
       [ 0.71379304,  4.91776713]])
>>> sampler.accepted_proposals
9932.0
>>> sampler.acceptance_rate
0.9932

R.Neal. Handbook of Markov Chain Monte Carlo, chapter 5: MCMC Using Hamiltonian Dynamics. CRC Press, 2011.

generate_sample(initial_pos, num_samples, trajectory_length, stepsize=None)[source]

Method returns a generator type object whose each iteration yields a sample using Hamiltonian Monte Carlo

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_samples: int
Number of samples to be generated
trajectory_length: int or float
Target trajectory length, stepsize * number of steps(L), where L is the number of steps taken per HMC iteration, and stepsize is step size for splitting time method.
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

genrator: yielding a 1d numpy.array type object for a sample

>>> 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])
>>> covariance = np.array([[3, 0.4], [0.4, 3]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = HMC(model=model, grad_log_pdf=GLPG)
>>> gen_samples = sampler.generate_sample(np.array([-1, 1]), num_samples = 10000,
...                                       trajectory_length=2, stepsize=0.25)
>>> samples_array = np.array([sample for sample in gen_samples])
>>> samples_array
array([[ 0.1467264 ,  0.27143857],
       [ 4.0371448 ,  0.15871274],
       [ 3.24656208, -1.03742621],
       ...,
       [ 6.45975905,  1.97941306],
       [ 4.89007171,  0.15413156],
       [ 5.9528083 ,  1.92983158]])
>>> np.cov(samples_array.T)
array([[ 2.95692642,  0.4379419 ],
       [ 0.4379419 ,  3.00939434]])
>>> sampler.acceptance_rate
0.9969
sample(initial_pos, num_samples, trajectory_length, stepsize=None)[source]

Method to return samples using Hamiltonian Monte Carlo

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_samples: int
Number of samples to be generated
trajectory_length: int or float
Target trajectory length, stepsize * number of steps(L), where L is the number of steps taken per HMC iteration, and stepsize is step size for splitting time method.
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

Returns 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)

>>> # 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
>>> mean = np.array([1, -1])
>>> covariance = np.array([[1, 0.2], [0.2, 1]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=ModifiedEuler)
>>> samples = sampler.sample(np.array([1, 1]), num_samples = 5,
...                          trajectory_length=6, stepsize=0.25)
>>> samples
               x              y
0   1.000000e+00   1.000000e+00
1   1.592133e+00   1.152911e+00
2   1.608700e+00   1.315349e+00
3   1.608700e+00   1.315349e+00
4   6.843856e-01   6.237043e-01
>>> mean = np.array([4, 1, -1])
>>> covariance = np.array([[1, 0.7 , 0.8], [0.7, 1, 0.2], [0.8, 0.2, 1]])
>>> model = JGD(['x', 'y', 'z'], mean, covariance)
>>> sampler = HMC(model=model, grad_log_pdf=GLPG)
>>> samples = sampler.sample(np.array([1, 1]), num_samples = 10000,
...                          trajectory_length=6, stepsize=0.25)
>>> np.cov(samples.values.T)
array([[ 1.00795398,  0.71384233,  0.79802097],
       [ 0.71384233,  1.00633524,  0.21313767],
       [ 0.79802097,  0.21313767,  0.98519017]])
class pgmpy.sampling.HMC.HamiltonianMCDA(model, grad_log_pdf=None, simulate_dynamics=<class 'pgmpy.sampling.base.LeapFrog'>, delta=0.65)[source]

Class for performing sampling in Continuous model using Hamiltonian Monte Carlo with dual averaging for adaptaion of parameter stepsize.

model: An instance pgmpy.models
Model from which sampling has to be done
grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
Class to compute the log and gradient log of distribution
simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
Class to propose future states of position and momentum in time by simulating HamiltonianDynamics
delta: float (in between 0 and 1), defaults to 0.65
The target HMC acceptance probability

sample() generate_sample()

>>> 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)
>>> samples = sampler.sample(np.array([0, 0, 0]), num_adapt=10000, num_samples = 10000, trajectory_length=7)
>>> samples_array = np.array([samples[var_name] for var_name in model.variables])
>>> np.cov(samples_array)
array([[ 1.83023816,  0.40449162,  0.51200707],
       [ 0.40449162,  2.85863596,  0.76747343],
       [ 0.51200707,  0.76747343,  3.87020982]])
>>> sampler.acceptance_rate
0.9929

Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15 (2014) 1351-1381 Algorithm 5 : Hamiltonian Monte Carlo with dual averaging

generate_sample(initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None)[source]

Method returns a generator type object whose each iteration yields a sample using Hamiltonian Monte Carlo

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_adapt: int
The number of interations to run the adaptation of stepsize
num_samples: int
Number of samples to be generated
trajectory_length: int or float
Target trajectory length, stepsize * number of steps(L), where L is the number of steps taken to propose new values of position and momentum per HMC iteration and stepsize is step size.
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

genrator: yielding a numpy.array type object for a sample

>>> from pgmpy.sampling import HamiltonianMCDA as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([1, 1])
>>> covariance = np.array([[1, 0.7], [0.7, 3]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = HMCda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
>>> gen_samples = sampler.generate_sample(np.array([1, 1]), num_adapt=10000,
...                                       num_samples = 10000, trajectory_length=2, stepsize=None)
>>> samples_array = np.array([sample for sample in gen_samples])
>>> np.cov(samples_array.T)
array([[ 0.98432155,  0.69517394],
       [ 0.69517394,  2.95449533]])
sample(initial_pos, num_adapt, num_samples, trajectory_length, stepsize=None)[source]

Method to return samples using Hamiltonian Monte Carlo

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_adapt: int
The number of interations to run the adaptation of stepsize
num_samples: int
Number of samples to be generated
trajectory_length: int or float
Target trajectory length, stepsize * number of steps(L), where L is the number of steps taken per HMC iteration, and stepsize is step size for splitting time method.
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

Returns 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)

>>> from pgmpy.sampling import HamiltonianMCDA as HMCda, GradLogPDFGaussian as GLPG, LeapFrog
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([1, 1])
>>> covariance = np.array([[1, 0.7], [0.7, 3]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = HMCda(model=model, grad_log_pdf=GLPG, simulate_dynamics=LeapFrog)
>>> samples = sampler.sample(np.array([1, 1]), num_adapt=10000,
...                          num_samples = 10000, trajectory_length=2, stepsize=None)
>>> samples_array = np.array([samples[var_name] for var_name in model.variables])
>>> np.cov(samples_array)
array([[ 0.98432155,  0.66517394],
       [ 0.66517394,  2.95449533]])

No U-Turn Sampler

class pgmpy.sampling.NUTS.NoUTurnSampler(model, grad_log_pdf, simulate_dynamics=<class 'pgmpy.sampling.base.LeapFrog'>)[source]

Class for performing sampling in Continuous model using No U Turn Sampler (a variant of Hamiltonian Monte Carlo)

model: An instance pgmpy.models
Model from which sampling has to be done
grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
Class to compute the log and gradient log of distribution
simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
Class to propose future states of position and momentum in time by simulating HamiltonianDynamics

sample() generate_sample()

>>> from pgmpy.sampling import NoUTurnSampler as NUTS, LeapFrog, GradLogPDFGaussian
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([1, 2, 3])
>>> covariance = np.array([[4, 0.1, 0.2], [0.1, 1, 0.3], [0.2, 0.3, 8]])
>>> model = JGD(['x', 'y', 'z'], mean, covariance)
>>> sampler = NUTS(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
>>> samples = sampler.sample(initial_pos=np.array([0.1, 0.9, 0.3]), num_samples=20000,stepsize=0.4)
>>> samples
rec.array([(0.1, 0.9, 0.3),
 (-0.27303886844752756, 0.5028580705249155, 0.2895768065049909),
 (1.7139810571103862, 2.809135711303245, 5.690811523613858), ...,
 (-0.7742669710786649, 2.092867703984895, 6.139480724333439),
 (1.3916152816323692, 1.394952482021687, 3.446906546649354),
 (-0.2726336476939125, 2.6230854954595357, 2.923948403903159)],
          dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])

Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15 (2014) 1351-1381 Algorithm 3 : Efficient No-U-Turn Sampler

generate_sample(initial_pos, num_samples, stepsize=None)[source]

Returns a generator type object whose each iteration yields a sample

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_samples: int
Number of samples to be generated
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

generator: yielding a numpy.array type object for a sample

>>> from pgmpy.sampling import NoUTurnSampler as NUTS, GradLogPDFGaussian
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([11, -6])
>>> covariance = np.array([[0.7, 0.2], [0.2, 14]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = NUTS(model=model, grad_log_pdf=GradLogPDFGaussian)
>>> samples = sampler.generate_sample(initial_pos=np.array([1, 1]), num_samples=10, stepsize=0.4)
>>> samples = np.array([sample for sample in samples])
>>> samples
array([[ 10.26357538,   0.10062725],
       [ 12.70600336,   0.63392499],
       [ 10.95523217,  -0.62079273],
       [ 10.66263031,  -4.08135962],
       [ 10.59255762,  -8.48085076],
       [  9.99860242,  -9.47096032],
       [ 10.5733564 ,  -9.83504745],
       [ 11.51302059,  -9.49919523],
       [ 11.31892143,  -8.5873259 ],
       [ 11.29008667,  -0.43809674]])
sample(initial_pos, num_samples, stepsize=None)[source]

Method to return samples using No U Turn Sampler

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_samples: int
Number of samples to be generated
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

Returns 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)

>>> # 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)
>>> samples
          x         y         z
0  1.000000  1.000000  1.000000
1  1.760756  0.271543 -0.613309
2  1.883387  0.990745 -0.611720
3  0.980812  0.340336 -0.916283
4  0.781338  0.647220 -0.948640
5  0.040308 -1.391406  0.412201
6  1.179549 -1.450552  1.105216
7  1.100320 -1.313926  1.207815
8  1.484520 -1.349247  0.768599
9  0.934942 -1.894589  0.471772
class pgmpy.sampling.NUTS.NoUTurnSamplerDA(model, grad_log_pdf, simulate_dynamics=<class 'pgmpy.sampling.base.LeapFrog'>, delta=0.65)[source]

Class for performing sampling in Continuous model using No U Turn sampler with dual averaging for adaptation of parameter stepsize.

model: An instance pgmpy.models
Model from which sampling has to be done
grad_log_pdf: A subclass of pgmpy.inference.continuous.GradientLogPDF
Class to compute the log and gradient log of distribution
simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
Class to propose future states of position and momentum in time by simulating HamiltonianDynamics
delta: float (in between 0 and 1), defaults to 0.65
The target HMC acceptance probability

sample() generate_sample()

>>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([-1, 12, -3])
>>> covariance = np.array([[-2, 7, 2], [7, 14, 4], [2, 4, -1]])
>>> model = JGD(['x', 'v', 't'], mean, covariance)
>>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian)
>>> samples = sampler.sample(initial_pos=np.array([0, 0, 0]), num_adapt=10, num_samples=10, stepsize=0.25)
>>> samples
rec.array([(0.0, 0.0, 0.0),
 (0.06100992691638076, -0.17118088764170125, 0.14048470935160887),
 (0.06100992691638076, -0.17118088764170125, 0.14048470935160887),
 (-0.7451883138013118, 1.7975387358691155, 2.3090698721374436),
 (-0.6207457594500309, 1.4611049498441024, 2.5890867012835574),
 (0.24043604780911487, 1.8660976216530618, 3.2508715592645347),
 (0.21509819341468212, 2.157760225367607, 3.5749582768731476),
 (0.20699150582681913, 2.0605044285377305, 3.8588980251618135),
 (0.20699150582681913, 2.0605044285377305, 3.8588980251618135),
 (0.085332419611991, 1.7556171374575567, 4.49985082288814)],
          dtype=[('x', '<f8'), ('v', '<f8'), ('t', '<f8')])

Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15 (2014) 1351-1381 Algorithm 6 : No-U-Turn Sampler with Dual Averaging

generate_sample(initial_pos, num_adapt, num_samples, stepsize=None)[source]

Returns a generator type object whose each iteration yields a sample

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_adapt: int
The number of interations to run the adaptation of stepsize
num_samples: int
Number of samples to be generated
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

genrator: yielding a numpy.array type object for a sample

>>> from pgmpy.sampling import NoUTurnSamplerDA as NUTSda, GradLogPDFGaussian
>>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
>>> import numpy as np
>>> mean = np.array([1, -100])
>>> covariance = np.array([[-12, 45], [45, -10]])
>>> model = JGD(['a', 'b'], mean, covariance)
>>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
>>> samples = sampler.generate_sample(initial_pos=np.array([12, -4]), num_adapt=10,
...                                   num_samples=10, stepsize=0.1)
>>> samples
<generator object NoUTurnSamplerDA.generate_sample at 0x7f4fed46a4c0>
>>> samples_array = np.array([sample for sample in samples])
>>> samples_array
array([[ 11.89963386,  -4.06572636],
       [ 10.3453755 ,  -7.5700289 ],
       [-26.56899659, -15.3920684 ],
       [-29.97143077, -12.0801625 ],
       [-29.97143077, -12.0801625 ],
       [-33.07960829,  -8.90440347],
       [-55.28263496, -17.31718524],
       [-55.28263496, -17.31718524],
       [-56.63440044, -16.03309364],
       [-63.880094  , -19.19981944]])
sample(initial_pos, num_adapt, num_samples, stepsize=None)[source]

Returns samples using No U Turn Sampler with dual averaging

initial_pos: A 1d array like object
Vector representing values of parameter position, the starting state in markov chain.
num_adapt: int
The number of interations to run the adaptation of stepsize
num_samples: int
Number of samples to be generated
stepsize: float , defaults to None
The stepsize for proposing new values of position and momentum in simulate_dynamics If None, then will be choosen suitably

Returns 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)

>>> # 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
>>> mean = np.array([10, -13])
>>> covariance = np.array([[16, -3], [-3, 13]])
>>> model = JGD(['x', 'y'], mean, covariance)
>>> sampler = NUTSda(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
>>> samples = sampler.sample(initial_pos=np.array([12, -4]), num_adapt=10, num_samples=10, stepsize=0.1)
>>> samples
           x          y
0  12.000000  -4.000000
1  11.864821  -3.696109
2  10.546986  -4.892169
3   8.526596 -21.555793
4   8.526596 -21.555793
5  11.343194  -6.353789
6  -1.583269 -12.802931
7  12.411957 -11.704859
8  13.253336 -20.169492
9  11.295901  -7.665058