<jupyter_start><jupyter_text>1. Learning Bayesian NetworksPrevious notebooks showed how Bayesian networks economically encode a probability distribution over a set of variables, and how they can be used to predict variable states, or to generate new samples from the joint distribution. This section will be about obtaining a Bayesian network, given a set of sample data. Learning a Bayesian network can be split into two problems:**Parameter learning**: Given a set of data samples and a DAG that captures the dependencies between the variables, estimate the (conditional) probability distributions of the individual variables.**Structure learning**: Given a set of data samples, estimate a DAG that captures the dependencies between the variables.This notebook aims to illustrate how parameter learning and structure learning can be done with pgmpy. Currently, the library supports:* Parameter learning for discrete nodes:    * Maximum Likelihood Estimation    * Bayesian Estimation* Structure learning for discrete, fully observed networks:    * Score-based structure estimation (BIC/BDeu/K2 score; exhaustive search, hill climb/tabu search)    * Constraint-based structure estimation (PC)    * Hybrid structure estimation (MMHC)     1.1 Parameter Learning Suppose we have the following data:<jupyter_code>import pandas as pd
data = pd.DataFrame(data={'fruit': ["banana", "apple", "banana", "apple", "banana","apple", "banana", 
                                    "apple", "apple", "apple", "banana", "banana", "apple", "banana",], 
                          'tasty': ["yes", "no", "yes", "yes", "yes", "yes", "yes", 
                                    "yes", "yes", "yes", "yes", "no", "no", "no"], 
                          'size': ["large", "large", "large", "small", "large", "large", "large",
                                    "small", "large", "large", "large", "large", "small", "small"]})
print(data)<jupyter_output>fruit   size tasty
0   banana  large   yes
1    apple  large    no
2   banana  large   yes
3    apple  small   yes
4   banana  large   yes
5    apple  large   yes
6   banana  large   yes
7    apple  small   yes
8    apple  large   yes
9    apple  large   yes
10  banana  large   yes
11  banana  large    no
12   apple  small    no
13  banana  small    no<jupyter_text>We know that the variable relate as follows:<jupyter_code>from pgmpy.models import BayesianModel

model = BayesianModel([('fruit', 'tasty'), ('size', 'tasty')])  # fruit -> tasty <- size<jupyter_output><empty_output><jupyter_text>Parameter learning is the task to estimate the values of the conditional probability distributions (CPDs), for the variables fruit, size, and tasty.**State counts**To make sense of the given data, we can start by counting how often each state of the variable occurs. If the variable is dependent on parents, the counts are done conditionally on the parents states, i.e. seperately for each parent configuration:<jupyter_code>from pgmpy.estimators import ParameterEstimator
pe = ParameterEstimator(model, data)
print("\n", pe.state_counts('fruit'))  # unconditional
print("\n", pe.state_counts('tasty'))  # conditional on fruit and size<jupyter_output>fruit
apple       7
banana      7

 fruit apple       banana      
size  large small  large small
tasty                         
no      1.0   1.0    1.0   1.0
yes     3.0   2.0    5.0   0.0<jupyter_text>We can see, for example, that as many apples as bananas were observed and that 5 large bananas were tasty, while only 1 was not. **Maximum Likelihood Estimation**A natural estimate for the CPDs is to simply use the relative frequencies, with which the variable states have occured. We observed 7 `apples` among a total of 14 `fruits`, so we might guess that about 50% of `fruits` are `apples`.This approach is *Maximum Likelihood Estimation (MLE)*. According to MLE, we should fill the CPDs in such a way, that $P(\text{data}|\text{model})$ is maximal. This is achieved when using the relative frequencies. pgmpy supports MLE as follows:<jupyter_code>from pgmpy.estimators import MaximumLikelihoodEstimator
mle = MaximumLikelihoodEstimator(model, data)
print(mle.estimate_cpd('fruit'))  # unconditional
print(mle.estimate_cpd('tasty'))  # conditional<jupyter_output>╒═══════════════╤═════╕
│ fruit(apple)  │ 0.5 │
├───────────────┼─────┤
│ fruit(banana) │ 0.5 │
╘═══════════════╧═════╛
╒════════════╤══════════════╤════════════════════╤═════════════════════╤═══════════════╕
│ fruit      │ fruit(apple) │ fruit(apple)       │ fruit(banana)       │ fruit(banana) │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ size       │ size(large)  │ size(small)        │ size(large)         │ size(small)   │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ tasty(no)  │ 0.25         │ 0.3333333333333333 │ 0.16666666666666666 │ 1.0           │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ tasty(yes) │ 0.75         │ 0.6666666666666666 │ 0.8333333333333334  │ 0.0           │
╘════════════╧══════════════╧════════════════════╧═════════════════════╧═══════════════╛<jupyter_text>`mle.estimate_cpd(variable)` computes the state counts and divides each cell by the (conditional) sample size. The `mle.get_parameters()`-method returns a list of CPDs for all variables of the model.The built-in `fit()`-method of `BayesianModel` provides more convenient access to parameter estimators:<jupyter_code># Calibrate all CPDs of `model` using MLE:
model.fit(data, estimator=MaximumLikelihoodEstimator)<jupyter_output><empty_output><jupyter_text>While very straightforward, the ML estimator has the problem of overfitting to the data. In above CPD, the probability of a large banana being tasty is estimated at 0.833, because 5 out of 6 observed large bananas were tasty. Fine. But note that the probability of a small banana being tasty is estimated at 0.0, because we observed only one small banana and it happened to be not tasty. But that should hardly make us certain that small bananas aren't tasty! We simply do not have enough observations to rely on the observed frequencies. If the observed data is not representative for the underlying distribution, ML estimations will be extremly far off.When estimating parameters for Bayesian networks, lack of data is a frequent problem. Even if the total sample size is very large, the fact that state counts are done conditionally for each parents configuration causes immense fragmentation. If a variable has 3 parents that can each take 10 states, then state counts will be done seperately for 10^3 = 1000 parents configurations. This makes MLE very fragile and unstable for learning Bayesian Network parameters. A way to mitigate MLE's overfitting is *Bayesian Parameter Estimation*. **Bayesian Parameter Estimation**:The Bayesian Parameter Estimator starts with already existing prior CPDs, that express our beliefs about the variables before the data was observed. Those "priors" are then updated, using the state counts from the observed data.One can think of the priors as consisting in pseudo state counts, that are added to the actual counts before normalization. Unless one wants to encode specific beliefs about the distributions of the variables, one commonly chooses uniform priors, i.e. ones that deem all states equiprobable.A very simple prior is the so-called K2 prior, which simply adds 1 to the count of every single state. A somewhat more sensible choice of prior is BDeu (Bayesian Dirichlet equivalent uniform prior). For BDeu we need to specify an equivalent sample size N and then the pseudo-counts are the equivalent of having observed N uniform samples of each variable (and each parent configuration). In pgmpy:<jupyter_code>from  pgmpy.estimators  import BayesianEstimator
est = BayesianEstimator(model, data)

print(est.estimate_cpd('tasty', prior_type='BDeu', equivalent_sample_size=10))<jupyter_output>╒════════════╤═════════════════════╤════════════════════╤════════════════════╤═════════════════════╕
│ fruit      │ fruit(apple)        │ fruit(apple)       │ fruit(banana)      │ fruit(banana)       │
├────────────┼─────────────────────┼────────────────────┼────────────────────┼─────────────────────┤
│ size       │ size(large)         │ size(small)        │ size(large)        │ size(small)         │
├────────────┼─────────────────────┼────────────────────┼────────────────────┼─────────────────────┤
│ tasty(no)  │ 0.34615384615384615 │ 0.4090909090909091 │ 0.2647058823529412 │ 0.6428571428571429  │
├────────────┼─────────────────────┼────────────────────┼────────────────────┼─────────────────────┤
│ tasty(yes) │ 0.6538461538461539  │ 0.5909090909090909 │ 0.7352941176470589 │ 0.35714285714285715 │
╘════════════╧═════════════════════╧════════════════════╧════════════════════╧═════════════════════╛<jupyter_text>The estimated values in the CPDs are now more conservative. In particular, the estimate for a small banana being not tasty is now around 0.64 rather than 1.0. Setting equivalent_sample_size to 10 means that for each parent configuration, we add the equivalent of 10 uniform samples (here: +5 small bananas that are tasty and +5 that aren't).For a quick comparison, recall that the MLE parameters were as follows:<jupyter_code>print(mle.estimate_cpd('tasty'))  # conditional<jupyter_output>╒════════════╤══════════════╤════════════════════╤═════════════════════╤═══════════════╕
│ fruit      │ fruit(apple) │ fruit(apple)       │ fruit(banana)       │ fruit(banana) │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ size       │ size(large)  │ size(small)        │ size(large)         │ size(small)   │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ tasty(no)  │ 0.25         │ 0.3333333333333333 │ 0.16666666666666666 │ 1.0           │
├────────────┼──────────────┼────────────────────┼─────────────────────┼───────────────┤
│ tasty(yes) │ 0.75         │ 0.6666666666666666 │ 0.8333333333333334  │ 0.0           │
╘════════════╧══════════════╧════════════════════╧═════════════════════╧═══════════════╛<jupyter_text>And, if we were to keep increasing the `equivalent_sample_size`, we can see that our probabilities begin to approach 50/50.<jupyter_code>print(est.estimate_cpd('tasty', prior_type='BDeu', equivalent_sample_size=100))
print(est.estimate_cpd('tasty', prior_type='BDeu', equivalent_sample_size=1000))<jupyter_output>╒════════════╤═════════════════════╤═════════════════════╤═════════════════════╤════════════════════╕
│ fruit      │ fruit(apple)        │ fruit(apple)        │ fruit(banana)       │ fruit(banana)      │
├────────────┼─────────────────────┼─────────────────────┼─────────────────────┼────────────────────┤
│ size       │ size(large)         │ size(small)         │ size(large)         │ size(small)        │
├────────────┼─────────────────────┼─────────────────────┼─────────────────────┼────────────────────┤
│ tasty(no)  │ 0.46551724137931033 │ 0.48214285714285715 │ 0.43548387096774194 │ 0.5192307692307693 │
├────────────┼─────────────────────┼─────────────────────┼─────────────────────┼────────────────────┤
│ tasty(yes) │ 0.5344827586206896  │ 0.5178571428571429  │ 0.5645161290322581  │ 0.4807692307692308 │
╘════════════╧═════════════════════╧═════════════════════╧═════════════════════╧════════════════════╛
╒════════════╤═════════════════════╤════════════════════╤═══════════════╤═════════[...]<jupyter_text>`BayesianEstimator`, too, can be used via the `fit()`-method. For example:<jupyter_code>import numpy as np
import pandas as pd
from pgmpy.models import BayesianModel
from pgmpy.estimators import BayesianEstimator

# generate data
data = pd.DataFrame(np.random.randint(low=0, high=2, size=(5000, 4)), columns=['A', 'B', 'C', 'D'])
model = BayesianModel([('A', 'B'), ('A', 'C'), ('D', 'C'), ('B', 'D')])

model.fit(data, estimator=BayesianEstimator, prior_type="BDeu") # default equivalent_sample_size=5
for cpd in model.get_cpds():

  print(cpd)<jupyter_output>╒══════╤══════════╕
│ A(0) │ 0.499401 │
├──────┼──────────┤
│ A(1) │ 0.500599 │
╘══════╧══════════╛
╒══════╤══════╤══════╕
│ A    │ A(0) │ A(1) │
├──────┼──────┼──────┤
│ B(0) │ 0.5  │ 0.5  │
├──────┼──────┼──────┤
│ B(1) │ 0.5  │ 0.5  │
╘══════╧══════╧══════╛
╒══════╤════════════════════╤════════════════════╤════════════════════╤════════════════════╕
│ A    │ A(0)               │ A(0)               │ A(1)               │ A(1)               │
├──────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┤
│ D    │ D(0)               │ D(1)               │ D(0)               │ D(1)               │
├──────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┤
│ C(0) │ 0.5170641846034516 │ 0.5086758934104524 │ 0.4921951219512195 │ 0.5118439861139473 │
├──────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┤
│ C(1) │ 0.4829358153965484 │ 0.4913241065895476 │ 0.5078048780487805 │ 0.48815601388605[...]
{'fraction_non_alphanumeric': 0.269618320610687, 'fraction_numerical': 0.05419847328244275, 'mean_word_length': 4.725961538461538, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 8, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '634863', 'license': 'MIT', 'n_tokens_mistral': 5401, 'n_tokens_neox': 4082, 'n_words': 1339, 'path': 'notebooks/Machine-Learning-Perspective/Probabilistic Graphical Models/Python/Section 3 - Learning Bayesian Networks.ipynb', 'repo_name': 'NathanielDake/NathanielDake.github.io', 'size': '13100'}