Probabilistic Graph Model Basics
Janpu Hou
December 7, 2017
- Bayesian networks
- The chest clinic narrative: Step 1: Setup DAG model
- 1. Setup DAG-based models: Specify chest clinic network.
- 2. Conditional probability tables (CPTs)
- 3. Querying the network
- 4. Setting evidence
- 5. Querying the network after evidence
- 6. Set additional evidence and query again
- 7. Get joint dist of tub, lung, bronc given evidence
- 8. Get distribution of tub given lung, bronc and evidence
- 9. Remove evidence
- The chest clinic narrative: Step 2: Brute Force Calculations
- Message passing : a small example
Bayesian networks
Lauritzen and Spiegehalter (1988) present the following narrative:
Shortness-of-breath (dyspnoea ) may be due to tuberculosis, lung cancer or bronchitis, or none of them, or more than one of them.
A recent visit to Asia increases the chances of tuberculosis, while smoking is known to be a risk factor for both lung cancer and bronchitis.
The results of a single chest Xray do not discriminate between lung cancer and tuberculosis, as neither does the presence or absence of dyspnoea.
The narrative can be pictured as a DAG (Directed Acyclic Graph)
The chest clinic narrative: Step 1: Setup DAG model
1. Setup DAG-based models: Specify chest clinic network.
Directed acyclic graphs
# install.packages("gRbase", dependencies = TRUE)
# install.packages("gRain", dependencies = TRUE)
# install.packages("gRim", dependencies = TRUE)
# source("http://bioconductor.org/biocLite.R");
# biocLite(c("graph","RBGL","Rgraphviz"))[,libPaths()]
library(gRbase)## Warning: package 'gRbase' was built under R version 3.4.3library(gRain)## Warning: package 'gRain' was built under R version 3.4.3library(gRim)## Warning: package 'gRim' was built under R version 3.4.3library(igraph)## Warning: package 'igraph' was built under R version 3.4.3##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionyn <- c("yes","no")
a <- cptable(~asia, values=c(1,99),levels=yn)
t.a <- cptable(~tub+asia, values=c(5,95,1,99),levels=yn)
s <- cptable(~smoke, values=c(5,5), levels=yn)
l.s <- cptable(~lung+smoke, values=c(1,9,1,99), levels=yn)
b.s <- cptable(~bronc+smoke, values=c(6,4,3,7), levels=yn)
e.lt <- cptable(~either+lung+tub,values=c(1,0,1,0,1,0,0,1),levels=yn)
x.e <- cptable(~xray+either, values=c(98,2,5,95), levels=yn)
d.be <- cptable(~dysp+bronc+either, values=c(9,1,7,3,8,2,1,9), levels=yn)
cpt.list <- compileCPT(list(a, t.a, s, l.s, b.s, e.lt, x.e, d.be))
cpt.list## CPTspec with probabilities:
## P( asia )
## P( tub | asia )
## P( smoke )
## P( lung | smoke )
## P( bronc | smoke )
## P( either | lung tub )
## P( xray | either )
## P( dysp | bronc either )bnet <- grain(cpt.list)
plot(bnet)
cpt.list$asia## asia
## yes no
## 0.01 0.99
## attr(,"class")
## [1] "parray" "array"cpt.list$tub## asia
## tub yes no
## yes 0.05 0.01
## no 0.95 0.99
## attr(,"class")
## [1] "parray" "array"ftable(cpt.list$either, row.vars=1) # Notice: logical variable## lung yes no
## tub yes no yes no
## either
## yes 1 1 1 0
## no 0 0 0 12. Conditional probability tables (CPTs)
CPTs are just multiway arrays WITH dimnames attribute
# For example p(tja):
yn <- c("yes","no");
x <- c(5,95,1,99)
# Vanilla R
t.a <- array(x, dim=c(2,2), dimnames=list(tub=yn,asia=yn))
t.a## asia
## tub yes no
## yes 5 1
## no 95 99# Alternative (partial) specification
t.a <- cptable(~tub | asia, values=c(5,95,1,99), levels=yn)
t.a## {v,pa(v)} : chr [1:2] "tub" "asia"
## <NA> <NA>
## yes 5 1
## no 95 99# Create network from CPT list:
bnet <- grain(cpt.list)
# Compile network (details follow)
bnet <- compile(bnet)
bnet## Independence network: Compiled: TRUE Propagated: FALSE
## Nodes: chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"3. Querying the network
Query network to find marginal probabilities of diseases
querygrain(bnet, nodes=c("tub","lung","bronc"))## $tub
## tub
## yes no
## 0.0104 0.9896
##
## $lung
## lung
## yes no
## 0.055 0.945
##
## $bronc
## bronc
## yes no
## 0.45 0.554. Setting evidence
Set evidence and query network again
bnet.ev<-setEvidence(bnet, nodes = c("asia","dysp"),states = c("yes","yes"))
bnet.ev## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"
## Evidence:
## nodes is.hard.evidence hard.state
## 1 asia TRUE yes
## 2 dysp TRUE yes
## pEvidence: 0.004501# or
bnet.f <- setFinding(bnet, nodes=c('asia', 'dysp'), state=c('yes','yes'))
bnet.f## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"
## Evidence:
## nodes is.hard.evidence hard.state
## 1 asia TRUE yes
## 2 dysp TRUE yes
## pEvidence: 0.0045015. Querying the network after evidence
querygrain(bnet.ev, nodes=c("tub","lung","bronc"))## $tub
## tub
## yes no
## 0.08775096 0.91224904
##
## $lung
## lung
## yes no
## 0.09952515 0.90047485
##
## $bronc
## bronc
## yes no
## 0.8114021 0.1885979querygrain(bnet.f, nodes=c('lung', 'tub', 'bronc'), type='joint')## , , bronc = yes
##
## lung
## tub yes no
## yes 0.003149038 0.04183722
## no 0.059831718 0.70658410
##
## , , bronc = no
##
## lung
## tub yes no
## yes 0.001827219 0.04093749
## no 0.034717170 0.11111605
##
## attr(,"class")
## [1] "parray" "array"6. Set additional evidence and query again
bnet.ev<-setEvidence(bnet.ev, nodes = "xray", states = "yes")
querygrain(bnet.ev, nodes=c("tub","lung","bronc"))## $tub
## tub
## yes no
## 0.3917117 0.6082883
##
## $lung
## lung
## yes no
## 0.4442705 0.5557295
##
## $bronc
## bronc
## yes no
## 0.6288218 0.3711782# Probability of observing the evidence (the normalizing constant)
pEvidence(bnet.ev)## [1] 0.00098822687. Get joint dist of tub, lung, bronc given evidence
x<-querygrain(bnet.ev, nodes=c("tub","lung","bronc"), type="joint")
ftable(x, row.vars=1)## lung yes no
## bronc yes no yes no
## tub
## yes 0.014056997 0.008156529 0.186757240 0.182740955
## no 0.267082934 0.154974048 0.160924606 0.0253066928. Get distribution of tub given lung, bronc and evidence
x<-querygrain(bnet.ev, nodes=c("tub","lung","bronc"), type="conditional")
ftable(x, row.vars=1)## bronc yes no
## tub yes no yes no
## lung
## yes 0.0500000 0.9500000 0.0500000 0.9500000
## no 0.5371498 0.4628502 0.8783611 0.12163899. Remove evidence
bnet.ev<-retractEvidence(bnet.ev, nodes="asia")
bnet.ev## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"
## Evidence:
## nodes is.hard.evidence hard.state
## 1 dysp TRUE yes
## 2 xray TRUE yes
## pEvidence: 0.070670The chest clinic narrative: Step 2: Brute Force Calculations
In principle, we can find e.g. p(b|a+; d+) by brute force calculations.
Recall: We have a collection of conditional probability tables (CPTs) of the form p(v|pa(v)):{ p(a); p(t|a); p(s); p(l|s); p(b|s); p(e|t; l); p(d|e; b); p(x|e)}
Brute force computations: 1) Form the joint distribution p(V ) by multiplying the CPTs p(V) = p(a)p(t|a)p(s)p(l|s)p(b|s)p(e|t; l)p(d|e; b)p(x|e):
This gives p(V) represented by a table with giving a table with 2^8 = 256 entries.
- Find the marginal distribution p(a; b; d) by marginalizing p(V) = p(a; t; s; k; e; b; x; d) p(a; b; d) = sum(p(t; s; k; e; b; x; d))
This is table with 2^3 = 8 entries.
- Lastly notice that p(b|a+; d+) / p(a+; b; d+). Hence from p(a; b; d) we must extract those entries consistent with a = a+ and d = d+ and normalize the result.
Alternatively (and easier): Set all entries not consistent with a = a+ and d = d+ in p(a; b; d) equal to zero.
1. Form the joint distribution p(V)
## collection of CPTs: p(v|pa(v))
cpt.list## CPTspec with probabilities:
## P( asia )
## P( tub | asia )
## P( smoke )
## P( lung | smoke )
## P( bronc | smoke )
## P( either | lung tub )
## P( xray | either )
## P( dysp | bronc either )## form joint p(V)= prod p(v|pa(v))
joint <- cpt.list$asia
joint## asia
## yes no
## 0.01 0.99
## attr(,"class")
## [1] "parray" "array"for (i in 2:length(cpt.list)){joint <- tableMult( joint, cpt.list[[i]] )}
joint## , , either = yes, xray = yes, lung = yes, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 1.323e-05 6.86e-06
## no 1.470e-06 2.94e-06
##
## , , either = no, xray = yes, lung = yes, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 2.7e-07 1.4e-07
## no 3.0e-08 6.0e-08
##
## , , either = no, xray = no, lung = yes, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0.00011907 6.174e-05
## no 0.00001323 2.646e-05
##
## , , either = no, xray = yes, lung = no, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = no, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 2.43e-06 1.26e-06
## no 2.70e-07 5.40e-07
##
## , , either = no, xray = no, lung = no, tub = yes, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = yes, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0.00025137 0.00013034
## no 0.00002793 0.00005586
##
## , , either = no, xray = yes, lung = yes, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 5.13e-06 2.66e-06
## no 5.70e-07 1.14e-06
##
## , , either = no, xray = no, lung = yes, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = yes, lung = no, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 1.026e-04 8.550e-06
## no 2.565e-05 7.695e-05
##
## , , either = yes, xray = no, lung = no, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = no, lung = no, tub = no, smoke = yes, asia = yes
##
## bronc
## dysp yes no
## yes 0.00194940 0.00016245
## no 0.00048735 0.00146205
##
## , , either = yes, xray = yes, lung = yes, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 6.615e-07 1.2005e-06
## no 7.350e-08 5.1450e-07
##
## , , either = no, xray = yes, lung = yes, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 1.35e-08 2.45e-08
## no 1.50e-09 1.05e-08
##
## , , either = no, xray = no, lung = yes, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 6.54885e-05 0.0001188495
## no 7.27650e-06 0.0000509355
##
## , , either = no, xray = yes, lung = no, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = no, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 1.3365e-06 2.4255e-06
## no 1.4850e-07 1.0395e-06
##
## , , either = no, xray = no, lung = no, tub = yes, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = yes, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 1.25685e-05 2.28095e-05
## no 1.39650e-06 9.77550e-06
##
## , , either = no, xray = yes, lung = yes, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 2.565e-07 4.655e-07
## no 2.850e-08 1.995e-07
##
## , , either = no, xray = no, lung = yes, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = yes, lung = no, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 5.64300e-05 1.645875e-05
## no 1.41075e-05 1.481288e-04
##
## , , either = yes, xray = no, lung = no, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = no, lung = no, tub = no, smoke = no, asia = yes
##
## bronc
## dysp yes no
## yes 0.0010721700 0.0003127163
## no 0.0002680425 0.0028144462
##
## , , either = yes, xray = yes, lung = yes, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.000261954 0.000135828
## no 0.000029106 0.000058212
##
## , , either = no, xray = yes, lung = yes, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 5.346e-06 2.772e-06
## no 5.940e-07 1.188e-06
##
## , , either = no, xray = no, lung = yes, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.002357586 0.001222452
## no 0.000261954 0.000523908
##
## , , either = no, xray = yes, lung = no, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = no, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 4.8114e-05 2.4948e-05
## no 5.3460e-06 1.0692e-05
##
## , , either = no, xray = no, lung = no, tub = yes, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = yes, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.025933446 0.013446972
## no 0.002881494 0.005762988
##
## , , either = no, xray = yes, lung = yes, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.000529254 0.000274428
## no 0.000058806 0.000117612
##
## , , either = no, xray = no, lung = yes, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = yes, lung = no, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.01058508 0.00088209
## no 0.00264627 0.00793881
##
## , , either = yes, xray = no, lung = no, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = no, lung = no, tub = no, smoke = yes, asia = no
##
## bronc
## dysp yes no
## yes 0.20111652 0.01675971
## no 0.05027913 0.15083739
##
## , , either = yes, xray = yes, lung = yes, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 1.30977e-05 2.37699e-05
## no 1.45530e-06 1.01871e-05
##
## , , either = no, xray = yes, lung = yes, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 2.673e-07 4.851e-07
## no 2.970e-08 2.079e-07
##
## , , either = no, xray = no, lung = yes, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0.0012966723 0.002353220
## no 0.0001440747 0.001008523
##
## , , either = no, xray = yes, lung = no, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = no, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 2.64627e-05 4.80249e-05
## no 2.94030e-06 2.05821e-05
##
## , , either = no, xray = no, lung = no, tub = yes, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = yes, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0.0012966723 0.002353220
## no 0.0001440747 0.001008523
##
## , , either = no, xray = yes, lung = yes, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = no, lung = yes, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 2.64627e-05 4.80249e-05
## no 2.94030e-06 2.05821e-05
##
## , , either = no, xray = no, lung = yes, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = yes, xray = yes, lung = no, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = yes, lung = no, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0.005821794 0.001698023
## no 0.001455448 0.015282209
##
## , , either = yes, xray = no, lung = no, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0 0
## no 0 0
##
## , , either = no, xray = no, lung = no, tub = no, smoke = no, asia = no
##
## bronc
## dysp yes no
## yes 0.11061409 0.03226244
## no 0.02765352 0.29036198dim(joint)## dysp bronc either xray lung tub smoke asia
## 2 2 2 2 2 2 2 2head( as.data.frame.table( joint ) )## dysp bronc either xray lung tub smoke asia Freq
## 1 yes yes yes yes yes yes yes yes 1.323e-05
## 2 no yes yes yes yes yes yes yes 1.470e-06
## 3 yes no yes yes yes yes yes yes 6.860e-06
## 4 no no yes yes yes yes yes yes 2.940e-06
## 5 yes yes no yes yes yes yes yes 0.000e+00
## 6 no yes no yes yes yes yes yes 0.000e+002. Form marginal p(a,b,d) by marginalization
marg <- tableMargin(joint, ~asia+bronc+dysp)
dim( marg )## asia bronc dysp
## 2 2 2ftable( marg )## dysp yes no
## asia bronc
## yes yes 0.003652425 0.000847575
## no 0.000848950 0.004651050
## no yes 0.359932815 0.085567185
## no 0.071536410 0.4729635903. Set entries not consistent with evidence equal to zero
Evidence:(asia=yes and dysp=yes)
marg <- tableSetSliceValue(marg, c("asia","dysp"), c("yes","yes"), complement=T)
ftable(marg)## dysp yes no
## asia bronc
## yes yes 0.003652425 0.000000000
## no 0.000848950 0.000000000
## no yes 0.000000000 0.000000000
## no 0.000000000 0.0000000004. Obtain the results
result <- tableMargin(marg, ~bronc);
result <- result / sum( result ); result## bronc
## yes no
## 0.8114021 0.1885979Message passing : a small example
require(gRbase); require(Rgraphviz)## Loading required package: Rgraphviz## Warning: package 'Rgraphviz' was built under R version 3.4.2## Loading required package: graph## Warning: package 'graph' was built under R version 3.4.2## Loading required package: BiocGenerics## Warning: package 'BiocGenerics' was built under R version 3.4.2## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:igraph':
##
## normalize, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colMeans,
## colnames, colSums, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, Map, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
## rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min##
## Attaching package: 'graph'## The following objects are masked from 'package:igraph':
##
## degree, edges, intersection## Loading required package: gridd<-dag( ~smoke + bronc|smoke + dysp|bronc ); plot(d)
Direct Method: Calculate the Joint Distribution
Recall that the joint distribution is p(s; b; d) = p(s)p(b|s)p(d|b)
yn <- c("yes","no")
s <- parray("smoke", list(yn), c(.5, .5))
b.s <- parray(c("bronc","smoke"), list(yn,yn), c(6,4, 3,7))
d.b <- parray(c("dysp","bronc"), list(yn, yn), c(9,1, 2,8))
s; b.s; d.b## smoke
## yes no
## 0.5 0.5
## attr(,"class")
## [1] "parray" "array"## smoke
## bronc yes no
## yes 6 3
## no 4 7
## attr(,"class")
## [1] "parray" "array"## bronc
## dysp yes no
## yes 9 2
## no 1 8
## attr(,"class")
## [1] "parray" "array"joint <- tableMult( tableMult(s, b.s), d.b) ; ftable(joint)## smoke yes no
## dysp bronc
## yes yes 27.0 13.5
## no 4.0 7.0
## no yes 3.0 1.5
## no 16.0 28.0Indirect Method: Message Passing:
We no longer need the DAG. Instead we use an undirected graph to dictate the message passing: The moral graph is obtained by 1) marrying parents and 2) dropping directions.
The moral graph is (in this case) triangulated which means that the cliques can be organized in a tree called a junction tree.
Defne:
q1(s; b) = p(s)p(b|s) and q2(b; d) = p(d|b)
and we have
p(s; b; d) = p(s)p(b|s)p(d|b) = q1(s; b)q2(b; d)
We see that the q functions are defned on the cliques of the moral graph or - equivalently - on the nodes of the junction tree. The q functions are called potentials; they are non-negative functions but they are typically not probabilities and they are hence difficult to interpret.We can think of the q functions as interactions.
dm <-moralize(d);
jtree<-ug(~smoke.bronc:bronc.dysp);
par(mfrow=c(1,3)); plot(d); plot(dm); plot(jtree)
q1.sb <- tableMult(s, b.s); q1.sb## smoke
## bronc yes no
## yes 3 1.5
## no 2 3.5q2.bd <- d.b; q2.bd## bronc
## dysp yes no
## yes 9 2
## no 1 8
## attr(,"class")
## [1] "parray" "array"The factorization p(s; b; d) = q1(s; b)q2(b; d) is called a clique potential representation. Goal: We shall operate on qfunctions such that at the end they will contain the marginal distributions, i.e. q1(s; b) = p(s; b); q2(b; d) = p(b; d)
Collect Evidence
work inwards towards the root
# work inwards towards the root
q1.b <- tableMargin(q1.sb, "bronc"); q1.b## bronc
## yes no
## 4.5 5.5Therefore, if we update potentials as q1(s; b) q1(s; b)=q1(b); q2(b; d) q2(b; d)q1(b) and we obtain new potentials defined on the cliques of the junction tree.
We still have p(s; b; d) = q1(s; b)q2(b; d)
Updating of potentials
q1(s; b) <— q1(s; b)/q1(b);
q2(b; d) <— q2(b; d)q1(b)
# Updating of potentials
q2.bd <- tableMult(q2.bd, q1.b); q2.bd## dysp
## bronc yes no
## yes 40.5 4.5
## no 11.0 44.0q1.sb <- tableDiv(q1.sb, q1.b); q1.sb## smoke
## bronc yes no
## yes 0.6666667 0.3333333
## no 0.3636364 0.6363636Distribute Evidence
Next work outwards from the root.
p(s; b; d) = q1(s; b)q2(b; d)
We set q1(s; b)<—q1(s; b)q2(b) and have p(s; b; d) = q1(s; b)q2(b; d)
q2.b <- tableMargin(q2.bd, "bronc"); q2.b## bronc
## yes no
## 45 55q1.sb <- tableMult(q1.sb, q2.b); q1.sb## smoke
## bronc yes no
## yes 30 15
## no 20 35The form p(s; b; d) = q1(s; b)q2(b; d) is called the clique marginal representation and the main point is now that q1(s; b) = p(s; b); q2(b; d) = p(b; d) and q1 and q2 fit on their marginals, i.e. q1(b) = q2(b).
joint## , , smoke = yes
##
## bronc
## dysp yes no
## yes 27 4
## no 3 16
##
## , , smoke = no
##
## bronc
## dysp yes no
## yes 13.5 7
## no 1.5 28Claim: After these steps q1(s; b) = p(s; b) and q2(b; d) = p(b; d).
q1.sb## smoke
## bronc yes no
## yes 30 15
## no 20 35tableMargin(joint, c("smoke","bronc"))## bronc
## smoke yes no
## yes 30 20
## no 15 35q2.bd## dysp
## bronc yes no
## yes 40.5 4.5
## no 11.0 44.0tableMargin(joint, c("bronc","dysp"))## dysp
## bronc yes no
## yes 40.5 4.5
## no 11.0 44.0## Now we can obtain, e.g. p(b) as
tableMargin(q1.sb, "bronc")## bronc
## yes no
## 45 55tableMargin(q2.bd, "bronc")## bronc
## yes no
## 45 55