Approximating Joint Distributions by Inferring Their Causal Structure

July 2013

Project Description

Suppose we are interested in a probabilistic system characterized by some set of random variables {(Subscript[X, 1]), (Subscript[X, 2]),..., (Subscript[X, n])}.The joint distribution P((Subscript[X, 1]), (Subscript[X, 2]),...,(Subscript[X, n])) is the most complete available representation of the probabilistic system, as all other probabalistic quantities of interest can be obtained from it.

The project will consist of two parts. First, Mathematica lacks certain built-in functionality with respect to manipulating discrete probability distributions. For example, it lacks a CategoricalDistribution or ConditionalDistribution function. The former is the most general distribution over a random variable with K non-ordered outcomes. The latter takes in a distribution and condition and returns the conditioned distribution. These functions will be written.

The second part of the project deals with the major computational problem of using joint distributions: the exponential growth in memory usage with increase in joint distribution dimension. For high-dimensional distributions, some approximation scheme for the joint distribution is needed. Two general classes of non-parametric approximations exist in the literature: spectral methods and causal decompositions (or Bayesian networks). In the project, we will implement the causal decomposition approach, where the causal structure of the random variables can be inferred, leading to a more efficient representation of the probabilistic system.

Code

This code needs to be run.

Functions Defined in Project

IndexPosition::usage="IndexPosition[x_,variables_]: variables_ is a list of variable names. x is of form x ⊆ variables. 
												Returns the position of the elements of x in the list variables_";

ContingencyTableDistribution::usage= "ContingencyTableDistribution[variables_,ContTab_]: takes in a list of categories ('variables_')\\n					and a normalized contingency table ('ContTab_'). Returns an unevaluated ContingencyTableDistribution object.";

ContingencyNormalize::usage= "ContingencyNormalize[ContTab]: takes in an unormalized contingency table (ContTab) and returns a normalized one";

EqualityContingencyTable::usage= "EqualityContingencyTable[ContingencyTableDistribution[...],equalities_]: Takes in ContingencyTableDistribution oject and a list of 
								equalities_ = {{'sex',1},...}. Returns ContingencyTableDistribution oject.";

MarginalDistribution::usage= "MarginalDistribution[ContingencyTableDistribution[variables_,CrossTab_],margs_]: takes in a ContingencyTableDistribution 
							and a list margs_ ={'sex',...} of variables to remain after marginalization. It returns a ContingencyTableDistribution 
							with variables = margs_. Extends built-in MarginalDistribution[] function to work with ContingencyTableDistribution objects.";

ConditionalDistribution::usage= "ConditionalDistribution[ContingencyTableDistribution[...],CondVar_]: Takes in a ContingencyTableDistribution object and list of variables
								to condition over (CondVar_), returns ContingencyTableDistribution object representing conditional distribution P(not in CondVar | CondVar ).";

Probability::usage = "Probability[pred_,vars_List \[Distributed] ContingencyTableDistribution[variables_, crossTab_]]: takes in a list of variables vars={x,y,...} which assign variables 
					to same order as variables in (variables_), and some logical condition on variables (pred_). Returns a probability.";

Entropy::usage = "Entropy[ContingencyTableDistribution[...]]: extends built in entropy function to compute cross entropy of ContingencyTableDistribution object.";

ConditionalMutualInformation::usage= "ConditionalMutualInformation[X,Y,Z_List,ContingencyTableDistribution[...]]: calculates P(X,Y | Z) where X,Y,Z_List are 
									category names in ContingencyTableDistribution. Z_List can be the empty list. ";

CITest::usage = "CITest[pvalue_,samplesize_,X_, Y_, Z_List, ContingencyTableDistribution_]: Returns True if X and Y are independent given Z using the (G^2)-test with some p-value";

PCAlgorithm::usage="PCAlgorithm[pvalue,sampleSize,contingencyTableDistribution]: 
		Returns a list of {complete graph of variables, skeleton of causal graph, equivalence class of graphs with red lines of undetermined direction, member of equivalence
		class}.";

Basic Joint Distribution Functions

ContingencyTableDistribution

ContingencyNormalize[CrossTab_]:=CrossTab/Total[CrossTab,Infinity]

ContingencyTableDistribution[variables_,CrossTab_] := ContingencyTableDistribution[variables, CrossTab, dummy];

Plus method for two ContingencyTableDistributions. They must have the same variables.

ContingencyTableDistribution/:Plus[ContingencyTableDistribution[variables_,CrossTab1_],ContingencyTableDistribution[variables_,CrossTab2_]]:=ContingencyTableDistribution[variables,CrossTab1+CrossTab2]

IndexPosition[userVar_,PDFvariables_]:=Position[PDFvariables,#][[1,1]]&/@userVar

Equality

ClearAll[EqualityContingencyTable];
EqualityContingencyTable[ContingencyTableDistribution[variables_,CrossTab_],equalities_]:=Module[{positions,default,newVariables,newCont},
	positions=IndexPosition[Transpose[equalities][[1]],variables];
	default=ConstantArray[All,Length@variables];
	default[[positions]]=Transpose[equalities][[2]];
	newVariables=variables[[Complement[Range@Length@variables,positions]]];
	newCont=CrossTab[[Sequence@@default]];
	ContingencyTableDistribution[newVariables,newCont,dummy]
]

Marginal Distribution

Clear[BringToFront];
BringToFront[matrix_, depth_, frontdims_, names_] := With[
	{perm = FindPermutation[Join[frontdims, Complement[Range[depth], frontdims]]]},
	{Transpose[matrix, perm], Permute[names,PermutationList@perm]}
]

(* marginalizeVar returns the variables in distr not in margs, as these must be marginalized over *)
ContingencyTableDistribution /: MarginalDistribution[ContingencyTableDistribution[variables_,CrossTab_],varToKeep_]:=
	Module[{newcross,marginalVar,front,newVariables},
	marginalVar=Complement[variables,varToKeep]; (* variables that we are marginalizing over *)
	front=BringToFront[CrossTab,Length@variables,IndexPosition[marginalVar,variables],variables];
	newcross=Total[front[[1]],{1,Length@marginalVar}]; 
	newVariables=front[[2]][[Length@marginalVar+1;;]];
	ContingencyTableDistribution[newVariables,newcross,dummy]
	]

Conditional Distribution

We use P((Subscript[X, 1]),...,(Subscript[X, n])|(Subscript[Z, 1]),...(Subscript[Z, m]))=P((Subscript[X, 1]),...,(Subscript[X, n]),(Subscript[Z, 1]),...,(Subscript[Z, m]))/P((Subscript[Z, 1]),...,(Subscript[Z, m]))

Clear[ConditionalDistribution];
ConditionalDistribution[ContingencyTableDistribution[variables_,CrossTab_],CondVar_]:=
	Module[{marginalDist,newjoint,rules,bringfront,newDist},
	marginalDist=MarginalDistribution[ContingencyTableDistribution[variables,CrossTab],CondVar]; (* marginal distribution *)
	bringfront=BringToFront[CrossTab,Length@variables,IndexPosition[marginalDist[[1]],variables],variables]; 
	newDist=bringfront[[1]]/marginalDist[[2]];
	ContingencyTableDistribution[bringfront[[2]],newDist,dummy]
	]

Probability Function

(* Helper Function *)
Combinations[vars_List,condition_,tableCross_]:=Module[{DimSize,ranges,allCombinations,select1,select2,validComb},
	DimSize=Dimensions[tableCross];
	ranges=Range[#]&/@DimSize;
	allCombinations=Partition[Flatten@Outer[List,Sequence@@ranges],Length@DimSize];
	(* selection function *)
	select1[Sequence@@(Pattern[#,_]&/@vars)]=condition;
	select2[list_List]:=select1[Sequence@@list];
	validComb=Select[allCombinations,select2]
	]

(* Probability Function *)
ContingencyTableDistribution;
Probability;
Unprotect[Probability, Distributed];

Probability[pred_,vars_List \[Distributed] ContingencyTableDistribution[variables_, crossTab_]] :=
	Module[{validComb1,validComb2,validComb3},
	If[TrueQ[Head[pred]==Conditioned],
		validComb1=Combinations[vars,pred[[1]],crossTab];
		validComb2=Combinations[vars,pred[[2]],crossTab];	
		validComb3=Intersection[validComb1,validComb2];
		Total[crossTab[[Sequence@@#]]&/@validComb3]/Total[crossTab[[Sequence@@#]]&/@validComb2]
		,
		validComb1=Combinations[vars,pred,crossTab];
		Total[crossTab[[Sequence@@#]]&/@validComb1]
	]
]
Protect[Probability, Distributed];

Test for conditional Independence: (G^2)

Description

Basic idea. The PC algorithm needs as input whether P(X,Y| (Subscript[Z, 1]),(Subscript[Z, 2])...(Subscript[Z, n])) is conditionally independent. One such test is the (G^2) test. It involves calculating the conditional mutual information:
I(X,Y| (Subscript[Z, 1]),(Subscript[Z, 2])...(Subscript[Z, n]))  = ∑ p((Subscript[z, 1]))...p((Subscript[z, n])) p(x,y|(Subscript[z, 1]),(Subscript[z, 2])...(Subscript[z, n])) log[p(x,y| (Subscript[z, 1])...(Subscript[z, n])) /p(x|(Subscript[z, 1])...(Subscript[z, n])) p(y| (Subscript[z, 1])...(Subscript[z, n]))].

Tricky to program. Much easier seems to be to use identity: I( X, Y | Z )  = H(X,Z)+H(Y,Z)- H(X,Y,Z)-H(Z). Generalizing this identity, we find
I( X, Y | (Subscript[Z, 1])...(Subscript[Z, n]) ) = H(X,(Subscript[Z, 1])...(Subscript[Z, n]))+H(Y,(Subscript[Z, 1])...(Subscript[Z, n]))- H(X,Y,(Subscript[Z, 1])...(Subscript[Z, n]))-H((Subscript[Z, 1])...(Subscript[Z, n]))

(Is using the definition directly faster? Or simplifying definition to joint distributions faster?)

Standard frequentist paradigm: calculate test statistic and compare against distribution of test statistic ((χ^2) in our case). Calculate p-value and accept or reject.

Assumptions

Large enough sample sizes: According to the Handbook of Biological Statistics: “When the sample sizes are too small, alternatives to the chi-square test or G-test are recommended. However, how small is “too small”? The conventional rule of thumb is that if all of the expected numbers are greater than 5, it’s acceptable to use the chi-square or G-test; if an expected number is less than 5, you should use an alternative, such as an exact test or a randomization test for goodness-of-fit, or a Fisher’s exact test or randomization test of independence.”

Note: We are ignoring this issue (and other issues like Yates continuity correction) for now!

Function Definitions

First add functionality to existing entropy function in order to calculate H((Subscript[X, 1]),(Subscript[X, 2]),...,(Subscript[X, n])):

ContingencyTableDistribution /: Entropy[ContingencyTableDistribution[variables_, CrossTab_]] := 
			If[Length@CrossTab==0,0,-Total[CrossTab Log[N[CrossTab]/. {0. -> 1}], Infinity]]

ConditionalMutualInformation[X_, Y_, Z_List, contDist_] := Entropy@MarginalDistribution[contDist, {X, Sequence@@Z}] + 
			Entropy@MarginalDistribution[contDist, {Y, Sequence@@Z}] - Entropy@MarginalDistribution[contDist, {X, Y, Sequence@@Z}] - 
			Entropy@MarginalDistribution[contDist, {Sequence@@Z}]

Needs["HypothesisTesting`"]
(* returns a p-value *)
GtestConditionalIndependence[samplesize_,X_, Y_, Z_List, ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{Gstatistic,dof,dimensions},
	Gstatistic= 2 samplesize ConditionalMutualInformation[X, Y, Z, ContingencyTableDistribution[variables, CrossTab]];
	(* We needs dimensions of variables X,Y,Z_List for degrees of freedom *)
	dimensions=Dimensions[CrossTab][[IndexPosition[{X,Y,Sequence@@Z},variables]]];
	dof=(dimensions[[1]]-1)(dimensions[[2]]-1) Times@@dimensions[[3;;]];
	ChiSquarePValue[Gstatistic,dof]
	]

CITest[pvalue_,samplesize_,X_, Y_, Z_List, ContingencyTableDistribution_]:=
		If[GtestConditionalIndependence[samplesize,X, Y, Z, ContingencyTableDistribution][[2]]>pvalue,True,False]

Test for conditional Independence: Bayesian Dirichlet

BayesianDirichletScore[x_,parents_,ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{s,sp,spp,marginalParents,marginalParentsX,term1,term2,term3},
	s=1; (* global sample size *)
	sp=1/Times@@Part[Dimensions@CrossTab,IndexPosition[parents,variables]];
	spp=sp/Times@@Part[Dimensions@CrossTab,IndexPosition[{x},variables]];
	marginalParents=MarginalDistribution[ContingencyTableDistribution[variables, CrossTab],parents];
	marginalParentsX=MarginalDistribution[ContingencyTableDistribution[variables, CrossTab],Join[parents,{x}]];
	term1=N@Log[((Gamma[sp])/(Gamma[spp]))];
	term2=If[Length@parents==0,
			0,-N@Total[Log[Gamma[marginalParents[[2]]+sp]],Infinity]
			];
	term3=N@Total[Log[Gamma[marginalParentsX[[2]]+spp]],Infinity];
	term1+term2+term3
	]

BayesianIndependence[x_,y_,z_,contingencyDist_]:=Module[{score1,score2},
	score1=BayesianDirichletScore[x,z,contingencyDist];
	score2=BayesianDirichletScore[x,Join[z,{y}],contingencyDist];
	If[score1>score2,True,False]
	]

PC Algorithm

Description

Algorithm Description: We use the description on page 117 from the book “Causation, Prediction, and Search” by P. Spirtes et al. The algorithm given finds an equivalence class of Bayesian Networks. Once the equivalence class is found, a particular member of the equivalence class can be found (useful for factorization).

Oher useful links: See Kalisch et. al. Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm, 2007. Also useful: Lecture notes S. Moral. Moral et al http://www.utia.cas.cz/files/mtr/pgm06/41_paper.pdf

This algorithm has three major sub-pieces:

A test for conditional independence. We will use the (G^2) test (we could also use the Pearson (χ^2) test, though it seems that the (G^2) test is to be preffered).

An algorithm that takes in an unstructured graph with every node connected to every other node (obtained from joint distribution), and trims the edges using the conditional independence test.

Finding the directions of the undirected edges.

Assumptions

acyclicity

causal sufficiency

independent sources of variation

Faithfulness: The joint distribution has all of the conditional independence relations implied by the causal Markov property, and only those conditional independence relations.

“The point of the faithfulness condition is to rule out “conspiracies among the parameters”, where, say, two causes of a common effect, which would typically be dependent conditional on that effect, have their impact on the joint effect and their own distributions matched just so exactly that they remain conditionally independent.” Shalizi

The local Markov assumption states that each node (Subscript[X, i]) is independent of its nondescendants given its parents.

Helper Graph Functions

CompleteContingencyGraph[ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{rules},
	rules={};
	Do[AppendTo[rules,Rule[i,ToString@variables[[i]]~~","~~ToString[i]]],{i,Length@variables}];
	CompleteGraph[Length@variables,VertexLabels-> rules]
	]

Define a function that returns False iff |(Subscript[ADJ, X])| ≤ i for all X: Used in next subsection

AdjQ[graph_,i_]:=Module[{vertices,adj},
	vertices=VertexList[graph];
	adj=AdjacencyList[graph,#]&/@vertices;
	If[Max[Length[adj],{2}]≤i,False,True]
	]

Part 1: Trimming Graph

GraphTrim1[pvalue_,sampleSize_,numConditioning_,X_,Y_,graph_,sxy_,adjX_,ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{sxy2,subsets,sub,graph2},
	subsets=Subsets[Select[adjX,#≠Y&],{numConditioning}]; (* All subsets of (AdjX-Y) of cardinality numConditioning *)
	graph2=graph; (* needs to be mutated *)
	sxy2=sxy;(* needs to be mutated *)
	Do[sub=subsets[[nn]];(* select nn subset *)
		If[CITest[pvalue,sampleSize,variables[[X]], variables[[Y]], variables[[sub]],ContingencyTableDistribution[variables, CrossTab]],
			sxy2[[Sequence@@{X,Y}]]=sub;sxy2[[Sequence@@{Y,X}]]=sub; (* add subset to sxy *)
			If[EdgeQ[graph2,X\[UndirectedEdge]Y],graph2=EdgeDelete[graph2,X\[UndirectedEdge]Y]]; (* Delete edge: but first make sure that there is an edge. *)
			Break (* exit Do loop *)
			],
		{nn,1,Length[subsets]} (* sum over all subsets *)
		];
	{graph2,sxy2}
	]

GraphTrim2[pvalue_,sampleSize_,numConditioning_,X_,graph_,sxy_,adjX_,ContingencyTableDistribution_]:=Module[{Y,graph2,sxy2},
	graph2=graph;(* needs to be mutated *)
	sxy2=sxy;(* needs to be mutated *)
	Do[Y=adjX[[nn]];{graph2,sxy2}=GraphTrim1[pvalue,sampleSize,numConditioning,X,Y,graph2,sxy2,adjX,ContingencyTableDistribution],
		{nn,1,Length[adjX]}];
	{graph2,sxy2}
	]

GraphTrim3[pvalue_,sampleSize_,numConditioning_,graph_,sxy_,ContingencyTableDistribution_]:=Module[{X,vertices,graph2,sxy2,adjacency},
	graph2=graph;(* needs to be mutated *)
	sxy2=sxy;(* needs to be mutated *)
	vertices=VertexList[graph2];
	Do[X=vertices[[nn]];adjacency=AdjacencyList[graph2,X];{graph2,sxy2}=GraphTrim2[pvalue,sampleSize,numConditioning,X,graph2,sxy2,adjacency,ContingencyTableDistribution],
		{nn,1,Length@vertices}];
	{graph2,sxy2}
	]

GraphTrim4[pvalue_,sampleSize_,ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{sxy,graph,numConditioning},
	graph=CompleteContingencyGraph[ContingencyTableDistribution[variables, CrossTab]];
	sxy=ConstantArray[{},{Length@variables,Length@variables}]; (* create an empty matrix to keep the values of sxy *)
	numConditioning=0;
	While[AdjQ[graph,numConditioning],{graph,sxy}=GraphTrim3[pvalue,sampleSize,numConditioning,graph,sxy,ContingencyTableDistribution[variables, CrossTab]];
		numConditioning++
		];
	{graph,sxy}
	]

Part 2: Finding directed edges

Note: Mathematica does not support mixed graphs!!! So we will keep track of the skeleton graph and directed graph seperately (maybe use GraphProperties?)

Idea Rule 0: Collider detection. Find X - Z  - Y (here X and Y are not adjacent). If Z ∉ (Subscript[S, XY]), then X - Z  - Y is a collider and hence X -> Z  ← Y.
One point: This procedure will always work IF no edges are missed when testing for conditional independence. One can get doubled edges if there are mistakes in the statistical tests. We will prevent any edge doubling.

HelperCollider[rulesSet_,rule_]:=If[Length@EdgeList[Graph[rulesSet],rule],True,False]

Colliders[graphDirected_,graphUndirected_,sxy_]:=Module[{subsets,rulesAdd,adj,graphDirected2,ruleGraph},
	rulesAdd={};
	graphDirected2=graphDirected;(* need to mutate *)
	Do[adj=AdjacencyList[graphUndirected,i]; (* loop through each vertex *)
		subsets=Select[Subsets[adj,{2}],Not@EdgeQ[graphUndirected,#[[1]]\[UndirectedEdge]#[[2]]]&]; (* find subsets of neighbours to vertex i that don't share edge *)	
		subsets=Select[subsets,Not@MemberQ[sxy[[Sequence@@#]],i]&]; (* for each pair (x,y) check whether i is in Sxy. Remove pair from subsets if it is. *)
		(* Statistical test failure: prevent doubling of edges. *)
		ruleGraph=UndirectedGraph@Graph@rulesAdd; 
		subsets=Select[subsets,((Not@EdgeQ[ruleGraph,i\[UndirectedEdge]#[[1]]])&&(Not@EdgeQ[ruleGraph,i\[UndirectedEdge]#[[2]]]))&]; 
		rulesAdd=Union[rulesAdd,Flatten[{#[[1]]\[DirectedEdge]i,#[[2]]\[DirectedEdge]i}&/@subsets]]; 
			,{i,Length@VertexList[graphUndirected]}];
	EdgeAdd[graphDirected2,rulesAdd]  (* returns a set of rules for directed edges *)
					]

The rest of the rules come from the fact that 1) all colliders are accounted for and 2) the graph must be acyclic. There are three seperate rules that are looped.

Idea Rule 1: Known non-colliders. We have already found all colliders in Rule 1. So if X -> Z  - Y then X -> Z -> Y (here X and Y are not adjacent, and Z has no other parents than X).

Rule1[graphSkeleton_,graphDirected_]:=Module[{adj,graphDirected2,edges,edge,graphDirected2temp},
	graphDirected2=graphDirected;
	edges=EdgeList[graphDirected2];
	Do[edge=edges[[i]];
		If[Length@EdgeList[graphDirected2,_\[DirectedEdge] edge[[2]]]==1, (* make sure edge[[2]] has no other parents than edge[[1]] *)
			adj=AdjacencyList[graphSkeleton,edge[[2]]]; (* all arrows point left to right *)
			adj=Select[adj,(Not@EdgeQ[UndirectedGraph@graphDirected2,edge[[2]]\[UndirectedEdge]#]&&Not@EdgeQ[graphSkeleton,edge[[1]]\[UndirectedEdge]#])&];
			Do[graphDirected2=EdgeAdd[graphDirected2,edge[[2]]\[DirectedEdge]adj[[nn]]],{nn,Length@adj}]
			];
			,{i,Length@edges}];
	graphDirected2
					]

Idea Rule 2: Cycle avoidance. Find all undirected edges X - Y and then ask if there exists a directed path from X to Y. If there is a path, then X -> Y

Rule2[graphSkeleton_,graphDirected_]:=Module[{graphDirected2,edges,edge},
	graphDirected2=graphDirected; (* to mutate *)
	edges=Select[EdgeList[graphSkeleton],(Not@EdgeQ[UndirectedGraph@graphDirected2,#[[2]]\[UndirectedEdge]#[[1]]])&];
	Do[edge=edges[[i]];
		If[GraphDistance[graphDirected2,edge[[1]],edge[[2]]]≠Infinity,
			graphDirected2=EdgeAdd[graphDirected2,edge[[1]]\[DirectedEdge]edge[[2]]],
			If[GraphDistance[graphDirected2,edge[[2]],edge[[1]]]≠Infinity,
				graphDirected2=EdgeAdd[graphDirected2,edge[[2]]\[DirectedEdge]edge[[1]]];
				]]
		,{i,Length@edges}];
	graphDirected2
	]

Find member of equivalence class: Fill in the red edges without 1) creating new colliders 2) create no new cycles.

numColliders[graph_]:=Total@((#[[2]]&/@Tally[#[[2]]&/@EdgeList@graph])-1)

There exist some configurations such that a collider must be added. This is a sure sign that something has gone wrong!

EquivalenceExample[missingEdges_,graphDirected_]:=Module[{i,j,directedEdges,numbercolliders,edgeList,graphDirected3},
	(* missing edges are undirected *)
	directedEdges={#[[1]]\[DirectedEdge]#[[2]],#[[2]]\[DirectedEdge]#[[1]]}&/@missingEdges;
	directedEdges=Partition[Flatten@Outer[List,Sequence@@directedEdges], Length@missingEdges]; (* all possible combinations of missing arrow directions *)
	numbercolliders=numColliders[graphDirected];
	i=1;j=1;
	While[i≤Length@directedEdges,
			edgeList=directedEdges[[i]];graphDirected3=EdgeAdd[graphDirected,edgeList];
			If[AcyclicGraphQ[graphDirected3],
				If[numColliders[graphDirected3]==numbercolliders,i=Length@directedEdges+1]
				]
			;i++;j++];
	graphDirected3
	]

Final PC algorithm.

PCAlgorithm[pvalue_, samplesize_, contingencyTableDistribution_] := Quiet @ iPCAlgorithm[pvalue, samplesize, contingencyTableDistribution]

iPCAlgorithm[pvalue_, samplesize_, ContingencyTableDistribution_]:=
	Module[{variables,rules,trimmedGraph,graphDirected,graphDirected2,missingEdges,equivalenceExample,equivalence},
	(* Step 1: Trimmed Graph *)
	trimmedGraph=GraphTrim4[pvalue,samplesize,ContingencyTableDistribution];
	(* Step 2: Create directed graph object with only vertices. Then find all colliders *)
	graphDirected=Graph[VertexList[trimmedGraph[[1]]],{},VertexLabels->"Name"];
	graphDirected2=Colliders[graphDirected,trimmedGraph[[1]],trimmedGraph[[2]]];
	(* Step 3: Apply rules 1 and 2 until no change in graph. Find equivalence class *)
	While[graphDirected≠graphDirected2,
		graphDirected=graphDirected2;
		graphDirected2=Rule1[trimmedGraph[[1]],graphDirected2];
		graphDirected2=Rule2[trimmedGraph[[1]],graphDirected2] 
			];
	missingEdges=(Complement[EdgeList[trimmedGraph[[1]]],EdgeList[UndirectedGraph@graphDirected]])/.UndirectedEdge-> DirectedEdge;
	variables=ContingencyTableDistribution[[1]];
	rules={}; (* rules to label graph *)
	Do[AppendTo[rules,Rule[i,ToString@variables[[i]]~~","~~ToString[i]]],{i,Length@variables}];
	graphDirected=Graph[VertexList[trimmedGraph[[1]]],
		Union[missingEdges,EdgeList@graphDirected2],EdgeStyle-> (Rule[#,Red]&/@missingEdges),VertexLabels-> rules];
	(* Step 4: find explicit example of member of equivalence class *)
	equivalence=EquivalenceExample[missingEdges,graphDirected2];
	equivalenceExample=Graph[VertexList@equivalence,EdgeList@equivalence,VertexLabels-> rules];
	(* return results *)
	{CompleteContingencyGraph[ContingencyTableDistribution],trimmedGraph[[1]],graphDirected,equivalenceExample}
	]

Generative Models Examples

Probabilistic Model 1

Model taken from page 15 of http://www.seas.upenn.edu/~taskar/pubs/gms-srl07.pdf. Some probabilities are not exactly the same, but probabilistic structure is.

fprob[p_]:=If[RandomReal[]<p,1,0];
RandomModel[a_]:=Module[{p,t,i,x,s},
p=fprob[0.5];
t=fprob[0.2];
If[p==1&&t==1,i=fprob[0.8],
If[p==1&&t==0,i=fprob[0.6],
If[p==0&&t==1,i=fprob[0.2],
If[p==0&&t==0,i=fprob[0.1]
	]]]];
If[i==1,x=fprob[0.8],x=fprob[0.6]];
If[t==1,s=fprob[0.8],s=fprob[0.6]];
{p,t,i,x,s}
]

RandomModelDistribution1[sample_]:=Module[{},
temp={};
Do[AppendTo[temp,RandomModel[1]],{sample}];
tal=Tally[temp];
array=ConstantArray[0,{2,2,2,2,2}];
(array[[Sequence@@(Part[#,1]+1)]]=#[[2]])&/@tal;
ContingencyTableDistribution[{"Pneumonia","Tuberculosis","Lung Infiltrates","X-ray","Sputum smear"},ContingencyNormalize[array]]
]

Probabilistic Model 2

http://www.ke.tu-darmstadt.de/lehre/archiv/ss09/ki/bayesian-networks1.pdf Page 24.
Probability assignments somewhat random, but structure is important.

fprob[p_]:=If[RandomReal[]<p,1,0];
RandomModel2[a_]:=Module[{s,va,b,lc,tb,toc,d,x},
va=fprob[0.2];
s=fprob[0.5];
If[va==1,tb=fprob[0.5],tb=fprob[0.1]];
If[s==1,lc=fprob[0.5],lc=fprob[0.1]];
If[s==1,b=fprob[0.3],b=fprob[0.02]];
If[tb==1∨lc==1,toc=1,toc=0];
If[toc==1,x=fprob[0.3],x=fprob[0.7]];
If[toc==1&&b==1,d=fprob[0.8],
	If[toc==1&&b==0,d=fprob[0.6],
		If[toc==0&&b==1,d=fprob[0.2],
			If[toc==0&&b==0,d=fprob[0.4]
	]]]];
{s,va,b,lc,tb,toc,d,x}
]

RandomModelDistribution2[sample_]:=Module[{},
temp={};
Do[AppendTo[temp,RandomModel2[1]],{sample}];
tal=Tally[temp];
array=ConstantArray[0,{2,2,2,2,2,2,2,2}];
(array[[Sequence@@(Part[#,1]+1)]]=#[[2]])&/@tal;
ContingencyTableDistribution[{"smoking","visit to asia","bronchitis","lung cancer","tuberculosis","tuberculosis or cancer","dyspnea","x-ray result"},ContingencyNormalize[array]]
]

Joint Distribution Approximation

The following function gives the (approximated) joint distribution of some graphical model in the form of a ContingencyTableDistribution. This takes exactly the same storage space as the non-approximated joint distribution, but this form is useful for testing the goodness of the approximation.

JointApproxFullForm[DAG_,ContingencyTableDistribution[variables_, CrossTab_]]:=Module[{jointDistVariables,newDist,randomVarName,condDist,approxJointDist,margDist,randomVar,randomVars,condVar},
	randomVars=VertexList[DAG];
	approxJointDist=ConstantArray[1,Dimensions@CrossTab];
	jointDistVariables=variables;
	Do[randomVar=randomVars[[i]];
		randomVarName=variables[[randomVar]];
		condVar=variables[[#[[1]]&/@EdgeList[DAG,_\[DirectedEdge]randomVar]]]; (* we condition over the parents of randomVar *)
		margDist=MarginalDistribution[ContingencyTableDistribution[variables, CrossTab],Join[condVar,{randomVarName}]]; (* we will marginalize over these variables *)
		condDist=ConditionalDistribution[margDist,condVar]; 
		newDist=BringToFront[approxJointDist, Length@variables, IndexPosition[condDist[[1]],jointDistVariables], jointDistVariables];
		approxJointDist=N@newDist[[1]] condDist[[2]];
		jointDistVariables=newDist[[2]]
		,{i,Length@randomVars}];
	newDist=BringToFront[approxJointDist, Length@variables, IndexPosition[variables,jointDistVariables], jointDistVariables];
	ContingencyTableDistribution[newDist[[2]], newDist[[1]]]
	]

Kullback-Leibler divergence, (Subscript[D, KL])(P||Q).

KLDivergence[contingencyTableDist_,contingencyTableDistApprox_]:=Module[{dist1,dist2},
	dist1=N@contingencyTableDist[[2]]/. {0. -> 1};
	dist2=N@contingencyTableDistApprox[[2]]/. {0. -> 1};
	Total[dist1 Log[dist1]- dist1 Log[dist2],Infinity]
	]

Main Results

Probability Functions: 4d Contingency Tables

We make use of the fact that the Probability function is completely independent of the various distribution functions (such as ConditionalDistribution). Thus we check them against each other.

SeedRandom[1234]
Crosstabulation=RandomInteger[{1,100},{2,10,2,2}]; (* normalized already *)
CrosstabulationN=ContingencyNormalize[Crosstabulation];
titles={"sex","age","TB","alcoholic"};
distribution=ContingencyTableDistribution[titles,CrosstabulationN];

P ("sex" = 1 && age = 7&& TB ==1 && alcoholic ==2)

CrosstabulationN[[1,7,1,2]]
Probability[x==1∧y==7∧z==1∧a==2,{x,y,z,a}\[Distributed]distribution]
EqualityContingencyTable[distribution,{{"age",7},{"TB",1},{"sex",1},{"alcoholic",2}}][[2]]

P (“sex” = 1 && age = 7 )

Probability[x==1∧y==7,{x,y,z,a}\[Distributed]distribution]
t11=MarginalDistribution[distribution,{"sex","age"}];
Probability[sex==1∧age==7,{sex,age}\[Distributed]t11]

P (“sex” = 1 && age = 7 | TB ==1 && alcoholic ==2)

Probability[x==1∧y==7\[Conditioned]z==1∧a==2,{x,y,z,a}\[Distributed]distribution]
cond=ConditionalDistribution[distribution,{"TB","alcoholic"}];
EqualityContingencyTable[cond,{{"TB",1},{"age",7},{"sex",1},{"alcoholic",2}}][[2]]

Model 1

Look at the model from page 15 of http://www.seas.upenn.edu/~taskar/pubs/gms-srl07.pdf. The PC algorithm starts with the complete graph:

sample=10000;pvalue=0.05;
contingency=RandomModelDistribution1[sample];
(* PC algorithm *)
PCAlgorithm[pvalue,sample,contingency]

Model 2

DAG from: http://www.ke.tu-darmstadt.de/lehre/archiv/ss09/ki/bayesian-networks1.pdf Page 24

sample2=10000;pvalue2=0.05;
contingency2=RandomModelDistribution2[sample2];
(* PC algorithm. *)
pc2=PCAlgorithm[pvalue2,sample2,contingency2]

Test the goodness of the approximation

(* Simple test: All variables equal to 1 *)
JointApproxFullForm[pc2[[4]],contingency2][[2]][[1,1,1,1,1,1,1,1]] (* approximation *)
N@contingency2[[2]][[1,1,1,1,1,1,1,1]] (* exact *)

(* Simple Test: All variables equal to 2 *)
JointApproxFullForm[pc2[[4]],contingency2][[2]][[2,2,2,2,2,2,2,2]](* approximation *)
N@contingency2[[2]][[2,2,2,2,2,2,2,2]](* exact *)

(* KL divergence. The smaller the closer the distributions are. *)

KLDivergence[contingency2,JointApproxFullForm[pc2[[4]],contingency2]]

Extensions to Project

Implement additional tests for conditional independence, such as the Bayesian Dirichlet score (See Section 3.2 http://www.utia.cas.cz/files/mtr/pgm06/41_paper.pdf)

For continuous variables, there exist kernel-based conditional independence tests (http://arxiv.org/ftp/arxiv/papers/1202/1202.3775.pdf)

There exist approaches to increase the robustness of the PC algorithm, for example Markus Kalischa & Peter Bühlmanna,  “Robustification of the PC-Algorithm for Directed Acyclic Graphs”, Journal of Computational and Graphical Statistics, Volume 17, Issue 4, 2008. (ftp://ftp.stat.math.ethz.ch/Manuscripts/buhlmann/pcrobMay08.pdf)

Implement other algorithms for Bayesian Network discovery, such as the work “The max-min hill-climbing Bayesian network structure learning algorithm” by I. Tsamardinos·et al. Machine Learning (2006) 65:31–78 http://goo.gl/t210tx.

Compare the Bayesian Network causal decomposition approach to approximating the joint distribution to other approximation schemes. See for example “Approximation of Discrete Multivariate Probability Distributions: Recursive and Hierarchical Approaches” by D. Fass