def get_distribution(samples, variables=None):
    """
    For marginal distribution, P(A): get_distribution(samples, variables=['A'])
    For joint distribution, P(A, B): get_distribution(samples, variables=['A', 'B'])
    """
    if variables is None:
        raise ValueError("variables must be specified")

    return samples.groupby(variables).size() / samples.shape[0]

# Define a model
from pgmpy.factors.discrete import TabularCPD

from pgmpy.utils import get_example_model

alarm = get_example_model("alarm")

Normal Bayesian Network (no time variation)

Normal simulation

samples = alarm.simulate(n_samples=int(1e4))
get_distribution(samples, variables=["HISTORY"])  # P(HISTORY)

HISTORY
FALSE    0.9452
TRUE     0.0548
dtype: float64

Simulation with some evidence

samples = alarm.simulate(n_samples=int(1e4), evidence={"CVP": "NORMAL", "HR": "HIGH"})
get_distribution(samples, variables=["HISTORY"])  # P(HISTORY)

HISTORY
FALSE    0.9847
TRUE     0.0153
dtype: float64

Simulation with soft/virtual evidence

soft_evidence = [
    TabularCPD(
        "CVP", 3, [[0.2], [0.3], [0.5]], state_names={"CVP": ["NORMAL", "LOW", "HIGH"]}
    )
]
samples = alarm.simulate(n_samples=int(1e4), virtual_evidence=soft_evidence)
get_distribution(samples, variables=["HISTORY"])  # P(HISTORY)

HISTORY
FALSE    0.9609
TRUE     0.0391
dtype: float64

Simulation with intervention

samples = alarm.simulate(n_samples=int(1e4), do={"CVP": "NORMAL", "HR": "HIGH"})
get_distribution(samples, variables=["HISTORY"])  # P(HISTORY)

HISTORY
FALSE    0.9488
TRUE     0.0512
dtype: float64

Simulation with soft/virtual intervention

soft_evidence = [
    TabularCPD(
        "CVP", 3, [[0.2], [0.3], [0.5]], state_names={"CVP": ["NORMAL", "LOW", "HIGH"]}
    )
]
samples = alarm.simulate(n_samples=int(1e4), virtual_intervention=soft_evidence)
get_distribution(samples, variables=["HISTORY"])  # P(HISTORY)

HISTORY
FALSE    0.9508
TRUE     0.0492
dtype: float64