Skip to content

Tutorial pyAgrum

Creative Commons LicenseaGrUMinteractive online version

Creating your first Bayesian network with pyAgrum

Section titled “Creating your first Bayesian network with pyAgrum”

(This example is based on an OpenBayes [closed] website tutorial)

A Bayesian network (BN) is composed of random variables (nodes) and their conditional dependencies (arcs) which, together, form a directed acyclic graph (DAG). A conditional probability table (CPT) is associated with each node. It contains the conditional probability distribution of the node given its parents in the DAG:

PyAgrum inline image

Such a BN allows to manipulate the joint probability P(C,S,R,W)P(C,S,R,W)   using this decomposition :

P(C,S,R,W)=XP(XParentsX)=P(C)P(SC)P(RC)P(WS,R)P(C,S,R,W)=\prod_X P(X | Parents_X) = P(C) \cdot P(S | C) \cdot P(R | C) \cdot P(W | S,R)

Imagine you want to create your first Bayesian network, say for example the ‘Water Sprinkler’ network. This is an easy example. All the nodes are Boolean (only 2 possible values). You can proceed as follows.

import pyagrum as gum

The next line creates an empty BN network with a ‘name’ property.

bn = gum.BayesNet("WaterSprinkler")
print(bn)
BN{nodes: 0, arcs: 0, domainSize: 1, dim: 0, mem: 0o}

pyAgrum(aGrUM) provides 5 types of variables :

  • LabelizedVariable {low|medium|high} : a discrete variable that can take 3 values.
  • RangeVariable [3,10] : a discrete variable that can take 8 values: from 3 to 10.
  • IntegerVariable {1|3|4|7} : a discrete variable that can take 4 values
  • NumericalDiscreteVariable {1.0|3.14|4.6|7} : another discrete
  • variable that can take 4 values
  • DiscretizedVariable [1.0,3.0,6.0] :a continuous variable discretized in 2 values 1-3 and 3-6.

In this tutorial, we will use LabelizedVariable, which is a variable whose domain is a finite set of labels. The next line will create a variable named ‘c’, with 2 values and described as ‘cloudy?’, and it will add it to the BN. The value returned is the id of the node in the graphical structure (the DAG). pyAgrum actually distinguishes the random variable (here the labelizedVariable) from its node in the DAG: the latter is identified through a numeric id. Of course, pyAg rum provides functions to get the id of a node given the corresponding variable and conversely.

id_c = bn.add(gum.LabelizedVariable("c", "cloudy ?", 2))
print(id_c)
print(bn.variable(id_c))
0
c:Labelized({0|1})

You can go on adding nodes in the network this way. Let us use python to compact a little bit the code:

id_s, id_r, id_w = [
bn.add(name, 2) for name in "srw"
] # bn.add(name, 2) === bn.add(gum.LabelizedVariable(name, name, 2))
print(id_s, id_r, id_w)
print(bn)
1 2 3
BN{nodes: 4, arcs: 0, domainSize: 16, dim: 4, mem: 64o}

When you add a variable to a graphical model, this variable is “owned” by the model but you still have (read*only) accessors between the 3 related objects : the variable itself, its name and its id in the BN :

print(f"{bn.variable(id_s)=}")
print(f"{bn['s']=}")
print(f"{id_s=}")
print(f"{bn.idFromName('s')=}")
print(f"{bn.variable(id_s).name()=}")
bn.variable(id_s)=(pyagrum.DiscreteVariable@0xb70759020) s:Range([0,1])
bn['s']=(pyagrum.DiscreteVariable@0xb70759020) s:Range([0,1])
id_s=1
bn.idFromName('s')=1
bn.variable(id_s).name()='s'

Now we have to connect nodes, i.e., to add arcs linking the nodes. Remember that id_c and id_s are ids for nodes. But we use directly access by name "c" or "s"

bn.addArc("c", "s")

For adding multiple arcs, o nce again, python can help us :

for link in [(id_c, id_r), ("s", "w"), ("r", "w")]:
bn.addArc(*link)
print(bn)
BN{nodes: 4, arcs: 4, domainSize: 16, dim: 9, mem: 144o}

pyAgrum provides tools to display bn in more user-fri endly fashions.
Notably, pyagrum.lib is a set of tools written in pyAgrum to help using aGrUM in python. pyagrum.lib.notebook adds dedicated functions for iPython notebook.

import pyagrum.lib.notebook as gnb
bn
G c c r r c->r s s c->s w w r->w s->w

The functions fast[model] encode the structure of the graphical model and the type of the variables in a concise language somehow derived from the dot language for graphs (see the doc for the underlying method : fastPrototype).

bn2 = gum.fastBN("c->r->w<-s<-c", "{No|Yes}") # by default, variable is labelized {No|Yes}
bn2
G c c r r c->r s s c->s w w r->w s->w
print(bn2.variable("c"))
print(bn.variable("c"))
bn = bn2
c:Labelized({No|Yes})
c:Labelized({0|1})

Once the network topology is constructed, we must initialize the conditional probability tables (CPT) distributions. Each CPT is considered as a Tensor object in pyagrum. There are several ways to fill such an object.

To get the CPT of a variable, use the cpt method of your BayesNet instance with the variable’s id as parameter.

Now we are ready to fill in the parameters of each node in our network. There are several ways to add these parameters.

bn.cpt(id_c).fillWith([0.4, 0.6]) # remember : c= 0
c
No
Yes
0.40000.6000

Most of the methods using a node id will also work with name of the random variable.

bn.cpt("c").fillWith([0.5, 0.5])
c
No
Yes
0.50000.5000
bn.cpt("s").names
('s', 'c')
bn.cpt("s")[:] = [[0.5, 0.5], [0.9, 0.1]]

Then P(SC=0)=[0.5,0.5]P(S | C=0)=[0.5,0.5]
and P(SC=1)=[0.9,0.1]P(S | C=1)=[0.9,0.1].

print(bn.cpt("s")[1])
[0.9 0.1]

The same process can be performed in several steps:

bn.cpt("s")[0, :] = 0.5 # equivalent to [0.5,0.5]
bn.cpt("s")[1, :] = [0.9, 0.1]
print(bn.cpt("w").names)
bn.cpt("w")
('w', 'r', 's')
w
s
r
No
Yes
No
No
0.89610.1039
Yes
0.54280.4572
Yes
No
0.30540.6946
Yes
0.82550.1745
bn.cpt("w")[0, 0, :] = [1, 0] # r=0,s=0
bn.cpt("w")[0, 1, :] = [0.1, 0.9] # r=0,s=1
bn.cpt("w")[1, 0, :] = [0.1, 0.9] # r=1,s=0
bn.cpt("w")[1, 1, :] = [0.01, 0.99] # r=1,s=1

This is probably the most convenient way:

bn.cpt("w")[{"r": 0, "s": 0}] = [1, 0]
bn.cpt("w")[{"r": 0, "s": 1}] = [0.1, 0.9]
bn.cpt("w")[{"r": 1, "s": 0}] = [0.1, 0.9]
bn.cpt("w")[{"r": 1, "s": 1}] = [0.01, 0.99]
bn.cpt("w")
w
s
r
No
Yes
No
No
1.00000.0000
Yes
0.10000.9000
Yes
No
0.10000.9000
Yes
0.01000.9900

Subscripting with dictionaries is a feature borrowed from OpenBayes. It makes it easier to use and avoids common mistakes that happen when introducing data into the wrong places.

bn.cpt("r")[{"c": 0}] = [0.8, 0.2]
bn.cpt("r")[{"c": 1}] = [0.2, 0.8]

Now our BN is complete. It can be saved in different format :

print(gum.availableBNExts())
bif|dsl|net|bifxml|o3prm|uai|xdsl|pkl
## the preferred format is bifxml
gum.saveBN(bn, "out/WaterSprinkler.bifxml")

But we can save a BN using BIF format

gum.saveBN(bn, "out/WaterSprinkler.bif")
with open("out/WaterSprinkler.bif", "r") as out:
print(out.read())
network "unnamedBN" {
// written by aGrUM 2.3.0.9
}
variable c {
type discrete[2] {No, Yes};
}
variable r {
type discrete[2] {No, Yes};
}
variable w {
type discrete[2] {No, Yes};
}
variable s {
type discrete[2] {No, Yes};
}
probability (c) {
table 0.5 0.5;
}
probability (r | c) {
(No) 0.8 0.2;
(Yes) 0.2 0.8;
}
probability (w | r, s) {
(No, No) 1 0;
(Yes, No) 0.1 0.9;
(No, Yes) 0.1 0.9;
(Yes, Yes) 0.01 0.99;
}
probability (s | c) {
(No) 0.5 0.5;
(Yes) 0.9 0.1;
}
bn2 = gum.loadBN("out/WaterSprinkler.bif")

We can also save and load it in other formats

gum.saveBN(bn, "out/WaterSprinkler.net")
with open("out/WaterSprinkler.net", "r") as out:
print(out.read())
bn3 = gum.loadBN("out/WaterSprinkler.net")
net {
name = unnamedBN;
software = "aGrUM 2.3.0.9";
node_size = (50 50);
}
node c {
states = (No Yes );
label = "c";
ID = "c";
}
node r {
states = (No Yes );
label = "r";
ID = "r";
}
node w {
states = (No Yes );
label = "w";
ID = "w";
}
node s {
states = (No Yes );
label = "s";
ID = "s";
}
potential (c) {
data = ( 0.5 0.5);
}
potential ( r | c ) {
data =
(( 0.8 0.2) % c=No
( 0.2 0.8)); % c=Yes
}
potential ( w | r s ) {
data =
((( 1 0) % s=No r=No
( 0.1 0.9)) % s=Yes r=No
(( 0.1 0.9) % s=No r=Yes
( 0.01 0.99))); % s=Yes r=Yes
}
potential ( s | c ) {
data =
(( 0.5 0.5) % c=No
( 0.9 0.1)); % c=Yes
}
## you can even use pickle
import pickle
with open("out/bn.pkl", "wb") as mypickle:
pickle.dump(bn, file=mypickle)
## or even
gum.saveBN(bn, "out/bn.pkl")

Note: Bayesian network objects in pyAgrum have some specific features that are not supported by standard file formats (e.g., the different types of discrete variables). Only BIFXML and Pickle formats guarantee that all these specific features will be preserved when saving the network.

We have to choose an inference engine to perform calculations for us. Many inference engines are currently available in pyAgrum:

  • Exact inference, particularly :
    • gum.LazyPropagation : an exact inference method that transforms the Bayesian network into a hypergraph called a join tree or a junction tree. This tree is constructed in order to optimize inference computations.
    • others: gum.VariableElimination, gum.ShaferShenoy, …
  • Samplig Inference : approximate inference engine using sampling algorithms to generate a sequence of samples from the joint probability distribution (gum.GibbSSampling, etc.)
  • Loopy Belief Propagation : approximate inference engine using inference algorithm exact for trees but not for DAG
ie = gum.LazyPropagation(bn)
ie.makeInference()
print(ie.posterior("w"))
w |
No |Yes |
---------|---------|
0.3529 | 0.6471 |
from IPython.core.display import HTML
HTML(f"In our BN, $P(W)=${ie.posterior('w')[:]}")

In our BN, P(W)=P(W)=[0.3529 0.6471]

With notebooks, it can be viewed as an HTML table

ie.posterior("w")[:]
array([0.3529, 0.6471])

Suppose now that you know that the sprinkler is on and that it is not cloudy, and you wonder what Is the probability of the grass being wet, i.e., you are interested in distribution P(WS=1,C=0)P(W|S=1,C=0).
The new knowledge you have (sprinkler is on and it is not cloudy) is called evidence. Evidence is entered using a dictionary. When you know precisely the value taken by a random variable, the evidence is called a hard evidence. This is the case, for instance, when I know for sure that the sprinkler is on. In this case, the knowledge is entered in the dictionary as ‘variable name’:label

ie.setEvidence({"s": 0, "c": 0})
ie.makeInference()
ie.posterior("w")
w
No
Yes
0.82000.1800

When you have incomplete knowledge about the value of a random variable, this is called a soft evidence. In this case, this evidence is entered as the belief you have over the possible values that the random variable can take, in other words, as P(evidence|true value of the variable). Imagine for instance that you think that if the sprinkler is off, you have only 50% chances of knowing it, but if it is on, you are sure to know it. Then, your belief about the state of the sprinkler is [0.5, 1] and you should enter this knowledge as shown below. Of course, hard evidence are special cases of soft evidence in which the beliefs over all the values of the random variable but one are equal to 0.

ie.setEvidence({"s": [0.5, 1], "c": [1, 0]})
ie.makeInference()
ie.posterior("w") # using gnb's feature
w
No
Yes
0.32800.6720

the pyagrum.lib.notebook utility proposes certain functions to graphically show distributions.

gnb.showProba(ie.posterior("w"))

svg

gnb.showPosterior(bn, {"s": 1, "c": 0}, "w")

svg

gnb.showInference(bn, evs={})

svg

gnb.showInference(bn, evs={"s": 1, "c": 0})

svg

gnb.showInference(bn, evs={"s": 1, "c": [0.3, 0.9]})

svg

gnb.showInference(bn, evs={"c": [0.3, 0.9]}, targets={"c", "w"})

svg

One of the strength of the Bayesian networks is to form a model that allows to read qualitative knowledge directly from the graph : the conditional independence. aGrUM/pyAgrum comes with a set of tools to query this qualitative knowledge.

## fast create a BN (random paramaters are chosen for the CPTs)
bn = gum.fastBN("A->B<-C->D->E<-F<-A;C->G<-H<-I->J")
bn
G J J A A B B A->B F F A->F C C G G C->G D D C->D C->B H H H->G E E D->E I I I->J I->H F->E

First, one can directly test independence

def testIndep(bn, x, y, knowing):
res = "" if bn.isIndependent(x, y, knowing) else " NOT"
giv = "." if len(knowing) == 0 else f" given {knowing}."
print(f"{x} and {y} are{res} independent{giv}")
testIndep(bn, "A", "C", [])
testIndep(bn, "A", "C", ["E"])
print()
testIndep(bn, "E", "C", [])
testIndep(bn, "E", "C", ["D"])
print()
testIndep(bn, "A", "I", [])
testIndep(bn, "A", "I", ["E"])
testIndep(bn, "A", "I", ["G"])
testIndep(bn, "A", "I", ["E", "G"])
A and C are independent.
A and C are NOT independent given ['E'].
E and C are NOT independent.
E and C are independent given ['D'].
A and I are independent.
A and I are independent given ['E'].
A and I are independent given ['G'].
A and I are NOT independent given ['E', 'G'].

Second, one can investigate the Markov Blanket of a node. The Markov blanket of a node XX is the set of nodes M ⁣B(X)M\!B(X) such that XX is independent from the rest of the nodes given M ⁣B(X)M\!B(X).

print(gum.MarkovBlanket(bn, "C").toDot())
gum.MarkovBlanket(bn, "C")
digraph "no_name" {
node [shape = ellipse];
0[label="A"];
1[label="B"];
2[label="C", color=red];
3[label="D"];
6[label="G"];
7[label="H"];
0 -> 1;
2 -> 3;
2 -> 1;
2 -> 6;
7 -> 6;
}
no_name 0 A 1 B 0->1 2 C 2->1 3 D 2->3 6 G 2->6 7 H 7->6
gum.MarkovBlanket(bn, "J")
no_name 8 I 9 J 8->9

Minimal conditioning set and evidence Impact using probabilistic inference

Section titled “Minimal conditioning set and evidence Impact using probabilistic inference”

For a variable and a list of variables, one can find the sublist that effectively impacts the variable if the list of variables was observed.

[bn.variable(i).name() for i in bn.minimalCondSet("B", ["A", "H", "J"])]
['A']
[bn.variable(i).name() for i in bn.minimalCondSet("B", ["A", "G", "H", "J"])]
['A', 'G', 'H']

This can be also viewed when using gum.LazyPropagation.evidenceImpact(target,evidence) which computes P(targetevidence)P(target|evidence) but reduces as much as possible the set of needed evidence for the result :

ie = gum.LazyPropagation(bn)
ie.evidenceImpact("B", ["A", "C", "H", "G"]) # H,G will be removed w.r.t the minimalCondSet above
B
C
A
0
1
0
0
0.09450.9055
1
0.34540.6546
1
0
0.91020.0898
1
0.75680.2432
ie.evidenceImpact("B", ["A", "G", "H", "J"]) # "J" is not necessary to compute the impact of the evidence
B
A
H
G
0
1
0
0
0
0.83470.1653
1
0.83310.1669
1
0
0.85310.1469
1
0.60700.3930
1
0
0
0.71880.2812
1
0.71800.2820
1
0
0.72800.2720
1
0.60390.3961

PS- the complete code to create the first image

Section titled “PS- the complete code to create the first image”
bn = gum.fastBN("Cloudy?->Sprinkler?->Wet Grass?<-Rain?<-Cloudy?")
bn.cpt("Cloudy?").fillWith([0.5, 0.5])
bn.cpt("Sprinkler?")[:] = [[0.5, 0.5], [0.9, 0.1]]
bn.cpt("Rain?")[{"Cloudy?": 0}] = [0.8, 0.2]
bn.cpt("Rain?")[{"Cloudy?": 1}] = [0.2, 0.8]
bn.cpt("Wet Grass?")[{"Rain?": 0, "Sprinkler?": 0}] = [1, 0]
bn.cpt("Wet Grass?")[{"Rain?": 0, "Sprinkler?": 1}] = [0.1, 0.9]
bn.cpt("Wet Grass?")[{"Rain?": 1, "Sprinkler?": 0}] = [0.1, 0.9]
bn.cpt("Wet Grass?")[{"Rain?": 1, "Sprinkler?": 1}] = [0.01, 0.99]
## the next line control the number of visible digits
gum.config["notebook", "tensor_visible_digits"] = 2
gnb.sideBySide(bn.cpt("Cloudy?"), captions=["$P(Cloudy)$"])
gnb.sideBySide(
bn.cpt("Sprinkler?"),
gnb.getBN(bn, size="3!"),
bn.cpt("Rain?"),
captions=["$P(Sprinkler|Cloudy)$", "", "$P(WetGrass|Sprinkler,Rain)$"],
)
gnb.sideBySide(bn.cpt("Wet Grass?"), captions=["$P(WetGrass|Sprinkler,Rain)$"])
Cloudy?
0
1
0.500.50

$P(Cloudy)$
Sprinkler?
Cloudy?
0
1
0
0.500.50
1
0.900.10

$P(Sprinkler|Cloudy)$
G Wet Grass? Wet Grass? Cloudy? Cloudy? Rain? Rain? Cloudy?->Rain? Sprinkler? Sprinkler? Cloudy?->Sprinkler? Rain?->Wet Grass? Sprinkler?->Wet Grass?
Rain?
Cloudy?
0
1
0
0.800.20
1
0.200.80

$P(WetGrass|Sprinkler,Rain)$
Wet Grass?
Rain?
Sprinkler?
0
1
0
0
1.000.00
1
0.100.90
1
0
0.100.90
1
0.010.99

$P(WetGrass|Sprinkler,Rain)$

(for more, see the notebook : config for pyAgrum)

bn = gum.fastBN("Cloudy?->Sprinkler?->Wet Grass?<-Rain?<-Cloudy?")
bn.cpt("Cloudy?").fillWith([0.5, 0.5])
bn.cpt("Sprinkler?")[:] = [[0.5, 0.5], [0.9, 0.1]]
bn.cpt("Rain?")[{"Cloudy?": 0}] = [0.8, 0.2]
bn.cpt("Rain?")[{"Cloudy?": 1}] = [0.2, 0.8]
bn.cpt("Wet Grass?")[{"Rain?": 0, "Sprinkler?": 0}] = [1, 0]
bn.cpt("Wet Grass?")[{"Rain?": 0, "Sprinkler?": 1}] = [0.1, 0.9]
bn.cpt("Wet Grass?")[{"Rain?": 1, "Sprinkler?": 0}] = [0.1, 0.9]
bn.cpt("Wet Grass?")[{"Rain?": 1, "Sprinkler?": 1}] = [0.01, 0.99]
## the next lines control the visualisation of proba as fraction
gum.config["notebook", "tensor_with_fraction"] = True
gum.config["notebook", "tensor_fraction_with_latex"] = True
gum.config["notebook", "tensor_fraction_limit"] = 100
gnb.sideBySide(bn.cpt("Cloudy?"), captions=["$P(Cloudy)$"])
gnb.sideBySide(
bn.cpt("Sprinkler?"),
gnb.getBN(bn, size="3!"),
bn.cpt("Rain?"),
captions=["$P(Sprinkler|Cloudy)$", "", "$P(WetGrass|Sprinkler,Rain)$"],
)
gnb.sideBySide(bn.cpt("Wet Grass?"), captions=["$P(WetGrass|Sprinkler,Rain)$"])
Cloudy?
0
1

12\frac{1}{2}

12\frac{1}{2}


$P(Cloudy)$
Sprinkler?
Cloudy?
0
1
0

12\frac{1}{2}

12\frac{1}{2}

1

910\frac{9}{10}

110\frac{1}{10}


$P(Sprinkler|Cloudy)$
G Wet Grass? Wet Grass? Cloudy? Cloudy? Rain? Rain? Cloudy?->Rain? Sprinkler? Sprinkler? Cloudy?->Sprinkler? Rain?->Wet Grass? Sprinkler?->Wet Grass?
Rain?
Cloudy?
0
1
0

45\frac{4}{5}

15\frac{1}{5}

1

15\frac{1}{5}

45\frac{4}{5}


$P(WetGrass|Sprinkler,Rain)$
Wet Grass?
Rain?
Sprinkler?
0
1
0
0

11

00

1

110\frac{1}{10}

910\frac{9}{10}

1
0

110\frac{1}{10}

910\frac{9}{10}

1

1100\frac{1}{100}

99100\frac{99}{100}


$P(WetGrass|Sprinkler,Rain)$