Potentials : named tensors¶

In pyAgrum, Potentials represent multi-dimensionnal arrays with (discrete) random variables attached to each dimension. This mathematical object have tensorial operators w.r.t. to the variables attached.

In [1]:

import pyAgrum as gum
import pyAgrum.lib.notebook as gnb

a,b,c=[gum.LabelizedVariable(s,s,2) for s in "abc"]

potential algebra¶

In [2]:

p1=gum.Potential().add(a).add(b).fillWith([1,2,3,4]).normalize()
p2=gum.Potential().add(b).add(c).fillWith([4,5,2,3]).normalize()

In [3]:

gnb.flow.row(p1,p2,p1+p2,
              captions=['p1','p2','p1+p2'])

	a
b	0	1
0	0.1000	0.2000
1	0.3000	0.4000

p1

	b
c	0	1
0	0.2857	0.3571
1	0.1429	0.2143

p2

		b
a	c	0	1
0	0	0.3857	0.6571
0	1	0.2429	0.5143
1	0	0.4857	0.7571
1	1	0.3429	0.6143

p1+p2

In [4]:

p3=p1+p2 
gnb.showPotential(p3/p3.margSumOut(["b"]))

		c
b	a	0	1
0	0	0.3699	0.3208
0	1	0.3908	0.3582
1	0	0.6301	0.6792
1	1	0.6092	0.6418

In [5]:

p4=gum.Potential()+p3
gnb.flow.row(p3,p4,
              captions=['p3','p4'])

		b
a	c	0	1
0	0	0.3857	0.6571
0	1	0.2429	0.5143
1	0	0.4857	0.7571
1	1	0.3429	0.6143

p3

		b
a	c	0	1
0	0	1.3857	1.6571
0	1	1.2429	1.5143
1	0	1.4857	1.7571
1	1	1.3429	1.6143

p4

Bayes' theorem¶

In [6]:

bn=gum.fastBN("a->c;b->c",3)
bn

Out[6]:

In such a small bayes net, we can directly manipulate $P(a,b,c)$. For instance : $$P(b|c)=\frac{\sum_{a} P(a,b,c)}{\sum_{a,b} P(a,b,c)}$$

In [7]:

pABC=bn.cpt("a")*bn.cpt("b")*bn.cpt("c")
pBgivenC=(pABC.margSumOut(["a"])/pABC.margSumOut(["a","b"]))

pBgivenC.putFirst("b") # in order to have b horizontally in the table

Out[7]:

	b
c	0	1	2
0	0.3267	0.1174	0.5559
1	0.4532	0.1154	0.4314
2	0.4143	0.2231	0.3626

Joint, marginal probability, likelihood¶

Let's compute the joint probability $P(A,B)$ from $P(A,B,C)$

In [8]:

pAC=pABC.margSumOut(["b"])
print("pAC really is a probability : it sums to {}".format(pAC.sum()))
pAC

pAC really is a probability : it sums to 1.0000000000000002

Out[8]:

	a
c	0	1	2
0	0.0156	0.0621	0.2078
1	0.0206	0.1758	0.1153
2	0.0337	0.1927	0.1765

Computing $p(A)$¶

In [9]:

pAC.margSumOut(["c"])

Out[9]:

a
0	1	2
0.0699	0.4306	0.4996

Computing $p(A |C=1)$¶

It is easy to compute $p(A, C=1)$

In [10]:

pAC.extract({"c":1})

Out[10]:

a
0	1	2
0.0206	0.1758	0.1153

Moreover, we know that $P(C=1)=\sum_A P(A,C=1)$

In [11]:

pAC.extract({"c":1}).sum()

Out[11]:

0.3116561031542595

Now we can compute $p(A|C=1)=\frac{P(A,C=1)}{p(C=1)}$

In [12]:

pAC.extract({"c":1}).normalize()

Out[12]:

a
0	1	2
0.0660	0.5640	0.3700

Computing $P(A|C)$¶

$P(A|C)$ is represented by a matrix that verifies $p(A|C)=\frac{P(A,C)}{P(C}$

In [13]:

pAgivenC=(pAC/pAC.margSumIn("c")).putFirst("a") 
# putFirst("a") : to correctly show a cpt, the first variable have to bethe conditionned one
gnb.flow.row(pAgivenC,pAgivenC.extract({'c':1}),
               captions=["$P(A|C)$","$P(A|C=1)$"])

	a
c	0	1	2
0	0.0547	0.2175	0.7278
1	0.0660	0.5640	0.3700
2	0.0836	0.4783	0.4381

$P(A|C)$

a
0	1	2
0.0660	0.5640	0.3700

$P(A|C=1)$

Likelihood $P(A=2|C)$¶

A likelihood can also be found in this matrix.

In [14]:

pAgivenC.extract({'a':2})

Out[14]:

c
0	1	2
0.7278	0.3700	0.4381

A likelihood does not have to sum to 1. It is not relevant to normalize it.

In [15]:

pAgivenC.margSumIn(["a"])

Out[15]:

a
0	1	2
0.2043	1.2598	1.5359

entropy of potential¶

In [16]:

%matplotlib inline
from pylab import *
import matplotlib.pyplot as plt
import numpy as np

In [17]:

p1=gum.Potential().add(a)
x = np.linspace(0, 1, 100)
plt.plot(x,[p1.fillWith([p,1-p]).entropy() for p in x])
plt.show()

In [18]:

t=gum.LabelizedVariable('t','t',3)
p1=gum.Potential().add(t)

def entrop(bc):
    """
    bc is a list [a,b,c] close to a distribution 
    (normalized just to be sure)
    """
    return p1.fillWith(bc).normalize().entropy()

import matplotlib.tri as tri

corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])

# Mid-points of triangle sides opposite of each corner
midpoints = [(corners[(i + 1) % 3] + corners[(i + 2) % 3]) / 2.0 \
             for i in range(3)]
def xy2bc(xy, tol=1.e-3):
    """
    From 2D Cartesian coordinates to barycentric.
    """
    s = [(corners[i] - midpoints[i]).dot(xy - midpoints[i]) / 0.75 \
         for i in range(3)]
    return np.clip(s, tol, 1.0 - tol)
    
def draw_entropy(nlevels=200, subdiv=6, **kwargs):
    import math

    refiner = tri.UniformTriRefiner(triangle)
    trimesh = refiner.refine_triangulation(subdiv=subdiv)
    pvals = [entrop(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]

    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)
    plt.axis('equal')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    plt.axis('off')
    
draw_entropy()
plt.show()