The ASIA network

The graph below corresponds to the ASIA network, a simple Bayesian model used extensively in educational settings. It was introduced by Lauritzen in 1988 ^{[lauritzen1988local]}.

┌───┐           ┌───┐
│ A │         ┌─┤ S ├─┐
└─┬─┘         │ └───┘ │
  │           │       │
  ▼           ▼       ▼
┌───┐       ┌───┐   ┌───┐
│ T │       │ L │   │ B │
└─┬─┘       └─┬─┘   └─┬─┘
  │   ┌───┐   │       │
  └──►│ E │◄──┘       │
      └─┬─┘           │
┌───┐   │   ┌───┐     │
│ X │◄──┴──►│ D │◄────┘
└───┘       └───┘

The table below explains the meanings of each random variable used in the ASIA network model.

We now demonstrate how to use the TensorInference.jl package for conducting a variety of inference tasks on the Asia network.

Import the TensorInference package, which provides the functionality needed for working with tensor networks and probabilistic graphical models.

using TensorInference

Load the ASIA network model from the asia.uai file located in the examples directory. See Model file format (.uai) for a description of the format of this file.

model = read_model_file(pkgdir(TensorInference, "examples", "asia-network", "model.uai"))

UAIModel(nvars = 8, nfactors = 8)
 cards : [2, 2, 2, 2, 2, 2, 2, 2]
 factors : 
  Factor(1), size = (2,)
  Factor(1, 2), size = (2, 2)
  Factor(3), size = (2,)
  Factor(3, 4), size = (2, 2)
  Factor(3, 5), size = (2, 2)
  Factor(2, 4, 6), size = (2, 2, 2)
  Factor(6, 7), size = (2, 2)
  Factor(5, 6, 8), size = (2, 2, 2)

Create a tensor network representation of the loaded model.

tn = TensorNetworkModel(model)

TensorNetworkModel{Int64, OMEinsum.DynamicNestedEinsum{Int64}, Array{Float64}}
variables: 1, 2, 3, 4, 5, 6, 7, 8
contraction time = 2^6.087, space = 2^2.0, read-write = 2^7.14

Calculate the partition function. Since the factors in this model are normalized, the partition function is the same as the total probability, $1$.

probability(tn) |> first

1.0000000000000002

Calculate the marginal probabilities of each random variable in the model.

marginals(tn)

Dict{Vector{Int64}, Vector{Float64}} with 8 entries:
  [8] => [0.435971, 0.564029]
  [3] => [0.5, 0.5]
  [1] => [0.01, 0.99]
  [5] => [0.45, 0.55]
  [4] => [0.055, 0.945]
  [6] => [0.064828, 0.935172]
  [7] => [0.11029, 0.88971]
  [2] => [0.0104, 0.9896]

Retrieve all the variables in the model.

get_vars(tn)

8-element Vector{Int64}:
 1
 2
 3
 4
 5
 6
 7
 8

Set the evidence: Assume that the "X-ray" result (variable 7) is negative. Since setting the evidence may affect the contraction order of the tensor network, recompute it.

tn = TensorNetworkModel(model, evidence = Dict(7 => 0))

TensorNetworkModel{Int64, OMEinsum.DynamicNestedEinsum{Int64}, Array{Float64}}
variables: 1, 2, 3, 4, 5, 6, 7 (evidence → 0), 8
contraction time = 2^6.0, space = 2^2.0, read-write = 2^7.066

Calculate the maximum log-probability among all configurations.

maximum_logp(tn)

0-dimensional Array{Float64, 0}:
-3.65222179200233

Generate 10 samples from the posterior distribution.

sample(tn, 10)

10-element TensorInference.Samples{Int64}:
 [1, 1, 1, 1, 1, 1, 0, 1]
 [1, 1, 0, 1, 0, 1, 0, 0]
 [1, 1, 0, 0, 0, 0, 0, 0]
 [1, 1, 0, 0, 1, 0, 0, 1]
 [1, 1, 0, 0, 1, 0, 0, 0]
 [1, 1, 1, 1, 1, 1, 0, 1]
 [1, 1, 1, 1, 1, 1, 0, 1]
 [1, 1, 0, 0, 1, 0, 0, 0]
 [1, 1, 0, 1, 0, 1, 0, 0]
 [1, 1, 0, 0, 0, 0, 0, 0]

Retrieve both the maximum log-probability and the most probable configuration.

logp, cfg = most_probable_config(tn)

(-3.65222179200233, [1, 1, 0, 0, 0, 0, 0, 0])

Compute the most probable values of certain variables (e.g., 4 and 7) while marginalizing over others. This is known as Maximum a Posteriori (MAP) estimation.

mmap = MMAPModel(model, evidence=Dict(7=>0), queryvars=[4,7])

MMAPModel{Int64, Array{Float64}}
variables: 4, 7 (evidence → 0)
query variables: 4
contraction time = 2^6.022, space = 2^2.0, read-write = 2^7.033

Get the most probable configurations for variables 4 and 7.

most_probable_config(mmap)

(-2.8754627318176693, [1, 0])

Compute the total log-probability of having lung cancer. The results suggest that the probability is roughly half.

log_probability(mmap, [1, 0]), log_probability(mmap, [0, 0])

(-2.8754627318176693, -2.9206248010671856)

This page was generated using Literate.jl.