Conditional Linear Gaussian models

Creative Commons License

aGrUM

interactive online version

In [1]:
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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 6
      3 import pyAgrum.lib.bn_vs_bn as gcm
      5 import pyAgrum.clg as gclg
----> 6 import pyAgrum.clg.notebook as gclgnb

File ~/.virtualenvs/devAgrum/lib/python3.13/site-packages/pyAgrum/clg/notebook.py:32
     29 import pyAgrum.lib.image as gimg
     31 import pyAgrum.clg.CLG as CLG
---> 32 import pyAgrum.clg.inference as clginference
     33 import pyAgrum.lib.notebook as gnb
     36 def CLG2dot(clg, *, size=None, nodeColor=None, arcWidth=None, arcColor=None, cmapNode=None, cmapArc=None, showMsg=None):

ModuleNotFoundError: No module named 'pyAgrum.clg.inference'

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)\)

In [ ]:
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.

In [ ]:
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)
In [ ]:
gnb.show(model2)

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

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

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

In [ ]:
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("====================")
In [ ]:
model3=gclg.SEM.loadCLG("out/model2.sem")
gnb.sideBySide(model2,model3,captions=["saved model","loaded model"])

input/output with pickle

In [ ]:
import pickle
with open("test.pkl","bw") as f:
  pickle.dump(model3,f)
model3
In [ ]:
with open("test.pkl","br") as f:
  copyModel3=pickle.load(f)
copyModel3

Exact or approximated Inference

Exact inference : Variable Elimination

Compute some posterior using difference exact inference

In [ ]:
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()}")
In [ ]:
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"])

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.

In [ ]:
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 :

In [ ]:
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} |")

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

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

Learning a CLG from data

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

In [ ]:
# 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.

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

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

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

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

In [ ]:
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"])

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

In [ ]:
# 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.

In [ ]:
# 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

In [ ]:
print(gum.config)
In [ ]:
gnb.flow.row(clg1,clg2,gcm.graphDiff(clg1,clg2),
             gcm.graphDiffLegend(),
             gcm.graphDiff(clg2,clg1))

We compare the CLG models.

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

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

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

Forward Sampling

In [ ]:
# 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
In [ ]:
fs = gclg.ForwardSampling(clg)
fs.makeSample(10)
In [ ]:
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))
In [ ]:
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'))
In [ ]:
fs.toarray()
In [ ]:
# export to dataframe
fs.topandas()
In [ ]:
# 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.

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

# Display the CLG
print(clg)

We then do the Forward Sampling and CLGLearner.

In [ ]:
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.

In [ ]:
# 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)
In [ ]:
# Create a Mixedgraph to display the skeleton
RAveL_MixGraph = gum.MixedGraph()

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

# 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
In [ ]:
# Create a BN with the same structure as the CLG
bn = gum.BayesNet()
# add variables
for name in clg.names():
  new_variable = gum.LabelizedVariable(name,'a labelized variable',2)
  bn.add(new_variable)
# add arcs
for arc in clg.arcs():
  bn.addArc(arc[0], arc[1])

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

Real_EssentialGraph
In [ ]:
# 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.CompareStructure(clg_PC)}")

We can also do the parameter learning.

In [ ]:
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 '-')}")
In [ ]: