Probabilistic Graphical Models

Currently only Bayesian Networks

Eventually others including Pearl's causal models.

Bayesian Networks

See the Wikipedia page on Bayesian networks

Example: Alarm

Define random variables

import axle.probability._

val bools = Vector(true, false)

val B = Variable[Boolean]("Burglary")
val E = Variable[Boolean]("Earthquake")
val A = Variable[Boolean]("Alarm")
val J = Variable[Boolean]("John Calls")
val M = Variable[Boolean]("Mary Calls")

Define Factor for each variable

import spire.math._
import cats.implicits._

val bFactor =
  Factor(Vector(B -> bools), Map(
    Vector(B is true) -> Rational(1, 1000),
    Vector(B is false) -> Rational(999, 1000)))

val eFactor =
  Factor(Vector(E -> bools), Map(
    Vector(E is true) -> Rational(1, 500),
    Vector(E is false) -> Rational(499, 500)))

val aFactor =
  Factor(Vector(B -> bools, E -> bools, A -> bools), Map(
    Vector(B is false, E is false, A is true) -> Rational(1, 1000),
    Vector(B is false, E is false, A is false) -> Rational(999, 1000),
    Vector(B is true, E is false, A is true) -> Rational(940, 1000),
    Vector(B is true, E is false, A is false) -> Rational(60, 1000),
    Vector(B is false, E is true, A is true) -> Rational(290, 1000),
    Vector(B is false, E is true, A is false) -> Rational(710, 1000),
    Vector(B is true, E is true, A is true) -> Rational(950, 1000),
    Vector(B is true, E is true, A is false) -> Rational(50, 1000)))

val jFactor =
  Factor(Vector(A -> bools, J -> bools), Map(
    Vector(A is true, J is true) -> Rational(9, 10),
    Vector(A is true, J is false) -> Rational(1, 10),
    Vector(A is false, J is true) -> Rational(5, 100),
    Vector(A is false, J is false) -> Rational(95, 100)))

val mFactor =
  Factor(Vector(A -> bools, M -> bools), Map(
    Vector(A is true, M is true) -> Rational(7, 10),
    Vector(A is true, M is false) -> Rational(3, 10),
    Vector(A is false, M is true) -> Rational(1, 100),
    Vector(A is false, M is false) -> Rational(99, 100)))

Arrange into a graph

import axle.pgm._
import axle.jung._
import edu.uci.ics.jung.graph.DirectedSparseGraph

// edges: ba, ea, aj, am

val bn: BayesianNetwork[Boolean, Rational, DirectedSparseGraph[BayesianNetworkNode[Boolean, Rational], Edge]] =
  BayesianNetwork.withGraphK2[Boolean, Rational, DirectedSparseGraph](
    Map(
      B -> bFactor,
      E -> eFactor,
      A -> aFactor,
      J -> jFactor,
      M -> mFactor))

Create an SVG visualization

import axle.visualize._

val bnVis  = BayesianNetworkVisualization(bn, 1000, 1000, 20)

Render as SVG file

import axle.web._
import cats.effect._

bnVis.svg[IO]("docwork/images/alarm_bayes.svg").unsafeRunSync()

alarm bayes network

The network can be used to compute the joint probability table:

import axle.math.showRational

val jpt = bn.jointProbabilityTable
jpt.show
// res2: String = """John Calls Alarm Earthquake Burglary Mary Calls
// true       true  true       true     true       1197/1000000000
// true       true  true       true     false      513/1000000000
// true       true  true       false    true       1825173/5000000000
// true       true  true       false    false      782217/5000000000
// true       true  false      true     true       1477539/2500000000
// true       true  false      true     false      633231/2500000000
// true       true  false      false    true       31405563/50000000000
// true       true  false      false    false      13459527/50000000000
// true       false true       true     true       1/20000000000
// true       false true       true     false      99/20000000000
// true       false true       false    true       70929/100000000000
// true       false true       false    false      7021971/100000000000
// true       false false      true     true       1497/50000000000
// true       false false      true     false      148203/50000000000
// true       false false      false    true       498002499/1000000000000
// true       false false      false    false      49302247401/1000000000000
// false      true  true       true     true       133/1000000000
// false      true  true       true     false      57/1000000000
// false      true  true       false    true       202797/5000000000
// false      true  true       false    false      86913/5000000000
// false      true  false      true     true       164171/2500000000
// false      true  false      true     false      70359/2500000000
// false      true  false      false    true       3489507/50000000000
// false      true  false      false    false      1495503/50000000000
// false      false true       true     true       19/20000000000
// false      false true       true     false      1881/20000000000
// false      false true       false    true       1347651/100000000000
// false      false true       false    false      133417449/100000000000
// false      false false      true     true       28443/50000000000
// false      false false      true     false      2815857/50000000000
// false      false false      false    true       9462047481/1000000000000
// false      false false      false    false      936742700619/1000000000000"""

Variables can be summed out of the factor:

import axle._

jpt.sumOut(M).sumOut(J).sumOut(A).sumOut(B).sumOut(E)
// res3: Factor[Boolean, Rational] = Factor(
//   variablesWithValues = Vector(),
//   probabilities = Map(Vector() -> 1)
// )

Also written as:

jpt.Σ(M).Σ(J).Σ(A).Σ(B).Σ(E)
// res4: Factor[Boolean, Rational] = Factor(
//   variablesWithValues = Vector(),
//   probabilities = Map(Vector() -> 1)
// )

Multiplication of factors also works:

import spire.implicits.multiplicativeSemigroupOps

val f = (bn.factorFor(A) * bn.factorFor(B)) * bn.factorFor(E)
f.show
// res5: String = """Burglary Earthquake Alarm
// true     true       true  19/10000000
// true     true       false 1/10000000
// true     false      true  23453/25000000
// true     false      false 1497/25000000
// false    true       true  28971/50000000
// false    true       false 70929/50000000
// false    false      true  498501/500000000
// false    false      false 498002499/500000000"""

Markov assumptions:

bn.markovAssumptionsFor(M).show
// res6: String = "I({Mary Calls}, {Alarm}, {John Calls, Earthquake, Burglary})"

Future Work

This is read as "M is independent of E, B, and J given A".