Repository: pgmpy/pgmpy

Issue:
Issue with dirichlet prior for bayesian learning
The issue is with Bayesian Estimator.
http://pgmpy.org/estimators.html#bayesian-estimator
 The prior type can be one of the following: ‘dirichlet’, ‘BDeu’, ‘K2’
My issue is with the dirichlet, and the parameter: "pseudo_counts". It seems like pgmpy only allows for ONE set of pseudo counts-(a0...ak) where k is the range of states for the target variable (the variable in the rows of the cpd). However this is NOT the accepted definition for the dirichlet pseudo counts-there should be a different set of pseudo counts (a0..ak) for every different combination of the targets variable's parents(for every column in the cpd)

For example-in the example you wrote in the documentation:
You only allow one set of pseudo counts(1,2) for each and every one of the columns(combination of C's parents)
╒══════╤══════╤══════╤══════╤════════════════════╕
│ A    │ A(0) │ A(0) │ A(1) │ A(1)               │
├──────┼──────┼──────┼──────┼────────────────────┤
│ B    │ B(0) │ B(1) │ B(0) │ B(1)               │
├──────┼──────┼──────┼──────┼────────────────────┤
│ C(0) │ 0.25 │ 0.25 │ 0.5  │ 0.3333333333333333 │
├──────┼──────┼──────┼──────┼────────────────────┤
│ C(1) │ 0.75 │ 0.75 │ 0.5  │ 0.6666666666666666 │
╘══════╧══════╧══════╧══════╧════════════════════╛
               (1,2)  (1,2)   (1,2)  (1,2)

But this is not the definition of dirichlet. The dirichlet definition allows a DIIFERENT set of hyper parameters(pseudo counts) for every column-like this

╒══════╤══════╤══════╤══════╤════════════════════╕
│ A    │ A(0) │ A(0) │ A(1) │ A(1)               │
├──────┼──────┼──────┼──────┼────────────────────┤
│ B    │ B(0) │ B(1) │ B(0) │ B(1)               │
├──────┼──────┼──────┼──────┼────────────────────┤
│ C(0) │ 0.25 │ 0.25 │ 0.5  │ 0.3333333333333333 │
├──────┼──────┼──────┼──────┼────────────────────┤
│ C(1) │ 0.75 │ 0.75 │ 0.5  │ 0.6666666666666666 │
╘══════╧══════╧══════╧══════╧════════════════════╛
               (1,5)   (1,2)  (7,9)   (8,3)

Of course, having the option to put a different set of hyper parameters for different columns(different combination of the target variable's parents) is a necessity because
 in reality you have different priors for different values of a the target variable's parents. For example, in the project I am working on one of the cpds is for a 
variable called "party_identification"(it has 3 states- republican, democrat, and none), and the parents are variables like age, gender, state etc. Obviously, there are different 
priors for P(party_identification-republican | age<30 and gender=female and state=Northeast) 
and for P(party_identification-republican | 50>age>30 and gender=female and state=MIDWEST), and so for those two columns I need 6 different sets of hyper parameters(pseudo counts):
(a_republican_age<30_gender=female,state=northeast), (a_democarat_age<30_gender=female,state=northeast), (a_none_age<30_gender=female,state=northeast)
(a_republican_50>age<30_gender=female,state=midwest), (a_democrat_50>age<30_gender=female,state=midwest), (a_none_50>age<30_gender=female,state=midwest)

but pgmpy-or at least it seems like-only allows 3 sets of hyper parameters(pseudo counts):
(a_republican), (a_democrat), (a_none)

In general, if I have 3 states for the "age" parent, 2 states for the gender variable and 3 states for the state variable, I will need 
3*2*3(parents states) * 3 (target variable state)=54 pseudo counts, or in different words, 3*2*3=18 triplets(a_rep, a_dem, a_none) of pseudo counts, one for each column.
However pgmpy only all allow 3 pseudo counts, or the same triplet (a_rep_a_dem_a_none) for each and every one of the columns.
In general, no matter how big my cpd is, and how many columns (target variable's parents combination) I have pgmpy will only allow the number of hyperparameters to be the number 
of rows (states of the target variable). But it should enable the number of hyperparameters to be number of rows(number of states of the target variables) * number of columns
(parents combination). Or, in other words, different set of pseudo counts (target_variable_state_1......target_variable_state_k) for EACH AND EVERY ONE of the columns.

As I said both from the documentation and the code it seems like pgmpy only allow one set of pseudo counts-a fixed set of pseudo counts-for all the columns.



Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=h1) Report
> Merging [#973](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=desc) into [dev](https://codecov.io/gh/pgmpy/pgmpy/commit/b85e9ca3b20296925dd4b5686356031dc1f80fec?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `92.85%`.

[![Impacted file tree graph](https://codecov.io/gh/pgmpy/pgmpy/pull/973/graphs/tree.svg?src=pr&token=UaJMCdHaEF&width=650&height=150)](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##              dev     #973   +/-   ##
=======================================
  Coverage   94.73%   94.73%           
=======================================
  Files         114      114           
  Lines       11216    11216           
=======================================
  Hits        10625    10625           
  Misses        591      591
```


| [Impacted Files](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pgmpy/models/DynamicBayesianNetwork.py](https://codecov.io/gh/pgmpy/pgmpy/pull/973/diff?src=pr&el=tree#diff-cGdtcHkvbW9kZWxzL0R5bmFtaWNCYXllc2lhbk5ldHdvcmsucHk=) | `91.8% <ø> (ø)` | :arrow_up: |
| [pgmpy/tests/test\_models/test\_BayesianModel.py](https://codecov.io/gh/pgmpy/pgmpy/pull/973/diff?src=pr&el=tree#diff-cGdtcHkvdGVzdHMvdGVzdF9tb2RlbHMvdGVzdF9CYXllc2lhbk1vZGVsLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [...py/tests/test\_estimators/test\_BayesianEstimator.py](https://codecov.io/gh/pgmpy/pgmpy/pull/973/diff?src=pr&el=tree#diff-cGdtcHkvdGVzdHMvdGVzdF9lc3RpbWF0b3JzL3Rlc3RfQmF5ZXNpYW5Fc3RpbWF0b3IucHk=) | `100% <100%> (ø)` | :arrow_up: |
| [pgmpy/estimators/BayesianEstimator.py](https://codecov.io/gh/pgmpy/pgmpy/pull/973/diff?src=pr&el=tree#diff-cGdtcHkvZXN0aW1hdG9ycy9CYXllc2lhbkVzdGltYXRvci5weQ==) | `91.89% <87.5%> (+2.7%)` | :arrow_up: |
| [pgmpy/models/BayesianModel.py](https://codecov.io/gh/pgmpy/pgmpy/pull/973/diff?src=pr&el=tree#diff-cGdtcHkvbW9kZWxzL0JheWVzaWFuTW9kZWwucHk=) | `95.63% <0%> (-0.4%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=footer). Last update [b85e9ca...c3f11cd](https://codecov.io/gh/pgmpy/pgmpy/pull/973?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

ankurankan: @khalibartan @yashu-seth Reviews please. Urgently if possible.

Patch:
diff --git a/pgmpy/estimators/BayesianEstimator.py b/pgmpy/estimators/BayesianEstimator.py
index 01689a50..4326ef35 100644
--- a/pgmpy/estimators/BayesianEstimator.py
+++ b/pgmpy/estimators/BayesianEstimator.py
@@ -27,8 +27,9 @@ class BayesianEstimator(ParameterEstimator):
         prior_type: 'dirichlet', 'BDeu', or 'K2'
             string indicting which type of prior to use for the model parameters.
             - If 'prior_type' is 'dirichlet', the following must be provided:
-                'pseudo_counts' = dirichlet hyperparameters; a dict containing, for each variable, a list
-                 with a "virtual" count for each variable state, that is added to the state counts.
+                'pseudo_counts' = dirichlet hyperparameters; a dict containing, for each variable, a 2-D
+                 array of the shape (node_card, product of parents_card) with a "virtual" count for each 
+                 variable state in the CPD, that is added to the state counts.
                  (lexicographic ordering of states assumed)
             - If 'prior_type' is 'BDeu', then an 'equivalent_sample_size'
                 must be specified instead of 'pseudo_counts'. This is equivalent to
@@ -65,7 +66,7 @@ class BayesianEstimator(ParameterEstimator):
         for node in self.model.nodes():
             _equivalent_sample_size = equivalent_sample_size[node] if isinstance(equivalent_sample_size, dict) else \
                                       equivalent_sample_size
-            _pseudo_counts = pseudo_counts[node] if isinstance(pseudo_counts, dict) else pseudo_counts
+            _pseudo_counts = pseudo_counts[node] if pseudo_counts else None
 
             cpd = self.estimate_cpd(node,
                                     prior_type=prior_type,
@@ -87,8 +88,9 @@ class BayesianEstimator(ParameterEstimator):
         prior_type: 'dirichlet', 'BDeu', 'K2',
             string indicting which type of prior to use for the model parameters.
             - If 'prior_type' is 'dirichlet', the following must be provided:
-                'pseudo_counts' = dirichlet hyperparameters; a list or dict
-                 with a "virtual" count for each variable state.
+                'pseudo_counts' = dirichlet hyperparameters; 2-D array of shape
+                 (node_card, product of parents_card) with a "virtual" count for 
+                 each variable state in the CPD.
                  The virtual counts are added to the actual state counts found in the data.
                  (if a list is provided, a lexicographic ordering of states is assumed)
             - If 'prior_type' is 'BDeu', then an 'equivalent_sample_size'
@@ -126,22 +128,23 @@ class BayesianEstimator(ParameterEstimator):
         node_cardinality = len(self.state_names[node])
         parents = sorted(self.model.get_parents(node))
         parents_cardinalities = [len(self.state_names[parent]) for parent in parents]
+        cpd_shape = (node_cardinality, np.prod(parents_cardinalities, dtype=int))
 
         if prior_type == 'K2':
-            pseudo_counts = [1] * node_cardinality
+            pseudo_counts = np.ones(cpd_shape, dtype=int)
         elif prior_type == 'BDeu':
             alpha = float(equivalent_sample_size) / (node_cardinality * np.prod(parents_cardinalities))
-            pseudo_counts = [alpha] * node_cardinality
+            pseudo_counts = np.ones(cpd_shape, dtype=float) * alpha
         elif prior_type == 'dirichlet':
-            if not len(pseudo_counts) == node_cardinality:
-                raise ValueError("'pseudo_counts' should have length {0}".format(node_cardinality))
-            if isinstance(pseudo_counts, dict):
-                pseudo_counts = sorted(pseudo_counts.values())
+            pseudo_counts = np.array(pseudo_counts)
+            if pseudo_counts.shape != cpd_shape:
+                raise ValueError("The shape of pseudo_counts must be: {shape}".format(
+                   shape=str(cpd_shape)))
         else:
             raise ValueError("'prior_type' not specified")
 
         state_counts = self.state_counts(node)
-        bayesian_counts = (state_counts.T + pseudo_counts).T
+        bayesian_counts = state_counts + pseudo_counts
 
         cpd = TabularCPD(node, node_cardinality, np.array(bayesian_counts),
                          evidence=parents,
diff --git a/pgmpy/models/DynamicBayesianNetwork.py b/pgmpy/models/DynamicBayesianNetwork.py
index a6ee9513..566dfba6 100644
--- a/pgmpy/models/DynamicBayesianNetwork.py
+++ b/pgmpy/models/DynamicBayesianNetwork.py
@@ -486,7 +486,7 @@ class DynamicBayesianNetwork(DirectedGraph):
         >>> student.add_nodes_from(['D', 'G', 'I', 'S', 'L'])
         >>> student.add_edges_from([(('D', 0),('G', 0)),(('I', 0),('G', 0)),(('D', 0),('D', 1)),(('I', 0),('I', 1))])
         >>> grade_cpd = TabularCPD(('G', 0), 3, [[0.3, 0.05, 0.9, 0.5],
-        ...                                      [0.4, 0.25, 0.8, 0.03],
+        ...                                      [0.4, 0.25, 0.08, 0.3],
         ...                                      [0.3, 0.7, 0.02, 0.2]],
         ...                        evidence=[('I', 0),('D', 0)],
         ...                        evidence_card=[2, 2])


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_BayesianEstimator.py b/pgmpy/tests/test_estimators/test_BayesianEstimator.py
index 86709d5f..fc85a46f 100644
--- a/pgmpy/tests/test_estimators/test_BayesianEstimator.py
+++ b/pgmpy/tests/test_estimators/test_BayesianEstimator.py
@@ -22,19 +22,21 @@ class TestBayesianEstimator(unittest.TestCase):
         self.est3 = BayesianEstimator(self.m1, self.d2)
 
     def test_estimate_cpd_dirichlet(self):
-        cpd_A = self.est1.estimate_cpd('A',  prior_type="dirichlet", pseudo_counts=[0, 1])
+        cpd_A = self.est1.estimate_cpd('A',  prior_type="dirichlet", pseudo_counts=[[0], [1]])
         self.assertEqual(cpd_A, TabularCPD('A', 2, [[0.5], [0.5]]))
 
-        cpd_B = self.est1.estimate_cpd('B',  prior_type="dirichlet", pseudo_counts=[9, 3])
+        cpd_B = self.est1.estimate_cpd('B',  prior_type="dirichlet", pseudo_counts=[[9], [3]])
         self.assertEqual(cpd_B, TabularCPD('B', 2, [[11.0/15], [4.0/15]]))
 
-        cpd_C = self.est1.estimate_cpd('C',  prior_type="dirichlet", pseudo_counts=[0.4, 0.6])
+        cpd_C = self.est1.estimate_cpd('C',  prior_type="dirichlet", pseudo_counts=[[0.4, 0.4, 0.4, 0.4],
+                                                                                    [0.6, 0.6, 0.6, 0.6]])
         self.assertEqual(cpd_C, TabularCPD('C', 2, [[0.2, 0.2, 0.7, 0.4],
                                                     [0.8, 0.8, 0.3, 0.6]],
                                            evidence=['A', 'B'], evidence_card=[2, 2]))
 
     def test_estimate_cpd_improper_prior(self):
-        cpd_C = self.est1.estimate_cpd('C',  prior_type="dirichlet", pseudo_counts=[0, 0])
+        cpd_C = self.est1.estimate_cpd('C',  prior_type="dirichlet", pseudo_counts=[[0, 0, 0, 0],
+                                                                                    [0, 0, 0, 0]])
         cpd_C_correct = (TabularCPD('C', 2, [[0.0, 0.0, 1.0, np.NaN],
                                              [1.0, 1.0, 0.0, np.NaN]],
                                     evidence=['A', 'B'], evidence_card=[2, 2],
@@ -64,7 +66,8 @@ class TestBayesianEstimator(unittest.TestCase):
         self.assertSetEqual(set(self.est3.get_parameters()), cpds)
 
     def test_get_parameters2(self):
-        pseudo_counts = {'A': [1, 2, 3], 'B': [4, 5], 'C': [6, 7]}
+        pseudo_counts = {'A': [[1], [2], [3]], 'B': [[4], [5]], 'C': [[6, 6, 6, 6, 6, 6],
+                                                                      [7, 7, 7, 7, 7, 7]]}
         cpds = set([self.est3.estimate_cpd('A', prior_type="dirichlet", pseudo_counts=pseudo_counts['A']),
                     self.est3.estimate_cpd('B', prior_type="dirichlet", pseudo_counts=pseudo_counts['B']),
                     self.est3.estimate_cpd('C', prior_type="dirichlet", pseudo_counts=pseudo_counts['C'])])
diff --git a/pgmpy/tests/test_models/test_BayesianModel.py b/pgmpy/tests/test_models/test_BayesianModel.py
index c2daf5db..dbc4d480 100644
--- a/pgmpy/tests/test_models/test_BayesianModel.py
+++ b/pgmpy/tests/test_models/test_BayesianModel.py
@@ -431,7 +431,11 @@ class TestBayesianModelFitPredict(unittest.TestCase):
     def test_bayesian_fit(self):
         print(isinstance(BayesianEstimator, BaseEstimator))
         print(isinstance(MaximumLikelihoodEstimator, BaseEstimator))
-        self.model2.fit(self.data1, estimator=BayesianEstimator, prior_type="dirichlet", pseudo_counts=[9, 3])
+        self.model2.fit(self.data1, estimator=BayesianEstimator, prior_type="dirichlet",
+                        pseudo_counts={'A': [[9], [3]],
+                                       'B': [[9], [3]],
+                                       'C': [[9, 9, 9, 9],
+                                             [3, 3, 3, 3]]})
         self.assertEqual(self.model2.get_cpds('B'), TabularCPD('B', 2, [[11.0 / 15], [4.0 / 15]]))
 
     def test_fit_missing_data(self):