Conditional Linear Gaussian models

import pyagrum as gum
import pyagrum.lib.notebook as gnb
import pyagrum.lib.bn_vs_bn as gcm

import pyagrum.clg as gclg
import pyagrum.clg.notebook as gclgnb

Build a CLG model

From scratch

Suppose we want to build a CLG with these specifications $A={\cal N}(5,1)$, $B={\cal N}(4,3)$ and $C=2.A+3.B+{\cal N}(3,2)$

model = gclg.CLG()
model.add(gclg.GaussianVariable("A", 5, 1))
model.add(gclg.GaussianVariable("C", 3, 2))
model.add(gclg.GaussianVariable("B", 4, 3))
model.addArc("A", "C", 2)
model.addArc("B", "C", 3)
model

From SEM (Structural Equation Model)

We can create a Conditional Linear Gaussian Bayesian networ(CLG model) using a SEM-like syntax.

A = 4.5 [0.3] means that the mean of the distribution for Gaussian random variable A is 4.5 and ist standard deviation is 0.3.

B = 3 + 0.8F [0.3] means that the mean of the distribution for the Gaussian random variable B is 3 and the standard deviation is 0.3.

pyagrum.CLG.SEM is a set of static methods to manipulate this kind of SEM.

sem2 = """
A=4.5 [0.3] # comments are allowed
F=7 [0.5]
B=3 + 1.2F [0.3]
C=9 +  2A + 1.5B [0.6]
D=9 + C + F[0.7]
E=9 + D [0.9]
"""

model2 = gclg.SEM.toclg(sem2)

gnb.show(model2)

One can of course build the SEM from a CLG using pyagrum.CLG.SEM.tosem :

gnb.flow.row(
  model,
  "<pre><div align='left'>" + gclg.SEM.tosem(model) + "</div></pre>",
  captions=["the first CLG model", "the SEM from the CLG"],
)

the first CLG model

B=4[3]
A=5[1]
C=3+2A+3B[2]

the SEM from the CLG

And this SEM allows of course input/output format for CLG

gclg.SEM.saveCLG(model2, "out/model2.sem")

print("=== file content ===")
with open("out/model2.sem", "r") as file:
  for line in file.readlines():
    print(line, end="")
print("====================")

=== file content ===
F=7.0[0.5]
B=3.0+1.2F[0.3]
A=4.5[0.3]
C=9.0+2.0A+1.5B[0.6]
D=9.0+F+C[0.7]
E=9.0+D[0.9]
====================

model3 = gclg.SEM.loadCLG("out/model2.sem")
gnb.sideBySide(model2, model3, captions=["saved model", "loaded model"])

saved model

loaded model

input/output with pickle for CLG

import pickle

with open("out/testCLG.pkl", "bw") as f:
  pickle.dump(model3, f)
model3

model.dag().sizeArcs()

2

with open("out/testCLG.pkl", "br") as f:
  copyModel3 = pickle.load(f)
copyModel3

Exact or approximated Inference

Exact inference : Variable Elimination

Compute some posterior using difference exact inference

ie = gclg.CLGVariableElimination(model2)
ie.updateEvidence({"D": 3})

print(ie.posterior("A"))
print(ie.posterior("B"))
print(ie.posterior("C"))
print(ie.posterior("D"))
print(ie.posterior("E"))
print(ie.posterior("F"))


v = ie.posterior("E")
print(v)
print(f"  - mean(E|D=3)={v.mu()}")
print(f"  - stdev(E|D=3)={v.sigma()}")

A:1.9327650111193468[0.28353638852446156]
B:-2.5058561897702[0.41002992170553515]
C:3.9722757598220895[0.5657771474513671]
D:3[0]
E:12.0[0.9]
F:-2.9836916234247597[0.32358490464094586]
E:12.0[0.9]
  - mean(E|D=3)=12.0
  - stdev(E|D=3)=0.9

gnb.sideBySide(
  model2,
  gclgnb.getInference(model2, evs={"D": 3}, size="3!"),
  gclgnb.getInference(model2, evs={"D": 3, "F": 1}),
  captions=["The CLG", "First inference", "Second inference"],
)

The CLG

First inference

Second inference

Approximated inference : MonteCarlo Sampling

When the model is too complex for exact infernece, we can use forward sampling to generate 5000 samples from the original CLG model.

fs = gclg.ForwardSampling(model2)
fs.makeSample(5000).tocsv("./out/model2.csv")

We will use the generated database to do learning. But before, we can also compute posterior but without evidence :

ie = gclg.CLGVariableElimination(model2)
print("| 'Exact' inference                        | Results from sampling                    |")
print("|------------------------------------------|------------------------------------------|")
for i in model2.names():
  print(f"| {str(ie.posterior(i)):40} | {str(gclg.GaussianVariable(i, fs.mean_sample(i), fs.stddev_sample(i))):40} |")

| 'Exact' inference                        | Results from sampling                    |
|------------------------------------------|------------------------------------------|
| A:4.499999999999998[0.3]                 | A:4.503776111579754[0.29947675934344437] |
| F:7.000000000000008[0.5000000000000002]  | F:6.987746094829205[0.5013179465519744]  |
| B:11.399999999999999[0.6708203932499367] | B:11.38781277159795[0.6623799495758246]  |
| C:35.099999999999994[1.3162446581088183] | C:35.10464359617556[1.3095831016940265]  |
| D:51.10000000000002[1.8364367672206963]  | D:51.095154564791166[1.8206151110413087] |
| E:60.100000000000016[2.0451161336217565] | E:60.096823722638106[2.020983807559692]  |

Now with the generated database and the original model, we can calculate the log-likelihood of the model.

print("log-likelihood w.r.t orignal model : ", model2.logLikelihood("./out/model2.csv"))

log-likelihood w.r.t orignal model :  -22019.357069680933

Learning a CLG from data

Use the generated database to do our RAvel Learning. This part needs some time to run.

## RAveL learning
learner = gclg.CLGLearner("./out/model2.csv")

We can get the learned_clg model with function learn_clg() which contains structure learning and parameter estimation.

learned_clg = learner.learnCLG()
gnb.sideBySide(model2, learned_clg, captions=["original CLG", "learned CLG"])

original CLG

learned CLG

Compare the learned model’s structure with that of the original model’.

cmp = gcm.GraphicalBNComparator(model2.asDiscreteBN(), learned_clg.asDiscreteBN())
print(f"F-score(original_clg,learned_clg) : {cmp.scores()['fscore']}")

F-score(original_clg,learned_clg) : 0.8

Get the learned model’s parameters and compare them with the original model’s parameters using the SEM syntax.

gnb.flow.row(
  "<pre><div align='left'>" + gclg.SEM.tosem(model2) + "</div></pre>",
  "<pre><div align='left'>" + gclg.SEM.tosem(learned_clg) + "</div></pre>",
  captions=["original sem", "learned sem"],
)

F=7.0[0.5]
B=3.0+1.2F[0.3]
A=4.5[0.3]
C=9.0+2.0A+1.5B[0.6]
D=9.0+F+C[0.7]
E=9.0+D[0.9]

original sem

E=60.096823722638106[2.020983807559692]
F=6.987746094829205[0.5013179465519744]
B=3.104657780785809+1.185382937273796F[0.2925913363477699]
D=5.825201280360922+1.063200862334291F+0.6296601594999158E[0.7033446428058315]
A=4.503776111579754[0.29947675934344437]
C=2.0509646739427296+1.2237610651514252A+0.6314112877840555B+0.3983105583477016D[0.46970348750580965]

learned sem

We can algo do parameter estimation only with function fitParameters() if we already have the structure of the model.

## We can copy the original CLG
copy_original = gclg.CLG(model2)

## RAveL learning again
RAveL_l = gclg.CLGLearner("./out/model2.csv")

## Fit the parameters of the copy clg
RAveL_l.fitParameters(copy_original)

copy_original

Compare two CLG models

We first create two CLG from two SEMs.

## TWO DIFFERENT CLGs

## FIRST CLG
clg1 = gclg.SEM.toclg("""
## hyper parameters
A=4[1]
B=3[5]
C=-2[5]

#equations
D=A[.2] # D is a noisy version of A
E=1+D+2B [2]
F=E+C+B+E [0.001]
""")

## SECOND CLG
clg2 = gclg.SEM.toclg("""
## hyper parameters
A=4[1]
B=3+A[5]
C=-2+2B+A[5]

#equations
D=A[.2] # D is a noisy version of A
E=1+D+2B [2]
F=E+C [0.001]
""")

This cell shows how to have a quick view of the differences

gnb.flow.row(clg1, clg2, gcm.graphDiff(clg1, clg2), gcm.graphDiffLegend(), gcm.graphDiff(clg2, clg1))

We compare the CLG models.

## We use the F-score to compare the two CLGs
cmp = gcm.GraphicalBNComparator(clg1.asDiscreteBN(), clg1.asDiscreteBN())
print(f"F-score(clg1,clg1) : {cmp.scores()['fscore']}")

cmp = gcm.GraphicalBNComparator(clg1.asDiscreteBN(), clg2.asDiscreteBN())
print(f"F-score(clg1,clg2) : {cmp.scores()['fscore']}")

F-score(clg1,clg1) : 1.0
F-score(clg1,clg2) : 0.7142857142857143

## The complete list of structural scores is :
print("score(clg1,clg2) :")
for score, val in cmp.scores().items():
  print(f"  - {score} : {val}")

score(clg1,clg2) :
  - count : {'tp': 5, 'tn': 6, 'fp': 3, 'fn': 1}
  - recall : 0.8333333333333334
  - precision : 0.625
  - fscore : 0.7142857142857143
  - dist2opt : 0.41036907507483766
  - sid : 3

Forward Sampling

## We create a simple CLG with 3 variables
clg = gclg.CLG()
## prog=« sigma=2;X=N(5);Y=N(3);Z=X+Y »
A = gclg.GaussianVariable(mu=2, sigma=1, name="A")
B = gclg.GaussianVariable(mu=1, sigma=2, name="B")
C = gclg.GaussianVariable(mu=2, sigma=3, name="C")

idA = clg.add(A)
idB = clg.add(B)
idC = clg.add(C)

clg.addArc(idA, idB, 1.5)
clg.addArc(idB, idC, 0.75)

## We can show it as a graph
original_clg = gclgnb.CLG2dot(clg)
original_clg

fs = gclg.ForwardSampling(clg)
fs.makeSample(10)

<pyagrum.clg.forwardSampling.ForwardSampling at 0x117dd5450>

print("A's sample_variance: ", fs.variance_sample(0))
print("B's sample_variance: ", fs.variance_sample("B"))
print("C's sample_variance: ", fs.variance_sample(2))

A's sample_variance:  2.162780323848712
B's sample_variance:  7.042666325597042
C's sample_variance:  12.627122883187766

print("A's sample_mean: ", fs.mean_sample("A"))
print("B's sample_mean: ", fs.mean_sample("B"))
print("C's sample_mean: ", fs.mean_sample("C"))

A's sample_mean:  2.2697860317737217
B's sample_mean:  4.712223796042158
C's sample_mean:  6.521127871421228

fs.toarray()

array([[ 5.43081542,  8.78960455,  6.26537806],
       [ 1.1723556 ,  0.50342667,  5.43421933],
       [ 1.02312139,  2.41771592,  4.89042683],
       [ 2.18967057,  6.63423781, 11.6027913 ],
       [ 4.50194765,  7.53038384, 13.46073451],
       [ 1.60924503,  7.07449439,  8.62891688],
       [ 0.93675031,  1.92544683,  3.3538191 ],
       [ 1.1766549 ,  2.5074888 ,  6.24443528],
       [ 2.81498398,  4.18787959,  4.09662382],
       [ 1.84231548,  5.55155957,  1.2339336 ]])

## export to dataframe
fs.topandas()

	A	B	C
0	5.430815	8.789605	6.265378
1	1.172356	0.503427	5.434219
2	1.023121	2.417716	4.890427
3	2.189671	6.634238	11.602791
4	4.501948	7.530384	13.460735
5	1.609245	7.074494	8.628917
6	0.936750	1.925447	3.353819
7	1.176655	2.507489	6.244435
8	2.814984	4.187880	4.096624
9	1.842315	5.551560	1.233934

## export to csv
fs.makeSample(10000)
fs.tocsv("./out/samples.csv")

PC-algorithm & Parameter Estimation

The module allows to investigale more deeply into the learning algorithm.

We first create a random CLG model with 5 variables.

## Create a new random CLG
clg = gclg.randomCLG(nb_variables=5, names="ABCDE")

## Display the CLG
print(clg)

A=4.238030191922599[6.823820213297951]
B=-4.132995496588267+7.547435155089987A[9.416753814775685]
C=-2.131132877989069+6.802147141657968B[7.362781603908117]
D=-0.7130383694753473-7.502000062332959C[3.7070850282956638]
E=2.0366077881169273+8.608799187031885A[2.0950479734141836]

We then do the Forward Sampling and CLGLearner.

n = 20  # n is the selected values of MC number n in n-MCERA
K = 10000  # K is the list of selected values of number of samples
Delta = 0.05  # Delta is the FWER we want to control

## Sample generation
fs = gclg.ForwardSampling(clg)
fs.makeSample(K).tocsv("./out/clg.csv")

## Learning
RAveL_l = gclg.CLGLearner("./out/clg.csv", n_sample=n, fwer_delta=Delta)

We use the PC algorithme to learn the structure of the model.

## Use the PC algorithm to get the skeleton
C = RAveL_l.PC_algorithm(order=clg.nodes(), verbose=False)
print("The final skeleton is:\n", C)

The final skeleton is:
 {0: set(), 1: {3}, 2: set(), 3: set(), 4: {0}}

## Create a Mixedgraph to display the skeleton
RAveL_MixGraph = gum.MixedGraph()

## Add variables
for i in range(len(clg.names())):
  RAveL_MixGraph.addNodeWithId(i)
  RAveL_MixGraph.setName(i, clg.variable(i).name())

## Add arcs and edges
for father, kids in C.items():
  for kid in kids:
    if father in C[kid]:
      RAveL_MixGraph.addEdge(father, kid)
    else:
      RAveL_MixGraph.addArc(father, kid)

RAveL_MixGraph

## Create a BN with the same structure as the CLG
bn = clg.asDiscreteBN()

## Compare the result above with the EssentialGraph
Real_EssentialGraph = gum.EssentialGraph(bn)

Real_EssentialGraph

## create a CLG from the skeleton of PC algorithm
clg_PC = gclg.CLG()
for node in clg.nodes():
  clg_PC.add(clg.variable(node))
for father, kids in C.items():
  for kid in kids:
    clg_PC.addArc(father, kid)

## Compare the structure of the created CLG and the original CLG
print(f"F-score : {clg.structuralFScore(clg_PC)}")

F-score : 0.6666666666666666

We can also do the parameter learning.

id2mu, id2sigma, arc2coef = RAveL_l.estimate_parameters(C)

for node in clg.nodes():
  print(f"Real Value: node {node} : mu = {clg.variable(node)._mu}, sigma = {clg.variable(node)._sigma}")
  print(f"Estimation: node {node} : mu = {id2mu[node]}, sigma = {id2sigma[node]}")


for arc in clg.arcs():
  print(f"Real Value: arc {arc} : coef = {clg.coefArc(*arc)}")
  print(f"Estimation: arc {arc} : coef = {(arc2coef[arc] if arc in arc2coef else '-')}")

Real Value: node 0 : mu = 2.0366077881169273, sigma = 2.0950479734141836
Estimation: node 0 : mu = 2.054266563309085, sigma = 2.123011948193108
Real Value: node 1 : mu = -2.131132877989069, sigma = 7.362781603908117
Estimation: node 1 : mu = 189.41549903074008, sigma = 354.93067649834893
Real Value: node 2 : mu = -4.132995496588267, sigma = 9.416753814775685
Estimation: node 2 : mu = 28.163831301604223, sigma = 52.16303364725259
Real Value: node 3 : mu = -0.7130383694753473, sigma = 3.7070850282956638
Estimation: node 3 : mu = -0.682196280556127, sigma = 3.7411505338945514
Real Value: node 4 : mu = 4.238030191922599, sigma = 6.823820213297951
Estimation: node 4 : mu = 4.2773930915293255, sigma = 6.802251293026495
Real Value: arc (4, 0) : coef = 8.608799187031885
Estimation: arc (4, 0) : coef = 8.609129056137842
Real Value: arc (1, 3) : coef = -7.502000062332959
Estimation: arc (1, 3) : coef = -7.5020983587072845
Real Value: arc (2, 1) : coef = 6.802147141657968
Estimation: arc (2, 1) : coef = -
Real Value: arc (4, 2) : coef = 7.547435155089987
Estimation: arc (4, 2) : coef = -