import itertools
from collections import namedtuple
import networkx as nx
import numpy as np
import pandas as pd
import torch
from joblib import Parallel, delayed
from tqdm.auto import tqdm
from pgmpy import config
from pgmpy.factors import factor_product
from pgmpy.models import BayesianNetwork, MarkovChain, MarkovNetwork
from pgmpy.sampling import BayesianModelInference, _return_samples
from pgmpy.utils.mathext import sample_discrete, sample_discrete_maps
State = namedtuple("State", ["var", "state"])
[docs]
class BayesianModelSampling(BayesianModelInference):
"""
Class for sampling methods specific to Bayesian Models
Parameters
----------
model: instance of BayesianNetwork
model on which inference queries will be computed
"""
def __init__(self, model):
super(BayesianModelSampling, self).__init__(model)
[docs]
def forward_sample(
self,
size=1,
include_latents=False,
seed=None,
show_progress=True,
partial_samples=None,
n_jobs=-1,
):
"""
Generates sample(s) from joint distribution of the Bayesian Network.
Parameters
----------
size: int
size of sample to be generated
include_latents: boolean
Whether to include the latent variable values in the generated samples.
seed: int (default: None)
If a value is provided, sets the seed for numpy.random.
show_progress: boolean
Whether to show a progress bar of samples getting generated.
partial_samples: pandas.DataFrame
A pandas dataframe specifying samples on some of the variables in the model. If
specified, the sampling procedure uses these sample values, instead of generating them.
n_jobs: int (default: -1)
The number of CPU cores to use. Default uses all cores.
Returns
-------
sampled: pandas.DataFrame
The generated samples
Examples
--------
>>> from pgmpy.models import BayesianNetwork
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianNetwork([('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(size=2)
rec.array([(0, 0, 1), (1, 0, 2)], dtype=
[('diff', '<i8'), ('intel', '<i8'), ('grade', '<i8')])
"""
sampled = pd.DataFrame(columns=list(self.model.nodes()))
if show_progress and config.SHOW_PROGRESS:
pbar = tqdm(self.topological_order)
else:
pbar = self.topological_order
if seed is not None:
np.random.seed(seed)
for node in pbar:
if show_progress and config.SHOW_PROGRESS:
pbar.set_description(f"Generating for node: {node}")
# If values specified in partial_samples, use them. Else generate the values.
if (partial_samples is not None) and (node in partial_samples.columns):
sampled[node] = partial_samples.loc[:, node].values
else:
cpd = self.model.get_cpds(node)
states = range(self.cardinality[node])
evidence = cpd.variables[1:]
if evidence:
evidence_values = np.vstack([sampled[i] for i in evidence])
unique, inverse = np.unique(
evidence_values.T, axis=0, return_inverse=True
)
unique = [tuple(u) for u in unique]
state_to_index, index_to_weight = self.pre_compute_reduce_maps(
variable=node, evidence=evidence, state_combinations=unique
)
if config.get_backend() == "numpy":
weight_index = np.array([state_to_index[u] for u in unique])[
inverse
]
else:
weight_index = torch.Tensor(
[state_to_index[u] for u in unique]
)[inverse]
sampled[node] = sample_discrete_maps(
states, weight_index, index_to_weight, size
)
else:
weights = cpd.values
sampled[node] = sample_discrete(states, weights, size)
samples_df = _return_samples(sampled, self.state_names_map)
if not include_latents:
samples_df.drop(self.model.latents, axis=1, inplace=True)
return samples_df
[docs]
def rejection_sample(
self,
evidence=[],
size=1,
include_latents=False,
seed=None,
show_progress=True,
partial_samples=None,
):
"""
Generates sample(s) from joint distribution of the Bayesian Network,
given the evidence.
Parameters
----------
evidence: list of `pgmpy.factor.State` namedtuples
None if no evidence
size: int
size of sample to be generated
include_latents: boolean
Whether to include the latent variable values in the generated samples.
seed: int (default: None)
If a value is provided, sets the seed for numpy.random.
show_progress: boolean
Whether to show a progress bar of samples getting generated.
partial_samples: pandas.DataFrame
A pandas dataframe specifying samples on some of the variables in the model. If
specified, the sampling procedure uses these sample values, instead of generating them.
Returns
-------
sampled: pandas.DataFrame
The generated samples
Examples
--------
>>> from pgmpy.models import BayesianNetwork
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.factors.discrete import State
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianNetwork([('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=evidence, size=2, return_type='dataframe')
intel diff grade
0 0 0 1
1 0 0 1
"""
if seed is not None:
np.random.seed(seed)
# If no evidence is given, it is equivalent to forward sampling.
if len(evidence) == 0:
return self.forward_sample(size=size, include_latents=include_latents)
# Setup array to be returned
sampled = pd.DataFrame()
prob = 1
i = 0
# Do the sampling by generating samples from forward sampling and rejecting the
# samples which do not match our evidence. Keep doing until we have enough
# samples.
if show_progress and config.SHOW_PROGRESS:
pbar = tqdm(total=size)
while i < size:
_size = int(((size - i) / prob) * 1.5)
# If partial_samples is specified, can only generate < partial_samples.shape[0] number of samples
# at a time. For simplicity, just generate the same size as partial_samples.shape[0].
if partial_samples is not None:
_size = partial_samples.shape[0]
_sampled = self.forward_sample(
size=_size,
include_latents=True,
show_progress=False,
partial_samples=partial_samples,
)
for var, state in evidence:
_sampled = _sampled[_sampled[var] == state]
prob = max(len(_sampled) / _size, 0.01)
sampled = pd.concat([sampled, _sampled], axis=0, join="outer").iloc[
:size, :
]
i += _sampled.shape[0]
if show_progress and config.SHOW_PROGRESS:
# Update at maximum to `size`
comp = _sampled.shape[0] if i < size else size - (i - _sampled.shape[0])
pbar.update(comp)
if show_progress and config.SHOW_PROGRESS:
pbar.close()
sampled = sampled.reset_index(drop=True)
if not include_latents:
sampled.drop(self.model.latents, axis=1, inplace=True)
return sampled
[docs]
def likelihood_weighted_sample(
self,
evidence=[],
size=1,
include_latents=False,
seed=None,
show_progress=True,
n_jobs=-1,
):
"""
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.
Parameters
----------
evidence: list of `pgmpy.factor.State` namedtuples
None if no evidence
size: int
size of sample to be generated
include_latents: boolean
Whether to include the latent variable values in the generated samples.
seed: int (default: None)
If a value is provided, sets the seed for numpy.random.
show_progress: boolean
Whether to show a progress bar of samples getting generated.
n_jobs: int (default: -1)
The number of CPU cores to use. Default uses all cores.
Returns
-------
sampled: A pandas.DataFrame
The generated samples with corresponding weights
Examples
--------
>>> from pgmpy.factors.discrete import State
>>> from pgmpy.models import BayesianNetwork
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.sampling import BayesianModelSampling
>>> student = BayesianNetwork([('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=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')])
"""
if seed is not None:
np.random.seed(seed)
# Convert evidence state names to number
evidence = [
(var, self.model.get_cpds(var).get_state_no(var, state))
for var, state in evidence
]
# Prepare the return dataframe
sampled = pd.DataFrame(columns=list(self.model.nodes()))
sampled["_weight"] = np.ones(size)
evidence_dict = dict(evidence)
if show_progress and config.SHOW_PROGRESS:
pbar = tqdm(self.topological_order)
else:
pbar = self.topological_order
# Do the sampling
for node in pbar:
if show_progress and config.SHOW_PROGRESS:
pbar.set_description(f"Generating for node: {node}")
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])
unique, inverse = np.unique(
evidence_values.T, axis=0, return_inverse=True
)
unique = [tuple(u) for u in unique]
state_to_index, index_to_weight = self.pre_compute_reduce_maps(
variable=node, evidence=evidence, state_combinations=unique
)
weight_index = np.array([state_to_index[tuple(u)] for u in unique])[
inverse
]
if node in evidence_dict:
evidence_value = evidence_dict[node]
sampled[node] = evidence_value
sampled.loc[:, "_weight"] *= np.array(
list(
map(
lambda i: index_to_weight[weight_index[i]][
evidence_value
],
range(size),
)
)
)
else:
sampled[node] = sample_discrete_maps(
states, weight_index, index_to_weight, size
)
else:
if node in evidence_dict:
sampled[node] = evidence_dict[node]
sampled.loc[:, "_weight"] *= np.array(
list(
map(lambda _: cpd.values[evidence_dict[node]], range(size))
)
)
else:
sampled[node] = sample_discrete(states, cpd.values, size)
# Postprocess the samples: Change state numbers to names, remove latents.
samples_df = _return_samples(sampled, self.state_names_map)
if not include_latents:
samples_df.drop(self.model.latents, axis=1, inplace=True)
return samples_df
[docs]
class GibbsSampling(MarkovChain):
"""
Class for performing Gibbs sampling.
Parameters
----------
model: BayesianNetwork or MarkovNetwork
Model from which variables are inherited and transition probabilities computed.
Examples
--------
Initialization from a BayesianNetwork object:
>>> from pgmpy.factors.discrete import TabularCPD
>>> from pgmpy.models import BayesianNetwork
>>> 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 = BayesianNetwork()
>>> student.add_nodes_from(['intel', 'sat'])
>>> student.add_edge('intel', 'sat')
>>> student.add_cpds(intel_cpd, sat_cpd)
>>> from pgmpy.sampling import GibbsSampling
>>> gibbs_chain = GibbsSampling(student)
>>> gibbs_chain.sample(size=3)
intel sat
0 0 0
1 0 0
2 1 1
"""
def __init__(self, model=None):
super(GibbsSampling, self).__init__()
if isinstance(model, BayesianNetwork):
self._get_kernel_from_bayesian_model(model)
elif isinstance(model, MarkovNetwork):
self._get_kernel_from_markov_model(model)
def _get_kernel_from_bayesian_model(self, model):
"""
Computes the Gibbs transition models from a Bayesian Network.
'Probabilistic Graphical Model Principles and Techniques', Koller and
Friedman, Section 12.3.3 pp 512-513.
Parameters
----------
model: BayesianNetwork
The model from which probabilities will be computed.
"""
self.variables = np.array(model.nodes())
self.latents = model.latents
self.cardinalities = {
var: model.get_cpds(var).variable_card for var in self.variables
}
for var in self.variables:
other_vars = [v for v in self.variables if var != v]
other_cards = [self.cardinalities[v] for v in other_vars]
kernel = {}
factors = [cpd.to_factor() for cpd in model.cpds if var in cpd.scope()]
factor = factor_product(*factors)
scope = set(factor.scope())
for tup in itertools.product(*[range(card) for card in other_cards]):
states = [State(v, s) for v, s in zip(other_vars, tup) if v in scope]
reduced_factor = factor.reduce(states, inplace=False)
kernel[tup] = reduced_factor.values / sum(reduced_factor.values)
self.transition_models[var] = kernel
def _get_kernel_from_markov_model(self, model):
"""
Computes the Gibbs transition models from a Markov Network.
'Probabilistic Graphical Model Principles and Techniques', Koller and
Friedman, Section 12.3.3 pp 512-513.
Parameters
----------
model: MarkovNetwork
The model from which probabilities will be computed.
"""
self.variables = np.array(model.nodes())
self.latents = model.latents
factors_dict = {var: [] for var in self.variables}
for factor in model.get_factors():
for var in factor.scope():
factors_dict[var].append(factor)
# Take factor product
factors_dict = {
var: factor_product(*factors) if len(factors) > 1 else factors[0]
for var, factors in factors_dict.items()
}
self.cardinalities = {
var: factors_dict[var].get_cardinality([var])[var] for var in self.variables
}
for var in self.variables:
other_vars = [v for v in self.variables if var != v]
other_cards = [self.cardinalities[v] for v in other_vars]
kernel = {}
factor = factors_dict[var]
scope = set(factor.scope())
for tup in itertools.product(*[range(card) for card in other_cards]):
states = [
State(first_var, s)
for first_var, s in zip(other_vars, tup)
if first_var in scope
]
reduced_factor = factor.reduce(states, inplace=False)
kernel[tup] = reduced_factor.values / sum(reduced_factor.values)
self.transition_models[var] = kernel
[docs]
def sample(self, start_state=None, size=1, seed=None, include_latents=False):
"""
Sample from the Markov Chain.
Parameters
----------
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.
seed: int (default: None)
If a value is provided, sets the seed for numpy.random.
include_latents: boolean
Whether to include the latent variable values in the generated samples.
Returns
-------
sampled: pandas.DataFrame
The generated samples
Examples
--------
>>> from pgmpy.factors.discrete import DiscreteFactor
>>> from pgmpy.sampling import GibbsSampling
>>> from pgmpy.models import MarkovNetwork
>>> model = MarkovNetwork([('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, return_tupe='dataframe')
A B C
0 0 1 1
1 1 0 0
2 1 1 0
3 1 1 1
"""
if start_state is None and self.state is None:
self.state = self.random_state()
elif start_state is not None:
self.set_start_state(start_state)
if seed is not None:
np.random.seed(seed)
types = [(str(var_name), "int") for var_name in self.variables]
sampled = np.zeros(size, dtype=types).view(np.recarray)
sampled[0] = tuple(st for var, st in self.state)
for i in tqdm(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] = tuple(st for var, st in self.state)
samples_df = _return_samples(sampled)
if not include_latents:
samples_df.drop(self.latents, axis=1, inplace=True)
return samples_df
[docs]
def generate_sample(
self, start_state=None, size=1, include_latents=False, seed=None
):
"""
Generator version of self.sample
Returns
-------
List of State namedtuples, representing the assignment to all variables of the model.
Examples
--------
>>> from pgmpy.factors.discrete import DiscreteFactor
>>> from pgmpy.sampling import GibbsSampling
>>> from pgmpy.models import MarkovNetwork
>>> model = MarkovNetwork([('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)]]
"""
if seed is not None:
np.random.seed(seed)
if start_state is None and self.state is None:
self.state = self.random_state()
elif start_state is not None:
self.set_start_state(start_state)
for i in range(size):
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)
if include_latents:
yield self.state[:]
else:
yield [s for s in self.state if i not in self.latents]
