% MIT MAS622j/1.126J, Problem Set 3
% Conditional Dependence, Explaining Away
% (Burglary-->Alarm<---Earthquake example)
% 
% (for Problem 2(g))

clear all
N = 3; 
dag = zeros(N,N);
B = 1; E = 2; A = 3; 
dag(B,A) = 1;
dag(E,A) = 1;

false= 1; true= 2;
ns = 2*ones(1,N); % binary nodes
names = {'B', 'E', 'A'};

% Set up the topology of the Bayes Net
bnet = mk_bnet(dag, ns, 'names', names, 'discrete', 1:N);

% Display the graph: download graph_draw from the link on the class webpage
graph_draw(bnet.dag, 'node_labels', names);

% Set up conditional probability tables for each node 
% ### Your code comes here! ###
% Hint: Fill in the row vectors with CPD values
bnet.CPD{B} = tabular_CPD(bnet, B, [0.9  0.1]);
bnet.CPD{E} = tabular_CPD(bnet, E, [0.99 0.01]);
bnet.CPD{A} = tabular_CPD(bnet, A, [1 0 0 0 0 1 1 1]);

% Choose an inference engine
engine = jtree_inf_engine(bnet);

% First, compute P(A)
evidence = cell(1,N);
[engine, ll] = enter_evidence(engine, evidence);

m = marginal_nodes(engine, A);
disp(sprintf('p(A) = %g', m.T(true))); 

% Second, compute P(B|A)
% ### Your code comes here! ###
% Hint: Enter the observation that Alarm = true
evidence = cell(1,N);
evidence{A} = true;
[engine, ll] = enter_evidence(engine, evidence);

m = marginal_nodes(engine, B);
disp(sprintf('p(B|A) = %g', m.T(true))); 

% Third, compute P(B|E,A)
% ### Your code comes here! ###
% Hint: Enter the observation that Earthquake = true and Alarm = true
evidence = cell(1,N);
evidence{E} = true;
evidence{A} = true;
[engine, ll] = enter_evidence(engine, evidence);

m = marginal_nodes(engine, B);
disp(sprintf('p(B|E,A) = %g', m.T(true))); 

% Fourth, compute P(E|A)
% ### Your code comes here! ###
% Hint: Enter the observation that Earthquake = true
evidence = cell(1,N);
evidence{A} = true;
[engine, ll] = enter_evidence(engine, evidence);

m = marginal_nodes(engine, E);
disp(sprintf('p(E|A) = %g', m.T(true))); 

% Fifth, compute P(E|B,A)
% ### Your code comes here! ###
% Hint: Enter the observation that Burglary = true and Alarm = true
evidence = cell(1,N);
evidence{B} = true;
evidence{A} = true;
[engine, ll] = enter_evidence(engine, evidence);

m = marginal_nodes(engine, E);
disp(sprintf('p(E|B,A) = %g', m.T(true))); 
