3 More Complex Cases:
Hybrid Bayesian Networks
The chapter where we learn to use JAGS in an unusual way
Light theme
Text in these colors
Code in this
font
Fewer memes More screenshots…
So far…. Single family BNs
Ch 1. Multinomial
Ch 2. Multivariate normal
This chapter…. BNs with a mix of probability distributions
What:
Rods (categorical and normal)
Rods (a mix of several distributions)
Crop loss due to pests (a mix of several distributions)
How:
Network defined as a model in JAGS. Conduct queries by providing the conditions of interest
as data.
Using JAGS
library(rjags)
model <- textConnection(“model { … }”)
distribution ~ dnorm(mean, variance)
defined_value <- a + b*x
jags.data <- list(var = values, …)
Fit with jags.model(model, data)
Update to remove burn-in update(model1, n.iter = 10000)
Sample values from the posterior coda.samples(model = model1, variable.names =
"csup",n.iter = 20000, thin = 20)
Creatively extract values of interest [[]], as.matrix()
One line per node
No overwriting/recycling (compiled)
Common practice to index values from
lists or matrices
All variables must appear on the left in
the model or be provided as data
Using JAGS
library(rjags)
model <- textConnection(“model { … }”)
distribution ~dnorm(mean, variance)
defined_value <- a + b*x
jags.data <- list(var = values, …)
Fit with model1 <- jags.model(model, data)
Update to remove burn-in update(model1, n.iter = 10000)
Sample values from the posterior sim <- coda.samples(model = model1,
variable.names = "myvar", n.iter =
20000, thin = 20)
Creatively extract values of interest sim[[]], as.matrix(sim)
Using JAGS to query BNs
library(rjags)
model <- textConnection(“model { … }”)
distribution ~dnorm(mean, variance)
defined_value <- a + b*x
jags.data <- list(var = values, …)
Fit with model1 <- jags.model(model, data)
Update to remove burn-in update(model1, n.iter = 10000)
Sample values from the posterior sim <- coda.samples(model = model1,
variable.names = "myvar", n.iter =
20000, thin = 20)
Creatively extract values of interest sim[[]], as.matrix(sim)
Provide the conditions in
the data list
Request the response
variable from
coda.samples
Key points:
There is no theoretical constraint on having multiple families of distributions in
the same BN.
But, BN-specific packages (bnlearn, ) cannot solve them analytically.
So we rely on MCMC sampling via rjags to approximate the conditional
probabilities.
Not covered: Fitting the BN from data; learning the BN structure from data.
Let’s look at some examples!
Rods
Two suppliers (S) with different mean rod diameters (D)
Rods
Two suppliers (S) with different mean rod diameters (D)
What is the probability that our 6.2 mm
rod was supplied by supplier 1?
Rods
Two suppliers (S) with different mean rod diameters (D)
model:
[S][D|S]
One line per node
Fancy indexing
csup, cdiam are defined, but mu, sp, and sigma are missing
Rods
Two suppliers (S) with different mean rod diameters (D)
Provide values for variables
Fix cdiam to query the BN
Rods
Two suppliers (S) with different mean rod diameters (D)
Provide values for variables
Fix cdiam to query the BN
Exclude burn-in, sample
Extract sample, count results
Rods
Two suppliers (S) with different mean rod diameters (D)
[1] 0.172
What is the probability that our 6.2 mm
rod was supplied by supplier 1?
Discrete rods
We handle the same scenario with a multinomial (categorical) BN if we make rod
diameter discrete
“thin” (S1) “average” (hard to tell) “thick” (S2)
(6.16, 6. 19)
Discrete rods
We handle the same scenario with a multinomial (categorical) BN if we make rod
diameter discrete
“thin” (S1) “average” (hard to tell) “thick” (S2)
(6.16, 6. 19)
Discrete rods
We handle the same scenario with a multinomial (categorical) BN if we make rod
diameter discrete
Conditional probability of S|D
Conditional probability of D|S
Joint,normalize
Discrete rods: Why?
We handle the same scenario with a multinomial (categorical) BN if we make rod
diameter discrete
“thin” (S1) “average” (hard to tell) “thick(S2)
(6.16, 6. 19)
Loss of information where the discrete
categories don’t completely approximate the
curve
But, benefit of fitting analytically ie. no loss of
accuracy due to sampling (and faster
computation)
Do:
Fit good categories
Balance the probable loss of accuracy due to
sampling vs. binning data
Mixed distributions
Pick 0 or 1
Just 4? of 27?
Mixed distributions
Values to provide to JAGS as data
list Produces a similar result to our
earlier setup
Pests -> crop loss
G1 pest population awakens in the spring, make G2 pests
G2 eat the crop and cause crop LOss
G1 population determined by PRevious crop, and winter CLimate
TReatment decisions are made based on G1 population size
(but here we will use different distributions)
Pest example
Pest example
Pest example
Pest example
What can we do with this model?
Sensitivity analysis: modify the
parameter values (dat0) to test their
effect on the model
Ask questions that are answered with a
distribution:
What is the effect of the last crop (PR)on
loss caused by the pest (LO)?
What is the effect of treatment (TR)on loss
(LO)?
What is the effect of the last crop (PR) on loss (LO)?
Query by fixing value of PR in the data list
3 possible values of PR (1:3), so use a loop to fix value as data
This could be done with from the distributions directly, or in JAGS
What is the effect of the last crop (PR) on loss (LO)?
Query by fixing value of PR in the data list
3 possible values of PR (1:3), so use a loop to fix value as data
This could be done with from the distributions directly, or in JAGS
What is the effect of the last crop (PR) on loss (LO)?
We know PR increases LO (we made the model)
More insight from other nodes: If CL is constant,
G1 and G2 increase in proportion to PR and TR
What is the effect of treatment (TR) on loss (LO)?
Loss is higher with treatment. But treatment doesn’t increase loss. It is conditioned
on all other nodes in the BN (TR = 1 occurs when PR, G1 and therefore G2 are
higher).
Why JAGS/BUGS?
vs. using distributions directly
More flexibility to query the BN - easier to condition on multiple possible node
combinations
It is more computationally efficient to use class-limited BNs, when possible
Why not JAGS/BUGS?
average
thin
thick
One step per category, representing
its value range
From the original model,
sp <-c(0.5, 0.5) # Supplier prob.
The discretized probabilities are captured in,
limits <-c(6.16, 6.19)
dsd <-matrix(c(diff(c(0, pnorm(limits, mu[1], sigma), 1)),
diff(c(0, pnorm(limits, mu[2], sigma), 1))),3, 2)
dimnames(dsd) <-list(D = c("thin", "average", "thick"),
S = c("s1", "s2"))
dsd # This is the conditional probability table
dsd[, 1] gives the probabilities of each rod size for supplier 1
thin average thick
0.88493033 0.07913935 0.03593032
New model:
model {
csup ~ dcat(sp);
cdiam ~ dcat(dsd[, csup]);
}
Text code in comments
Text code in comments
Sampling/numerical approximation vs
analytical solution
Fix w/many more samples?
1. Changing the marginals will of course change
the probability value (.182…)
2. Changing the marginals will also change how
well sampling approximates the analytical
result. If the new marginals pushes 6.2 closer to
a distribution mean, sampling better
approximates the true value. If the new
marginals push 6.2 closer to the tails, sampling
will produce a worse approximation of the true
values (or more samples will be required).
From the book with (6.16, 6.19)
With the new limits
The new categories are less clearly
delineated:
new: mostly S1, also mostly S1, mostly S2
vs.
old: mostly S1, coin toss, mostly S2