Quasi-continuous BN

aGrUM cannot (currently) deal with with continuous variables. However, a discrete variable with a large enough domain size is an approximation of such variables.

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

# nbr of states for quasi continuous variables. You can change the value
# but be careful of the quadratic behavior of both memory and time complexity
# in this example.
minB, maxB = -3, 3
minC, maxC = 4, 14
NB = 200

# the line with fastBN replace the commented ones.
# bn=gum.BayesNet()
# bn.add(gum.LabelizedVariable("A","A binary variable",2))
# bn.add(gum.NumericalDiscreteVariable("B","A range variable",minB,maxB,NB))
# bn.addArc("A","B")
bn = gum.fastBN(f"A[3]->B[{minB}:{maxB}:{NB}]")

gnb.showBN(bn)

bn.cpt("A")[:] = [0.4, 0.1, 0.5]
gnb.showProba(bn.cpt("A"))

CPT for quasi-continuous variables (with parents)

Using python (and scipy), it is easy to find pdf for continuous variable

from scipy.stats import norm

bn.cpt("B").fillFromDistribution(norm, loc="-2+A*2", scale="(5+A*4)/20")

gnb.flow.clear()
gnb.flow.add(gnb.getProba(bn.cpt("B").extract({"A": 0})), caption="P(B|A=0)")
gnb.flow.add(gnb.getProba(bn.cpt("B").extract({"A": 1})), caption="P(B|A=1)")
gnb.flow.add(gnb.getProba(bn.cpt("B").extract({"A": 2})), caption="P(B|A=1)")
gnb.flow.display()

P(B|A=0)

P(B|A=1)

P(B|A=1)

Quasi-continuous inference (with no evidence)

gnb.showPosterior(bn, target="B", evs={})
gnb.showInference(bn)

Quasi-continuous inference with numerical evidence expressed as logical propositions evEq or evIn,evLt, evGt; and boolean operators

gnb.showInference(bn, evs=[bn.evIn("B", -1, 2)])  # we observed B between -1 and 2

gnb.showInference(bn, evs=[~bn.evIn("B", -1, 0)])  # we observed B not being between -1 and 0

gnb.showInference(bn, evs=[bn.evLt("B", 1)])  # we observed B being less than 1

gnb.showInference(bn, evs=[bn.evEq("B", 0) | bn.evEq("B", -2)])  # we observed B being -1 or 2

Quasi-continuous variable with quasi-continuous parent

bn = gum.fastBN("A[3]->B->C", f"[{minB}:{maxB}:{NB}]")  # default type of variables (for B and C)
gnb.showBN(bn)  # B and C are quasi-continouous

Even if this BN is quite small (and linear), the size of nodes \(B\) et \(C\) are rather big and creates a complex model (NBxNB parameters in \(P(C|B)\)).

print("nombre de paramètres du bn : {0}".format(bn.dim()))
print("domaine du bn : 10^{0}".format(bn.log10DomainSize()))

nombre de paramètres du bn : 40399
domaine du bn : 10^5.079181246047625

from scipy.stats import gamma

bn.cpt("B").fillFromDistribution(norm, loc="-2+A*2", scale="(5+A*4)/20")
bn.cpt("C").fillFromDistribution(gamma, a="B+3.1", loc=-3, scale=5)


def showCgivenBequals(x: float):
  gnb.flow.add(gnb.getProba(bn.cpt("C").extract({"B": f"{x}"})), caption=f"P(C|B={x})")


gnb.flow.clear()
showCgivenBequals(0)
showCgivenBequals(3)
showCgivenBequals(-3)
# showB(NB-1)
gnb.flow.display()

P(C|B=0)

P(C|B=3)

P(C|B=-3)

Inference in quasi-continuous BN

import time

ts = time.time()
ie = gum.LazyPropagation(bn)
ie.makeInference()
q = ie.posterior("C")
te = time.time()
gnb.flow.add(
  gnb.getPosterior(bn, target="C", evs={}),
  caption=f"P(C) computed in {te - ts:2.5f} sec for a model with {bn.dim()} parameters",
)
gnb.flow.display()

P(C) computed in 0.00067 sec for a model with 40399 parameters

Changing prior

gnb.showInference(bn, size="10")

bn.cpt("A")[:] = [0.9, 0.1, 0.0]

gnb.showInference(bn, size="10")

inference with evidence in quasi-continuous BN

We want to compute

• \(P(A | C=3)\)

• \(P(B | C=3)\)

ie = gum.LazyPropagation(bn)
ie.setEvidence([bn.evEq("C", 3)])
ie.makeInference()
gnb.showProba(ie.posterior("B"))

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

gnb.showInference(bn, evs=[bn.evEq("C", 3)])

Multiple inference : MAP DECISION between complex distributions

What is the behaviour of \(P(A | C=i)\)   when \(i\) varies ? I.e. we perform a MAP decision between the two models (\(A=0\)  for the Gaussian distribution and \(A=1\)  for the generalized hyperbolic distribution).

import matplotlib.pyplot as plt
import numpy as np

bn.cpt("A")[:] = [0.1, 0.7, 0.2]
ie = gum.LazyPropagation(bn)
p0 = []
p1 = []
p2 = []
x = bn.variable("C").ticks()
for i in x:
  ie.setEvidence([bn.evEq("C", i)])
  ie.makeInference()
  p0.append(ie.posterior("A")[0])
  p1.append(ie.posterior("A")[1])
  p2.append(ie.posterior("A")[2])

plt.plot(x, p0)
plt.plot(x, p1)
plt.plot(x, p2)
plt.title(f"P( A | C=x) with prior p(A)={bn.cpt('A').tolist()}")
plt.legend(["A=0", "A=1", "A=2"], loc="best")
inters = (np.transpose(p0) > np.transpose(p1)).argmin()

plt.text(
  x[inters] + 0.2,
  p0[inters],
  "{0:5.4},{1:5.4f}  ".format(x[inters], p0[inters]),
  bbox=dict(facecolor="red", alpha=0.1),
  ha="left",
)
plt.show()
print("\n\n")
print("==========================================================")
print(f"  DECISION RULE : If C<{x[inters]:0.3f} Then A=0 else A=1")
print("==========================================================")

==========================================================
  DECISION RULE : If C<-1.950 Then A=0 else A=1
==========================================================

Same MAP with another \(P(A)\)

bn.cpt("A").fillWith([0.4, 0.3, 0.3])
ie = gum.LazyPropagation(bn)
p0 = []
p1 = []
p2 = []
x = bn.variable("C").ticks()
for i in x:
  ie.setEvidence([bn.evEq("C", i)])
  ie.makeInference()
  p0.append(ie.posterior("A")[0])
  p1.append(ie.posterior("A")[1])
  p2.append(ie.posterior("A")[2])

plt.plot(x, p0)
plt.plot(x, p1)
plt.plot(x, p2)
plt.title(f"P( A | C=x) with prior p(A)={bn.cpt('A').tolist()}")
plt.legend(["A=0", "A=1", "A=2"], loc="best")
inters1 = (np.transpose(p0) > np.transpose(p1)).argmin()
inters2 = (np.transpose(p1) > np.transpose(p2)).argmin()

plt.text(
  x[inters1] - 0.2,
  p0[inters1],
  "{0:5.3f},{1:5.4f}  ".format(x[inters1], p0[inters1]),
  bbox=dict(facecolor="red", alpha=0.1),
  ha="right",
)
plt.text(
  x[inters2] + 0.2,
  p1[inters2],
  "{0:5.3f},{1:5.4f}  ".format(x[inters2], p0[inters2]),
  bbox=dict(facecolor="red", alpha=0.1),
  ha="left",
)
plt.show()
print("\n\n")
print("==========================================================")
print(f"  DECISION RULE : If C<{x[inters1]:0.3f} Then A=0")
print(f"                  ElseIf C<{x[inters2]:0.3f} Then A=1")
print("                  Else A=2")
print("==========================================================")

==========================================================
  DECISION RULE : If C<0.540 Then A=0
                  ElseIf C<1.410 Then A=1
                  Else A=2
==========================================================