# Probabilistic Graphical Models

Currently only Bayesian Networks

Eventually others including Pearl's causal models.

## 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()``````

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".

• Test: start with `ABE.jointProbabilityTable` (monotype `tuple5[Boolean]`)
• Factor out each variable until original 5-note network is reached
• Basically the inverse of factor multiplication
• `bn.factorFor(B) * bn.factorFor(E)` should be defined? (It errors)
• `MonotypeBayesanNetwork.filter` collapase into a single BNN

• Review `InteractionGraph`, `EliminationGraph`, `JoinTree` and the functions they power

• Consider a "case" to be a `Map` vs a `Vector`

• Consider usefulness of `Factor` in terms of `Region`

• `MonotypeBayesanNetwork`.{`pure`, `map`, `flatMap`, `tailRecR`}

• Reconcile `MBN` `combine1` & `combine2`

• Monad tests for `MonotypeBayesanNetwork[Alarm-Burglary-Earthquake]`

• `Bayes[MonotypeBayesanNetwork]` -- could be viewed as "belief updating" (vs "conditioning")

• If it took a ProbabilityModel itself