The BayesianNetwork paclet

Andreas Hafver, DNV

#WolframTechConf

ABSTRACT

This presentation describes an implementation of Bayesian Networks (BNs) in the Wolfram Language. BNs are graphical probabilistic models for compact representation of multivariate categorical probability distributions. They can be used for reasoning under uncertainty and to evaluate the effects of evidence on beliefs, for example, in decision-support systems. Conveniently, the Wolfram Language allows one to specify BNs with either numerical, symbolic or mixed probabilities. This makes it possible to obtain symbolic answers to BN queries, which is useful for performing sensitivity analysis in cases where precise probability values cannot be assigned.

Outline

About me

A brief introduction to probabilistic modelling and Bayesian Networks

The BayesianNetwork paclet

Some implementation details

Concluding remarks and thoughts

About me

Physics education

Bachelor and Master degrees in Theoretical Physics, University of Stellenbosch

PhD in Computational Physics/Geoscience, University of Oslo

Attended the Wolfram Science Summer School 2013

Working as a researcher in DNV (since 2014)

Current research focus

Foundational issues in risk management and assurance

Assurance of AI-enabled systems

Trustworthy AI in industrial contexts

Emerging technologies (e.g., quantum technologies, scientific machine learning, Tsetlin machines)

Probabilistic modelling and Bayesian Networks

Probability theory

There are different interpretations of probabilities:

Frequentist: Probabilities represent a physical tendency of something to occur, measured as relative frequencies if a situation is repeated many times.
(Key figures: von Mises,  Fisher, Pearson)

Bayesian/subjective: Probabilities represent subjective degrees of belief in outcomes. 
(Key figures: Bayes, Laplace, Finetti,Savage)

Logical: Probability is an extension of logic and a tool for reasoning under uncertainty. 
(Key figures: Jaynes)

Irrespective of interpretation, probability theory provides a framework for reasoning under uncertainty:

Axiom 1: For any outcome A, P(A)  ≥ 0

Axiom 2: P(S) = 1, where S is the space of all possible outcomes.

Axiom 3: P(A∪B) = P(A) + P(B) if A and B are mutually exclusive.

Joint and conditional probabilities can be used to express when outcomes occur together:

Joint probability of A and B:    P(A∩B)
	
	Conditional probability of A given B:  P(A|B)= ((P(A∩B))/(P(B)))

Bayes theorem let’s us reverse the conditioning and update beliefs based on observations:
	
	     P(A|B)= ((P(B|A)P(A))/(P(B)))

Note: All probabilities are conditional (assigned based on choices and assumptions)

Different extensions of probability theory attempt to circumvent this by adding “more dimensions” to probabilities (Dempster-Shafer theory,  imprecise probabilities, subjective logic, fuzzy logic, etc.)

My preference: Stick to probabilities for consistent reasoning (inspired by Jaynes), but identify and assess the impact of your assumptions (inspired by Aven’s concepts of “knowledge strength” and “assumption-deviation risk”)

Bayesian networks

In practical applications, we are often interested in multivariate probability distributions:

Notation: P(A,B,C) = P(A∩B∩C)

We can factor multivariate probabilities, e.g., P(A,B,C) = P(A)P(B)P(C|A,B) and represent it graphically

There are many ways to factor a distribution, which can be represented graphically:

Grid[{{"P(A,B,C)=P(A)P(B)P(C|A,B)","P(A,B,C)=P(A)P(B|A)P(C|B)","P(A,B,C)=P(A)P(B|A)P(C|A)"},{,,}},Frame->All]

Note: The factorization need not reflect causal dependencies, but they often will when we build a model.

Bayesian networks (BNs) offer a compact representation of multivariate distributions by factoring them as directed acyclic graphs (DAGs).

For multivariate categorical probability distributions, the number of state combinations quickly becomes unmanageable as the number of variables increases.

Factorization (exploiting known dependencies and independencies) can greatly reduce the number of parameters needed to specify a joint distribution.

Each node in a BN represent a variable with an associated conditional probability table, for example:
           ({{({{, "X=True", "X=False"}, {"P(X)", 0.9`, 0.1`}})}, {({{, "Y=Early", "Y=Late"}, {"P(Y)", 0.8`, 0.2`}})}, {({{, "Z=Low", "Z=Medium", "Z=High"}, {"P(Z|X=True,Y=Early)", 0.7`, 0.2`, 0.1`}, {"P(Z|X=True,Y=Late)", 0.5`, 0.3`, 0.2`}, {"P(Z|X=False,Y=Early)", 0.3`, 0.4`, 0.3`}, {"P(Z|X=False,Y=Late)", 0.1`, 0.2`, 0.7`}})}})

By applying Bayes theorem, one can update beliefs about variables in the network given observations of other variables.

Note: Bayesian networks are usually restricted to categorical variables.

We can deal with continuous variables through discretization.

There are generalizations of BNs to continuous variables, but they rely on assumptions (e.g. Gaussian distributions and linear dependencies)

Probabilities in the Wolfram language

The WL has extensive support for probabilistic modelling.

The WL has a categorical distribution

Univariate distributions (variable X has outcomes A,B,C):

uni=CategoricalDistribution[{"A","B","C"},{0.3,0.2,0.5}]

RandomVariate[uni,10]

PDF[uni,"A"]

Probability[x=="A",x\[Distributed]uni]

Multivariate distributions (variable X has outcomes A,B,C, and variable Y has outcomes D,E):

multi=CategoricalDistribution[{{"A","B","C"},{"D","E"}},{{(1/21),(2/21)},{(1/7),(4/21)},{(5/21),(2/7)}}]

RandomVariate[multi,10]

PDF[multi,{"A","E"}]

Probability[x=="A"\[Conditioned]y=="E",{x,y}\[Distributed]multi]

MarginalDistribution[multi,1]

However, WL lacks BN functionality...

CategoricalDistribution only has a notion of outcome names, but not variable names for multivariate distributions

You cannot specify ConditionalDistribution in factorized form as products of other distributions

ProbabilityDistribution can create distributions from other distributions, but does not support products of CategoricalDistribution

The BayesianNetwork paclet

My motivations for making the paclet:

I have wished for many years that Wolfram would add BN support, but it has not happened yet.

I wanted to learn how to make a paclet

My basic design principles:

Imagine it was already part of WL, how would I want to use it with existing functions like Probability?

Try to make it “feel and look WL native”

Simple example usage:

The Monty Hall problem:

You are presented with three doors...
Behind one door is a car (the prize you want), and behind the other  two doors are goats...
You choose one of the doors....
The game show host, who knows what’s behind each door, opens another door with a goat behind it...
You are now given a choice: Stick with your original choice, or switch to the other unopened door?

We can formulate this as a BN:

s = Function[CategoricalDistribution[{"A","B","C"},{1/3,1/3,1/3}]];
c=Function[CategoricalDistribution[{"A","B","C"},{1/3,1/3,1/3}]];
o = Function[CategoricalDistribution[{"A","B","C"},
Boole[#Selected=="A"&&#Car=="A"]{0,1/2,1/2}
+Boole[#Selected=="A"&&#Car=="B"]{0,0,1}
+Boole[#Selected=="A"&&#Car=="C"]{0,1,0}
+Boole[#Selected=="B"&&#Car=="A"]{0,0,1}
+Boole[#Selected=="B"&&#Car=="B"]{1/2,0,1/2}
+Boole[#Selected=="B"&&#Car=="C"]{1,0,0}
+Boole[#Selected=="C"&&#Car=="A"]{0,1,0}
+Boole[#Selected=="C"&&#Car=="B"]{1,0,0}
+Boole[#Selected=="C"&&#Car=="C"]{1/2,1/2,0}
]];

BN = BayesianNetwork[{
<|"Name"->"Selected","Distribution"->s|>,
<|"Name"->"Car","Distribution"->c|>,
<|"Name"->"Opened","Distribution"->o|>
}]

We can plot it:

Graph[BN,VertexLabels->"Probabilities"]

We can also get some information about BN:

Information[BN]

Suppose you select door A and the host opens door C. What should you do?

Probability[#Car=="A"\[Conditioned]#Selected=="A"&&#Opened=="C"&,BN]

Probability[#Car=="B"\[Conditioned]#Selected=="A"&&#Opened=="C"&,BN]

We can also plot it :

Graph[BN1,#Selected=="A"&&#Opened=="C"&,VertexLabels->"Probabilities"]

Advanced functionality for free in WL

“Wet grass” example with numeric probabilities:

c = Function[CategoricalDistribution[{"True","False"},{0.5,0.5}]];
r= Function[CategoricalDistribution[{"True","False"},Boole[#Cloudy=="True"]{0.8,0.2}+Boole[#Cloudy=="False"]{0.2,0.8}]];
s=Function[CategoricalDistribution[{"On","Off"},Boole[#Cloudy=="True"]{0.1,0.9}+Boole[#Cloudy=="False"]{0.5,0.5}]];
g = Function[CategoricalDistribution[{"Wet","Dry"},Boole[#Rain=="True"&&#Sprinkler=="On"]{0.99,0.01}+ Boole[#Rain=="False"&&#Sprinkler=="On"]{0.9,0.1}+Boole[#Rain=="True"&&#Sprinkler=="Off"]{0.9,0.1}+Boole[#Rain=="False"&&#Sprinkler=="Off"]{0.0,1.0}]];

BN2 = BayesianNetwork[{
<|"Name"->"Cloudy","Distribution"->c|>,
<|"Name"->"Rain","Distribution"->r|>,
<|"Name"->"Sprinkler","Distribution"->s|>,
<|"Name"->"Grass","Distribution"->g|>
}]

Graph[BN2,VertexLabels->"Probabilities"]

Graph[BN2,#Cloudy=="True"&,VertexLabels->"Probabilities"]

“Wet grass” example with mixed probabilities

BN3= BayesianNetwork[{
<|"Name"->"Cloudy","Distribution"->Function[CategoricalDistribution[{"True","False"},{0.5,0.5}]]|>,
<|"Name"->"Rain","Distribution"->Function[CategoricalDistribution[{"True","False"},Boole[#Cloudy=="True"]{p,1-p}+Boole[#Cloudy=="False"]{q,1-q}]]|>,
<|"Name"->"Sprinkler","Distribution"->Function[CategoricalDistribution[{"On","Off"},Boole[#Cloudy=="True"]{0.1,0.9}+Boole[#Cloudy=="False"]{0.5,0.5}]]|>,
<|"Name"->"Grass","Distribution"->Function[CategoricalDistribution[{"Wet","Dry"},Boole[#Rain=="True"&&#Sprinkler=="On"]{0.99,0.01}+ Boole[#Rain=="False"&&#Sprinkler=="On"]{0.9,0.1}+Boole[#Rain=="True"&&#Sprinkler=="Off"]{0.9,0.1}+Boole[#Rain=="False"&&#Sprinkler=="Off"]{0.0,1.0}]]|>
}];

Graph[BN3,VertexLabels->"Probabilities"]

“Wet grass” example with all symbolic probabilities

BN4 = BayesianNetwork[{
<|"Name"->"Cloudy","Distribution"->Function[CategoricalDistribution[{"True","False"},{pC,1-pC}]]|>,
<|"Name"->"Rain","Distribution"->Function[CategoricalDistribution[{"True","False"},Boole[#Cloudy=="True"]{pRc,1-pRc}+Boole[#Cloudy=="False"]{pRnc,1-pRnc}]]|>,
<|"Name"->"Sprinkler","Distribution"->Function[CategoricalDistribution[{"On","Off"},Boole[#Cloudy=="True"]{pSc,1-pSc}+Boole[#Cloudy=="False"]{pSnc,1-pSnc}]]|>,
<|"Name"->"Grass","Distribution"->Function[CategoricalDistribution[{"Wet","Dry"},Boole[#Rain=="True"&&#Sprinkler=="On"]{pGrs,1-pGrs}+ Boole[#Rain=="False"&&#Sprinkler=="On"]{pGnrs,1-pGnrs}+Boole[#Rain=="True"&&#Sprinkler=="Off"]{pGrns,1-pGrns}+Boole[#Rain=="False"&&#Sprinkler=="Off"]{pGnrns,1-pGnrns}]]|>
}]

Probability[#Grass =="Wet"\[Conditioned] #Cloudy == "False"&,BN4]

We can do more complicated queries:

Probability[#Grass =="Wet"∧#Cloudy=="False"\[Conditioned] #Rain =="True"∨#Sprinkler == "On"&,BN4]

Implementation details

Key issues I needed to solve:

Making Probability work with named variables

The  Probability function in the WL works with multivariate distributions, however one has to remember the order of the variables, e.g., {x, y, z} \[Distributed] dist.

To overcome this, my BayesianNetwork keeps track of named variables using #name (Slot[“name”]).

I then overloaded the Probability function so that it works with constructs like Probability[#var2 \[Conditioned] #var1 && #var3 &, BN].

Making ProbabilityDistribution work with categorical variables

To solve the second issue, I noted that there are two ways to implement categorical distributions in the WL:

dist1 = CategoricalDistribution[{"A","B","C"},{0.1,0.2,0.7}]

dist2 = EmpiricalDistribution[{0.1,0.2,0.7}->{"A","B","C"}]

Note that they behave differently when we try to retrieve their PDF functions:

PDF[dist1,x]

PDF[dist2,x]

EmpiricalDistribution directly produces a symbolic PDF function in terms of Boole . This is exactly what I need to be able to construct multivariate distributions in my BN model. However, the domain is still categorical, which ProbabilityDistribution doesn’t like.

I solved it by introducing a domain mapping, e.g. <|”A”->1,”B”->2,”C”->3|>.

My implementation of BayesianNetworkaccepts both CategoricalDistribution and EmpiricalDistribution as input, but represent all distributions internally as DataDistribution with integer range domains, e.g.:

Code for creating  a BayesianNetwork

(* Create Bayesian network model *)
BayesianNetwork[vars_?(VectorQ[#,AssociationQ]&)] := Module[{names,distributions,name2slot,a,parents,edges,domainmaps,transformationrules,pdf,ranges,distribution,x},
	(* Extract names and distributions *)
	names = vars[[All,"Name"]];
	distributions = vars[[All,"Distribution"]];
	name2slot = Association@@Table[names[[i]]->Slot[i],{i,Length@names}];
	
	(* Create BN graph*)
	parents = Table[DeleteDuplicates[Cases[dist,Slot[a_]->a,-1]],{dist,distributions}];
	edges = Flatten[Table[Thread[parents[[i]]->names[[i]]],{i,Length@names}]];
	
	(* Transform all distributions to integer domains *)
	domainmaps = Table[AssociationThread[#->Range@Length@#]&@DistributionDomain[dist[[1]]],{dist,distributions}];
	transformationrules = Flatten@Table[
		Module[{slot,dom},
			dom = domainmaps[[i]];
			slot = Slot@names[[i]];
			MapThread[Equal[slot,#1]->Equal[slot,#2]&,{Keys[dom],Values[dom]}]
		],
	{i,Length@distributions}];
	distributions = distributions/.transformationrules;
	
	(* Replace named variables by dummy variables *)
	distributions = Table[
		Module[{args,dist,dom,probs},
			args = Thread[name2slot[[parents[[i]]]]];
			dist = distributions[[i]][args];
			dom = DistributionDomain[dist];
			probs = Table[Evaluate[PDF[dist,y]],{y,dom}];
			EmpiricalDistribution[probs->Range@Length@dom]
		]
	,{i,Length@distributions}];	

	(* Create joint distribution *)
	pdf = Function[Evaluate[Times@@Table[PDF[distributions[[i]],Slot[i]],{i,Length@distributions}]]];
	ranges = Table[{x[i],Sequence@@MinMax[DistributionDomain[distributions[[i]]]],1},{i,Length@distributions}];
	distribution = ProbabilityDistribution[pdf@@Table[x[i],{i,Length@distributions}],Sequence@@ranges];

	(* Return BN object *)
	BayesianNetwork[<|"Variables"->Table[<|"Name"->names[[i]],"Parents"->parents[[i]],"Domain"->domainmaps[[i]],"Distribution"->distributions[[i]]|>,{i,Length@names}],"Distribution"->distribution,"Edges"->edges|>]
]

(* Test *)
BayesianNetworkQ[asc_?AssociationQ] := AllTrue[{"Variables","Distribution","Edges"}, KeyExistsQ[asc, #]&]

(* Visual appearance of BN objects *)
BayesianNetworkObject/:MakeBoxes[obj:BayesianNetworkObject[asc_?BayesianNetworkQ],form:(StandardForm|TraditionalForm)]:=
Module[{info},
	info={{BoxForm`SummaryItem[{Graph[asc["Edges"],VertexLabels->None,VertexSize->Large,ImageSize->50]}]}};
	BoxForm`ArrangeSummaryBox[BayesianNetworkObject,obj,None,info,{},form,"Interpretable"->Automatic]
]

(* Define associated information object *)
Information`AddRegistry[BayesianNetwork, getBayesianNetworkInformation];
getBayesianNetworkInformation[BayesianNetwork[asc_?BayesianNetworkQ]] := <|
	"ObjectType" -> "Bayesian Network", <||>, 
	"Variable count" :> Length[asc["Variables"]],
	"Variables"->asc[["Variables",All,"Name"]],
	"Variable domains":>Map[Keys,asc[["Variables",All,"Domain"]]],
	"Edge count":>Length[asc["Edges"]],
	"Edges"->asc["Edges"]
|>

Code for computing probabilities (overloading Probability)

(* Calculation of probabilities *)
Unprotect[Probability];
BayesianNetwork/: Probability[expr_Function,BayesianNetwork[asc_?BayesianNetworkQ],opts:OptionsPattern[Probability]]:= Module[{variables,variableNames,variableSlots,variableAssociation,dom,dist,x,rules,newexpr,slot},
	variables = asc["Variables"];
	variableNames = variables[[All,"Name"]];
	variableSlots = Table[x[i],{i,Length@variableNames}];
	variableAssociation = Association[Thread[variableNames->variableSlots]];

	(* Transform expression to the BN dummy variables with integer domains *)
	newexpr = expr; 
	Do[
		dom = variables[[i]]["Domain"];
		slot = Slot@variableNames[[i]];
		rules = MapThread[Equal[slot,#1]->Equal[slot,#2]&,{Keys[dom],Values[dom]}];
		newexpr = (newexpr/.rules),
	{i,Length@variables}];

	(* Compute the probability *)
	Probability[newexpr[variableAssociation],variableSlots \[Distributed] asc["Distribution"],opts]
]
Protect[Probability];

Code for plotting BNs (overloading Graph)

(* Plotting *)
Unprotect[Graph];
BayesianNetwork/: Graph[BayesianNetwork[asc_?BayesianNetworkQ],evidence_Function:Function[True],opts:OptionsPattern[Graph]]:= Module[{variables,edges,evidencenodes,a,colour,labels,probs,domain,outcomes,barcharts},
	variables = asc["Variables"];
	edges = asc["Edges"];
	
	If[OptionValue[FilterRules[{opts},VertexLabels],VertexLabels] ==="Probabilities",
		If[evidence===Function[True],
			evidencenodes = {},
			evidencenodes = DeleteDuplicates[Cases[evidence,Slot[a_]->a,-1]]
		];
		labels = Table[
			domain = {Slot[var["Name"]],Sequence@@MinMax[var["Domain"]],1};
			probs = Table[Probability[Evaluate[Slot[var["Name"]]==val\[Conditioned] evidence[[1]]]&,BayesianNetwork[asc]],{val,Keys[var["Domain"]]}];
			If[VectorQ[probs,NumericQ],
				var["Name"]->Labeled[BarChart[probs,PlotRange->{0,1},BarOrigin->Left,ImageSize->Tiny,Frame->True,FrameTicks->{None,Automatic},ChartLabels->Keys[var["Domain"]],ChartStyle->If[MemberQ[evidencenodes,var["Name"]],Blue,Orange]],var["Name"],Top],
				var["Name"]->Labeled[TableForm[probs,TableHeadings->{Keys[var["Domain"]],None}],var["Name"],Top]
			],
			{var,variables}];
			Graph[edges,VertexLabels->labels,opts],
			Graph[edges,opts]
		]
]
Protect[Graph];

That’s all the code I needed to implement BNs in the WL!

Concluding remarks

Advantages of implementing BN in WL

We can use existing inbuilt probability functionality, so very little new code is required.

WL naturally supports numeric and symbolic probabilities (nice if you want to explore the effects of changing probability assignments)

We can do arbitrarily complicated queries in Probability.

Limitations of this paclet

It is slow for larger networks (because I have not implemented any of the algorithms used for effective BN inference in other languages).

The code for achieving the variable transformations feels a bit ‘hacky’ and could maybe be done more elegantly.

Possible next steps

Implementing efficient BN inference algorithms

Adding support for importing common BN file formats

Adding support for discretization of continuous variables

Thank you!

Get the paclet:

https://resources.wolframcloud.com/PacletRepository/resources/AndreasHafver/BayesianNetwork

https://github.com/ahafver/BayesianNetwork

Get in touch with me:

Andreas.Hafver@dnv.com

https://www.linkedin.com/in/andreas-hafver/

https://community.wolfram.com/web/andreashafver/