This PR deals with the basic continuous node representation feature. It will comprise a base class for Continuous node representation along with various methods to discretize the continuous variables into discrete factors. This involves the first three weeks of my GSoC project.

This PR deals with implementing HMC with dual averaging. The implementation is still open for discussion. If you find anything ambiguous please comment on the line. This PR is open for discussion.

This PR deals with

• the creation of a base class ContinuousFactor for multivariate representations.
• the creation of the class JointGaussainDistribution - a model to represent the gaussain random variables.

All functions are not implemented. only get_variable,get_states and get_property are implemented.

Creating this pull request for easy review. I have tested these methods on munin2.bif and dog-problem.bif and they are working fine. Implemented these functions in accordance with BIF v0.15 as given here . Time taken to run if not importing pgmpy and numpy was 0.06s on average with printing the complete variable_states for munin2.bif (PS: 1003 nodes) Check issue #506 .

@ankurankan I have removed cardinalities from the model attributes but it seems it is being used at other places as well. Should cardinalities be computed every single time it is required? Is there a particular problem in having it as an attribute?

Here, we finish the implementation of VE with few missing steps, namely computing good elimination orderings (with 4 heuristics - min neighbors, min fill, min weight, and weighted min fill) and safely removing irrelevant variables from the model (barren and independent by evidence nodes). In order to make these improvements, few new methods were necessary and modifying other members too. In our preliminary experimental results, queries that were taking up to 30 minutes using VE now take less than 2 minutes. Please, help us test this new code for a robust implementation and also looking for better ways of coding the algorithms. Thanks.

• Replaces the recursive version of fun with an iterative one, much easier to read in my opinion. This should be slightly fast as well because recursive calls are expensive in Python(results in new stack frame each time), plus Python has a default limit of 1000 for recursive calls(though this is highly unlikely to occur in our case).
• I am not sure why the library is using Python 2 based super() calls considering the fact that dependency includes Python 3.3. In Python 3 thanks to the cell variable __class__ we can simply use super().method_name(...)(PEP 3135 -- New Super).

Added a new method that gives probabilities of missing variables given a predict data #794

        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`

Also a new error test for predict that increases the coverage.

Subject of the issue

I want to reproduce the example here

Your environment

• pgmpy version: 0.1.12
• Python version: 3.6.9
• Operating System: Ubuntu 18.04.5 LTS

Steps to reproduce

import pandas as pd
import numpy as np
from pgmpy.estimators import HillClimbSearch, BicScore
data = pd.DataFrame(np.random.randint(0, 5, size=(5000, 9)), columns=list('ABCDEFGHI'))
data['J'] = data['A'] * data['B']
est = HillClimbSearch(data, scoring_method=BicScore(data))
best_model = est.estimate()
best_model.edges()

Expected behaviour

[('B', 'J'), ('A', 'J')]

Actual behaviour

[('A', 'B'), ('J', 'A'), ('J', 'B')]

@ankurankan I have send this PR to aid us in discussion. I was experimenting things with how to handle gradients (removing the gradient argument). I tested with two ways:

• First I tried handled grad_log_pdf argument on places itself depending upon how user passed the argument, if None was passed then I created a lambda function to call model.get_gradient_log_pdf otherwise I created a lambda function to use the custom class. But with this things were messy as I have to handle this parameter at two places, one in sampling class and other in BaseSimulateHamiltonianDynamics class.
• Second ( this PR implements it). Handle everything in model.get_gradient_log_pdf. This code is less messy, because every call is made to model.get_gradient_log_pdf and the method internally handles the rest so need of making suitable changes at different places.

How do you suggest I should handle the gradients ? You can look at the last commit to specifically see the changes I made https://github.com/pgmpy/pgmpy/pull/702/commits/748eb1fe13488bb8f0cf27a7064a67384ec3315e

After the discussion I'll close one of the PR.

A factor product following "Probabilistic Graphical Models" (Koller 09) on page 359, Algorithm 10.A.1 - Efficient implementation of a factor product operation. Koller's algorithm was modified to fit the configuration used in pgmpy. For example, in pgmpy the configurations of Factor are supposed to be like (0,0,0) (0,0,1) (0,1,0) (1,0,0) and so on, instead of (0,0,0) (1,0,0) (0,1,0) (0,0,1) as expected for Koller's algorithm.

Koller's implementation is around 98% faster than the current one in pgmpy. This benchmark was done by using a simple python script as follows:

from pgmpy.factors import Factor
from pgmpy.factors import factor_product
from time import time

phi = Factor(['x1', 'x2'], [2, 2], range(4))
phi1 = Factor(['x3', 'x4'], [2, 2], range(4))
t0 = time()
prod = factor_product(phi, phi1)
t1 = time()
print(t1-t0)

After running 6 time each implementation, here is the results:

Unfortunately, we don't know how to use JobLib. But we leave this TODO with the hope that using parallel computation can improve this implementation even further.

This is a question rather than an issue. I was wondering if it is possible to update the edges and dag in a training procedure! I could not find any related method if there is one.

Would you please help me with this question?

I appreciate it, Thank you so much in advance.

Subject of the issue

A Dynamic Bayesian network is built, but the network is always displayed: CPD associated with (B, 1) doesn't have proper parents associated with it.

Your environment

• pgmpy version

Steps to reproduce

from pgmpy.factors.discrete import TabularCPD from pgmpy.models import DynamicBayesianNetwork as DBN from pgmpy.inference import DBNInference

dbnet = DBN() dbnet.add_edges_from( [(('A', 0), ('B', 0)), (('A', 0), ('C', 0)), (('C', 0), ('D', 0)), (('B', 0), ('B', 1)), (('B', 0), ('D', 1 )) ] ) a_cpds = TabularCPD(('A', 0), 2, [[0.7], [0.3]]) b_start_cpds = TabularCPD( ('B', 0), 2, [[0.3, 0.6], [0.7, 0.4]], evidence=[('A', 0)], evidence_card=[2] ) b_trans_cpds = TabularCPD( ('B',1), 2, [[0.1, 0.3, 0.8, 0.6], [0.9, 0.7, 0.2, 0.4]], evidence=[('A', 0), ('B', 0)], evidence_card=[2, 2] ) c_cpds = TabularCPD( ('C', 0), 2, [[0.3, 0.1], [0.7,0.9]], evidence=[('A', 0)], evidence_card=[2] ) d_start_cpds = TabularCPD( ('D',0), 2, [[0.4, 0.2], [0.6, 0.8]], evidence=[('C', 0)], evidence_card=[2] ) d_trans_cpds = TabularCPD( ('D',1), 2, [[0.3, 0.4, 0.8, 0.9], [0.7, 0.6, 0.2, 0.1]], evidence=[('B', 0),('C',0)], evidence_card=[2, 2] ) dbnet.add_cpds(a_cpds, b_start_cpds, b_trans_cpds, c_cpds, d_start_cpds, d_trans_cpds) dbnet.initialize_initial_state() dbn_inf = DBNInference(dbnet) temp = dbn_inf.query([('D', 1)], {('D', 0):0, ('D', 2):1})['D', 1].values print(temp)

Expected behaviour

When the evidence under the first time slice of D is true and the evidence under the third time slice is false, the probability of the second time slice is calculated by reasoning. The final calculation result is D2, True: 0.7365902, False:0.2634098

Actual behaviour

CPD associated with (B, 1) doesn't have proper parents associated with it. Is this a problem of network model building? Specific network diagram, see [https://www.bilibili.com/video/BV1W3411M7cp/?spm_id_from=333.337.search-card.all.click&vd_source=1a42ddbfb403e da9de41b20ccdca8523]

Subject of the issue

If I try to run the GibbsSampling for a dataset, it sometimes succeeds, and sometimes fails. This seems to be because DiscreteFactor uses set() in the product() which was introduced in #1548 in July by @ankurankan and the order of set() changes for every execution, which causes phi.variable to not equal phi.variables[0] sometimes.

I think this can be solved by replacing

    new_variables = list(set(phi.variables).union(phi1.variables))

by something like

    new_variables = phi.variables + [var for var in phi1.variables if var not in phi.variables]

Your environment

• pgmpy version : 0.1.20 and also dev branch (as of november 2022)
• Python version : 3.10.8
• Operating System : MacOS

Steps to reproduce

The following unit test sometimes succeeds completely, and sometimes fails completely.

import unittest

from pgmpy.factors import factor_product
from pgmpy.factors.discrete import TabularCPD
from pgmpy.models import BayesianNetwork
from pgmpy.sampling import GibbsSampling


class TestGibbsSamplingIssue(unittest.TestCase):
    def setUp(self):
        self.cpt_cloudy = TabularCPD(variable='Cloudy', variable_card=2, values=[[0.5], [0.5]])
        self.cpt_sprinkler = TabularCPD(variable='Sprinkler', variable_card=2,
                                        values=[[0.5, 0.9], [0.5, 0.1]],
                                        evidence=['Cloudy'], evidence_card=[2])
        self.cpt_rain = TabularCPD(variable='Rain', variable_card=2,
                                   values=[[0.8, 0.2], [0.2, 0.8]],
                                   evidence=['Cloudy'], evidence_card=[2])
        self.cpt_wet_grass = TabularCPD(variable='Wet_Grass', variable_card=2,
                                        values=[[1, 0.1, 0.1, 0.01],
                                                [0, 0.9, 0.9, 0.99]],
                                        evidence=['Sprinkler', 'Rain'],
                                        evidence_card=[2, 2])

        self.dag = BayesianNetwork(
            [('Cloudy', 'Sprinkler'), ('Cloudy', 'Rain'), ('Sprinkler', 'Wet_Grass'), ('Rain', 'Wet_Grass')]
        )
        self.dag.add_cpds(self.cpt_cloudy, self.cpt_sprinkler, self.cpt_rain, self.cpt_wet_grass)

    def test_sampling_gibbs(self):
        # Initial issue
        GibbsSampling(self.dag).sample()

    def test_factor_product(self):
        # Trying to get to the root cause
        factor_product(self.cpt_cloudy, self.cpt_sprinkler, self.cpt_rain)

    def test_product(self):
        self.cpt_cloudy * self.cpt_sprinkler * self.cpt_rain

    def test_mul(self):
        prod = self.cpt_cloudy * self.cpt_sprinkler
        # the order or variables get mixed up, which I believe should not happen
        assert prod.variable == prod.variables[0]

Expected behaviour

It should not throw an error sometimes

Actual behaviour

Traceback (most recent call last):

  File "./pgmpy/pgmpy/factors/discrete/CPD.py", line 299, in copy
    return TabularCPD(
  File "./pgmpy/pgmpy/factors/discrete/CPD.py", line 142, in __init__
    super(TabularCPD, self).__init__(
  File "./pgmpy/pgmpy/factors/discrete/DiscreteFactor.py", line 99, in __init__
    raise ValueError("Variable names cannot be same")
ValueError: Variable names cannot be same

List of upcoming enhancements to the Causal Inference class:

• [ ] estimate_ate and query methods should work with frontdoor adjustment sets.
• [ ] Remove averaging over different adjustment sets for computing ATE.
• [ ] Extend instrumental variables (IVs), conditional IVs, instrumental sets, and conditional instrumental sets from SEMs to work with DAGs as well.
• [ ] Extend estimation methods to work with IVs and Conditional IVs.

Subject of the issue (with proposed fix)

When creating multiple MaximumLikelihoodEstimator objects and calling estimate_cpd for all the model nodes in each, we seem to get a memory leak. This is what happens with repeated calls to BayesianNetwork.fit.

The state_counts function of the BaseEstimator class has a decorator ([https://github.com/pgmpy/pgmpy/blob/dev/pgmpy/estimators/base.py#L66])

@lru_cache(maxsize=2048)

which seems to be the culprit. When I remove it, the memory leak goes away. I haven't bothered to dig deeper on the issue, but it's something that may be worth considering if repeatedly training models.

Your environment

• Google Colab Notebook
• PGMPY version: pgmpy-0.1.20
• Python version: 3.7.15
• Operating System: Ubuntu 18.04.6 LTS (Bionic Beaver)
• uname -a output: Linux 5.10.133+

Steps to reproduce

Try the following code out. You should see the memory slowly get eaten up.

import pandas as pd
from pgmpy.models import BayesianNetwork
from pgmpy.estimators import MaximumLikelihoodEstimator

data = pd.read_csv('https://raw.githubusercontent.com/mwaskom/seaborn-data/master/iris.csv')

cols = list(data.columns)
edgelist = [(a,b) for a in cols[0:3] for b in cols[3:] if a != b]

model = BayesianNetwork()
model.add_nodes_from(cols)
model.add_edges_from(edgelist)


for i in range(1000):
  model.fit(data, MaximumLikelihoodEstimator)

Expected behaviour

Repeatedly calling the model.fit function with the MLE should not consume memory ad infinitum, the model gets re-trained and so old data should not persist.

Actual behaviour

During the 1000 iterations additional RAM is consumed and not freed EVER.

Subject of the issue

Belief Propagation tries to allocate exponential memory on factor trees. This shouldnt be happening since trees are still without cycles and therefore efficient for exact marginalization.

Your environment

• pgmpy version 0.1.19
• Python version 3.8.10
• Operating System Ubuntu 20.04

Steps to reproduce

import pgmpy.models
import numpy as np
import pgmpy.factors.discrete
import networkx
import matplotlib.pyplot as plt
from networkx.drawing.nx_agraph import graphviz_layout
import pgmpy.inference


# parameters for the example. Adjust this to see when it breaks
number_of_timesteps = 20
cardinality_of_timesteps = 5

# initialize random transition model
transition_model = np.random.uniform(low=0., high=1., size=(pow(cardinality_of_timesteps,2),))
transition_model /= sum(transition_model)


# create factorgraph
factor_graph = pgmpy.models.FactorGraph()

# add variable nodes for timesteps
timesteps = ["t%s" % t for t in range(number_of_timesteps)]
factor_graph.add_nodes_from(timesteps)

# create transition factors
factors = []

# for each transition
for idx in range(len(timesteps)-1):

    # get the variable names
    state_names = {"t%s" % idx: list(range(cardinality_of_timesteps)),
                    "t%s" % (idx+1): list(range(cardinality_of_timesteps))}

    # create factor with values from transition model
    factor = pgmpy.factors.discrete.DiscreteFactor(list(state_names.keys()),
                                                    [cardinality_of_timesteps,
                                                    cardinality_of_timesteps], 
                                                    transition_model, 
                                                    state_names)
    factors.append(factor)

# add factors
factor_graph.add_factors(*factors)

# add edges for state variables and transition variables
for idx, factor in enumerate(factors):
    factor_graph.add_edges_from([("t%s" % idx, factor),
                                    (factor, "t%s" % (idx+1))])

# create prior factors
for timestep in timesteps:

    # create values of current variable
    state_names = {timestep: list(range(cardinality_of_timesteps))}

    # create random evidence
    evidence = np.random.uniform(low=0.1, high=1., size=(cardinality_of_timesteps,))
    evidence /= sum(evidence)

    # create a factor from it
    factor = pgmpy.factors.discrete.DiscreteFactor([timestep], [cardinality_of_timesteps],
                                                    evidence, state_names)

    # add factor and edge from variable to prior
    factor_graph.add_factors(factor)
    factor_graph.add_edge(timestep, factor)


# plot the structure
g = networkx.Graph(factor_graph.edges())
pos = graphviz_layout(g, prog='dot')
networkx.draw(g, pos, with_labels=False, arrows=False)
# plt.show()

bp = pgmpy.inference.BeliefPropagation(factor_graph)
bp.calibrate()
independent_marginals = bp.query(factor_graph.get_variable_nodes(), joint=False)
print(independent_marginals)

Expected behaviour

Calculation of the independent marginals in polynomial time and memory consumption.

Actual behaviour

Trying to allocate exponential big memory numpy.core._exceptions.MemoryError: Unable to allocate 45.5 GiB for an array with shape (5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5) and data type float64