Scores, Chi2, etc. with BNLearner

Creative Commons License

aGrUM

interactive online version

In [1]:
import os

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

Generating the database for scoring

In [2]:
bn=gum.loadBN("res/asia.bif")
bn
Out[2]:
G tuberculosis tuberculosis tuberculos_or_cancer tuberculos_or_cancer tuberculosis->tuberculos_or_cancer bronchitis bronchitis dyspnoea dyspnoea bronchitis->dyspnoea tuberculos_or_cancer->dyspnoea positive_XraY positive_XraY tuberculos_or_cancer->positive_XraY lung_cancer lung_cancer lung_cancer->tuberculos_or_cancer visit_to_Asia visit_to_Asia visit_to_Asia->tuberculosis smoking smoking smoking->bronchitis smoking->lung_cancer
In [3]:
# we create a quite large database
gum.generateSample(bn,500000,"out/sample_score.csv",False)
Out[3]:
-1613207.0068915116

Testing d-separations using chi2 in the database

In [4]:
# do not forget that the generation process above is random : from time to time, the tests my not be correct...
def isIndep(pvalue):
    return pvalue>=0.05

def testIndepFromChi2(learner,var1,var2,kno=[]):
    """
    Just prints the resultat of the chi2
    """
    stat,pvalue=learner.chi2(var1,var2,kno)
    if len(kno)==0:
        print("From Chi2 tests, is '{}' indep from '{}' ==> {}".format(var1,var2,isIndep(pvalue)))
    else:
        print("From Chi2 tests, is '{}' indep from '{}' given {} : {}".format(var1,var2,kno,isIndep(pvalue)))

learner=gum.BNLearner("out/sample_score.csv")

testIndepFromChi2(learner,"visit_to_Asia","smoking")
testIndepFromChi2(learner,"visit_to_Asia","smoking",['tuberculos_or_cancer'])
testIndepFromChi2(learner,"visit_to_Asia","smoking",['positive_XraY'])
testIndepFromChi2(learner,"dyspnoea","smoking")
testIndepFromChi2(learner,"dyspnoea","smoking",["lung_cancer","bronchitis"])

From Chi2 tests, is 'visit_to_Asia' indep from 'smoking' ==> True
From Chi2 tests, is 'visit_to_Asia' indep from 'smoking' given ['tuberculos_or_cancer'] : False
From Chi2 tests, is 'visit_to_Asia' indep from 'smoking' given ['positive_XraY'] : False
From Chi2 tests, is 'dyspnoea' indep from 'smoking' ==> False
From Chi2 tests, is 'dyspnoea' indep from 'smoking' given ['lung_cancer', 'bronchitis'] : True

Evolution of chi2 p-values w.r.t the size of the database (in Asia)

In [5]:
def consolidationIndepFromChi2(bn,size,lindep,nbr=20):
    """
    Using $nbr$ generated databases of size $size$ from the bn $bn$,
    consolidate the p-value for a list $lindep$ of conditional independence to test.

    return the list of consolidated pValues
    """
    pvalue_cumul=[0.0]*len(lindep)
    for i in range(nbr):
        gum.generateSample(bn,size,"out/sample_score.csv",False)
        learner=gum.BNLearner("out/sample_score.csv")
        for i,(var1,var2,kno) in enumerate(lindep):
            stat,pvalue=learner.chi2(var1,var2,kno)
            pvalue_cumul[i]+=pvalue
    return [p/nbr for p in pvalue_cumul]

sizes=[50,100,500,1000,2000,5000,10000,20000,50000,100000,200000]
pvalues1,pvalues2,pvalues3,pvalues4,pvalues5,pvalues6=zip(*[consolidationIndepFromChi2(bn,siz,
                                     [("visit_to_Asia","smoking",['tuberculos_or_cancer']),
                                      ("visit_to_Asia","smoking",[]),
                                      ("dyspnoea","smoking",[]),
                                      ("dyspnoea","smoking",["lung_cancer","bronchitis"]),
                                      ("tuberculosis","bronchitis",[]),
                                      ("tuberculosis","bronchitis",["dyspnoea"])])
                                          for siz in sizes])

In [6]:
%matplotlib inline
from pylab import *
import matplotlib.pyplot as plt
import matplotlib.patches as patches

fig=figure(figsize=(10,6))
ax  = fig.add_subplot(1, 1, 1)

ax.plot(sizes,pvalues1,label="NOT(A indep S given TorC)", linestyle='dashed')
ax.plot(sizes,pvalues2,label="A indep S")
ax.plot(sizes,pvalues3,label="NOT(D indep S)", linestyle='dashed')
ax.plot(sizes,pvalues4,label="D indep S given L,B")
ax.plot(sizes,pvalues5,label="T indep B")
ax.plot(sizes,pvalues6,label="NOT(T indep B given D)", linestyle='dashed')

ax.tick_params(rotation=90)
ax.set_xlabel("data size")
ax.set_ylabel("pValue")
ax.legend(bbox_to_anchor=(0.15, 0.88, 0.7, .102), loc=3,ncol=3, mode="expand", borderaxespad=0.)

rect = patches.Rectangle((0,0),max(sizes),0.05,linewidth=1,edgecolor='#FF8888',facecolor='#FF8888')
ax.add_patch(rect)
ax.annotate("Critical region",xytext=(190000,0.2),xy=(190000,0.05),
            ha="right", va="center",
            arrowprops=dict(arrowstyle="->",
                            connectionstyle="arc3,rad=-0.15"
                           ),
            bbox=dict(boxstyle="square", fc="w"))

ax.set_title("Chi2 pvalue=f(datasize)")
gnb.flow.add(fig)
gnb.flow.add(gnb.getBN(gum.fastBN("A->T->TorC->X;S->C->TorC->D<-B<-S")))
gnb.flow.new_line()

fig
Out[6]:
../_images/notebooks_36-Learning_Chi2AndScoresFromBNLearner_10_0.svg

Testing d-separations using G2 in the database

In [7]:
def testIndepFromG2(learner,var1,var2,kno=[]):
    """
    Just prints the resultat of the G2
    """
    stat,pvalue=learner.G2(var1,var2,kno)
    if len(kno)==0:
        print("From G2 tests, is '{}' indep from '{}' ==> {}".format(var1,var2,isIndep(pvalue)))
    else:
        print("From G2 tests, is '{}' indep from '{}' given {} : {}".format(var1,var2,kno,isIndep(pvalue)))

learner=gum.BNLearner("out/sample_score.csv")

testIndepFromG2(learner,"visit_to_Asia","smoking")
testIndepFromG2(learner,"visit_to_Asia","smoking",['tuberculos_or_cancer'])
testIndepFromG2(learner,"visit_to_Asia","smoking",['positive_XraY'])
testIndepFromG2(learner,"dyspnoea","smoking")
testIndepFromG2(learner,"dyspnoea","smoking",["lung_cancer","bronchitis"])

From G2 tests, is 'visit_to_Asia' indep from 'smoking' ==> True
From G2 tests, is 'visit_to_Asia' indep from 'smoking' given ['tuberculos_or_cancer'] : False
From G2 tests, is 'visit_to_Asia' indep from 'smoking' given ['positive_XraY'] : False
From G2 tests, is 'dyspnoea' indep from 'smoking' ==> False
From G2 tests, is 'dyspnoea' indep from 'smoking' given ['lung_cancer', 'bronchitis'] : True

Evolution of G2 p-values w.r.t the size of the database (in Asia)

In [8]:
def consolidationIndepFromG2(bn,size,lindep,nbr=20):
    """
    Using $nbr$ generated databases of size $size$ from the bn $bn$,
    consolidate the p-value for a list $lindep$ of conditional independence to test.

    return the list of consolidated pValues
    """
    pvalue_cumul=[0.0]*len(lindep)
    for i in range(nbr):
        gum.generateSample(bn,size,"out/sample_chi2.csv",False)
        learner=gum.BNLearner("out/sample_chi2.csv")
        for i,(var1,var2,kno) in enumerate(lindep):
            stat,pvalue=learner.G2(var1,var2,kno)
            pvalue_cumul[i]+=pvalue
    return [p/nbr for p in pvalue_cumul]

sizes=[50,100,500,1000,2000,5000,10000,20000,50000,100000,200000]
pvalues1,pvalues2,pvalues3,pvalues4,pvalues5,pvalues6=zip(*[consolidationIndepFromG2(bn,siz,
                                                                                     [("visit_to_Asia","smoking",['tuberculos_or_cancer']),
                                                                                      ("visit_to_Asia","smoking",[]),
                                                                                      ("dyspnoea","smoking",[]),
                                                                                      ("dyspnoea","smoking",["lung_cancer","bronchitis"]),
                                                                                      ("tuberculosis","bronchitis",[]),
                                                                                      ("tuberculosis","bronchitis",["dyspnoea"])])
                                                            for siz in sizes])
In [9]:
fig=figure(figsize=(10,6))
ax  = fig.add_subplot(1, 1, 1)

ax.plot(sizes,pvalues1,label="NOT(A indep S given TorC)", linestyle='dashed')
ax.plot(sizes,pvalues2,label="A indep S")
ax.plot(sizes,pvalues3,label="NOT(D indep S)", linestyle='dashed')
ax.plot(sizes,pvalues4,label="D indep S given L,B")
ax.plot(sizes,pvalues5,label="T indep B")
ax.plot(sizes,pvalues6,label="NOT(T indep B given D)", linestyle='dashed')

ax.tick_params(rotation=90)
ax.set_xlabel("data size")
ax.set_ylabel("pValue")
ax.legend(bbox_to_anchor=(0.15, 0.83, 0.7, .102), loc=3,ncol=2, mode="expand", borderaxespad=0.)

rect = patches.Rectangle((0,0),max(sizes),0.05,linewidth=1,edgecolor='#FF8888',facecolor='#FF8888')
ax.add_patch(rect)
ax.annotate("Critical region",xytext=(190000,0.2),xy=(190000,0.05),
            ha="right", va="center",
            arrowprops=dict(arrowstyle="->",
                            connectionstyle="arc3,rad=-0.15"
                           ),
            bbox=dict(boxstyle="square", fc="w"))

ax.set_title("G2 pvalue=f(datasize)")
gnb.flow.add(fig)
gnb.flow.add(gnb.getBN(gum.fastBN("A->T->TorC->X;S->C->TorC->D<-B<-S")))
gnb.flow.new_line()
Out[9]:
<pyAgrum.lib.notebook.FlowLayout at 0x7f479413e710>

Conditional joint log-likelihood

With BNLearner, you can also check the joint (condtional) log-likelihood in the base

In [10]:
bn
Out[10]:
G tuberculosis tuberculosis tuberculos_or_cancer tuberculos_or_cancer tuberculosis->tuberculos_or_cancer bronchitis bronchitis dyspnoea dyspnoea bronchitis->dyspnoea tuberculos_or_cancer->dyspnoea positive_XraY positive_XraY tuberculos_or_cancer->positive_XraY lung_cancer lung_cancer lung_cancer->tuberculos_or_cancer visit_to_Asia visit_to_Asia visit_to_Asia->tuberculosis smoking smoking smoking->bronchitis smoking->lung_cancer
In [11]:
siz=10000
gum.generateSample(bn,siz,"out/sample_score.csv",False)
learner=gum.BNLearner("out/sample_score.csv")

def affLL(learner,s1,s2=[]):
    if len(s2)==0:
        print("{} : {}".format(s1,learner.logLikelihood(s1)))
    else:
        print("{}|{} : {}".format(s1,s2,learner.logLikelihood(s1,s2)))

def dsepByLL(learner,x,y,z): # is X indep of Y given Z ?
    lxy_z=learner.logLikelihood([x,y],[z])
    lx_z=learner.logLikelihood([x],[z])
    ly_z=learner.logLikelihood([y],[z])
    print("{} indep {} given {} : {}".format(x,y,z,lxy_z-lx_z-ly_z))

print("Condional Joint LogLikelihood")
affLL(learner,["lung_cancer","bronchitis","smoking"])
affLL(learner,["smoking"])
affLL(learner,["lung_cancer","bronchitis"],["smoking"])

print("--------------")
print("LL indep test")
dsepByLL(learner,"lung_cancer","bronchitis","smoking")
dsepByLL(learner,"tuberculos_or_cancer","bronchitis","dyspnoea")
Condional Joint LogLikelihood
['lung_cancer', 'bronchitis', 'smoking'] : -21960.28716457141
['smoking'] : -9999.95123685268
['lung_cancer', 'bronchitis']|['smoking'] : -11960.335927718732
--------------
LL indep test
lung_cancer indep bronchitis given smoking : 1.9743569250294968
tuberculos_or_cancer indep bronchitis given dyspnoea : 188.35824992442758

Evolution of conditional log-likelihood w.r.t the size of the database (in Asia)

In [12]:
def consolidationIndepFromLL(bn,size,lindep,nbr=20):
    """
    Using $nbr$ generated databases of size $size$ from the bn $bn$,
    consolidate the logLikelihoos for a list $lindep$ of conditional independence to test.

    return the list of consolidated pValues
    """
    LL_cumul=[0.0]*len(lindep)
    for i in range(nbr):
        gum.generateSample(bn,size,"out/sample_score.csv",False)
        learner=gum.BNLearner("out/sample_score.csv")
        for i,(var1,var2,kno) in enumerate(lindep):
            LL12=learner.logLikelihood([var1,var2],kno)
            LL1=learner.logLikelihood([var1],kno)
            LL2=learner.logLikelihood([var2],kno)
            LL_cumul[i]+=(LL12-LL1-LL2)/size
    return [p/nbr for p in LL_cumul]

sizes=[50,100,500,1000,2000,5000,10000,20000,50000,100000,200000]
LL1,LL2,LL3,LL4,LL5,LL6=zip(*[consolidationIndepFromLL(bn,siz,
                                     [("visit_to_Asia","smoking",['tuberculos_or_cancer']),
                                      ("visit_to_Asia","smoking",[]),
                                      ("dyspnoea","smoking",[]),
                                      ("dyspnoea","smoking",["lung_cancer","bronchitis"]),
                                      ("tuberculosis","bronchitis",[]),
                                      ("tuberculosis","bronchitis",["dyspnoea"])])
                                          for siz in sizes])



In [13]:
%matplotlib inline
from pylab import *
import matplotlib.pyplot as plt
import matplotlib.patches as patches

fig=figure(figsize=(10,6))
ax  = fig.add_subplot(1, 1, 1)

ax.plot(sizes,LL1,label="NOT(A indep S given TorC)", linestyle='dashed')
ax.plot(sizes,LL2,label="A indep S")
ax.plot(sizes,LL3,label="NOT(D indep S)", linestyle='dashed')
ax.plot(sizes,LL4,label="D indep S given C,B")
ax.plot(sizes,LL5,label="T indep B")
ax.plot(sizes,LL6,label="NOT(T indep B given D)", linestyle='dashed')
ax.tick_params(rotation=90)
ax.set_xlabel("data size")
ax.set_ylabel("LL/size")
ax.semilogy()
ax.legend(bbox_to_anchor=(0.15, 0.8, 0.8, .102), loc=3,ncol=3, mode="expand", borderaxespad=0.)

ax.set_title("logLikelihood=f(datasize)")
gnb.flow.add(fig)
gnb.flow.add(gnb.getBN(gum.fastBN("A->T->TorC->X;S->C->TorC->D<-B<-S")))
gnb.flow.new_line();

Comparing the scores

In [14]:
gnb.flow.display()
G A A T T A->T D D TorC TorC TorC->D X X TorC->X T->TorC C C C->TorC B B B->D S S S->C S->B

G A A T T A->T D D TorC TorC TorC->D X X TorC->X T->TorC C C C->TorC B B B->D S S S->C S->B

G A A T T A->T D D TorC TorC TorC->D X X TorC->X T->TorC C C C->TorC B B B->D S S S->C S->B

In [ ]: