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'])
|
| |
---|---|---|
0.1000 | 0.2000 | |
0.3000 | 0.4000 |
|
| |
---|---|---|
0.2857 | 0.3571 | |
0.1429 | 0.2143 |
|
| ||
---|---|---|---|
| 0.3857 | 0.6571 | |
0.2429 | 0.5143 | ||
| 0.4857 | 0.7571 | |
0.3429 | 0.6143 |
In [4]:
p3=p1+p2
gnb.showPotential(p3/p3.sumOut(["b"]))
|
| ||
---|---|---|---|
| 0.3699 | 0.3208 | |
0.3908 | 0.3582 | ||
| 0.6301 | 0.6792 | |
0.6092 | 0.6418 |
In [5]:
p4=gum.Potential()+p3
gnb.flow.row(p3,p4,
captions=['p3','p4'])
|
| ||
---|---|---|---|
| 0.3857 | 0.6571 | |
0.2429 | 0.5143 | ||
| 0.4857 | 0.7571 | |
0.3429 | 0.6143 |
|
| ||
---|---|---|---|
| 1.3857 | 1.6571 | |
1.2429 | 1.5143 | ||
| 1.4857 | 1.7571 | |
1.3429 | 1.6143 |
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.sumOut(["a"])/pABC.sumOut(["a","b"]))
pBgivenC.putFirst("b") # in order to have b horizontally in the table
Out[7]:
|
|
| |
---|---|---|---|
0.0022 | 0.9786 | 0.0192 | |
0.0058 | 0.9733 | 0.0209 | |
0.0038 | 0.9752 | 0.0210 |
Joint, marginal probability, likelihood¶
Let's compute the joint probability $P(A,B)$ from $P(A,B,C)$
In [8]:
pAC=pABC.sumOut(["b"])
print("pAC really is a probability : it sums to {}".format(pAC.sum()))
pAC
pAC really is a probability : it sums to 0.9999999999999998
Out[8]:
|
|
| |
---|---|---|---|
0.1422 | 0.1240 | 0.0582 | |
0.1299 | 0.0704 | 0.1512 | |
0.0666 | 0.1115 | 0.1459 |
Computing $p(A)$¶
In [9]:
pAC.sumOut(["c"])
Out[9]:
|
|
|
---|---|---|
0.3387 | 0.3060 | 0.3553 |
Computing $p(A |C=1)$¶
It is easy to compute $p(A, C=1)$
In [10]:
pAC.extract({"c":1})
Out[10]:
|
|
|
---|---|---|
0.1299 | 0.0704 | 0.1512 |
Moreover, we know that $P(C=1)=\sum_A P(A,C=1)$
In [11]:
pAC.extract({"c":1}).sum()
Out[11]:
0.35154423108981625
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]:
|
|
|
---|---|---|
0.3694 | 0.2004 | 0.4302 |
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.sumIn("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)$"])
|
|
| |
---|---|---|---|
0.4384 | 0.3823 | 0.1793 | |
0.3694 | 0.2004 | 0.4302 | |
0.2056 | 0.3442 | 0.4502 |
|
|
|
---|---|---|
0.3694 | 0.2004 | 0.4302 |
Likelihood $P(A=2|C)$¶
A likelihood can also be found in this matrix.
In [14]:
pAgivenC.extract({'a':2})
Out[14]:
|
|
|
---|---|---|
0.1793 | 0.4302 | 0.4502 |
A likelihood does not have to sum to 1. It is not relevant to normalize it.
In [15]:
pAgivenC.sumIn(["a"])
Out[15]:
|
|
|
---|---|---|
1.0133 | 0.9269 | 1.0598 |
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()