% BNT_example.m
% CS 8803 PGM  Intro to Probabilistic Graphical Models
% Spring 2005
% Jim Rehg
%
% INSTRUCTIONS: At the end of the file are a number of questions that you must answer for HW 6.
%
% To do this exercise you should use a standard package for inference in
% Bayes nets. There are a number of free programs on the web and you free
% to use whatever you want. However I recommend using the Bayes Net Toolkit
% (BNT) by Kevin Murphy. It is a Matlab package and can be downloaded here:
% http://www.cs.ubc.ca/~murphyk/Software/BNT/bnt.html
% Also useful is the tutorial:
% http://www.cs.ubc.ca/~murphyk/Software/BNT/usage.html
% All of the examples that follow are specialized for BNT, which will be
% the basis for the majority of code-oriented assignments in this class.
%
% Follow Kevin Murphy's instructions to install and test BNT.  You may want to edit your startup.m
% file, which goes in $matlabroot/toolbox/local, where 'matlabroot' is the variable that gives the
% path location where Matlab is installed, which is C:\MATLAB6p5 on my PC. A standard startup.m
% (due to Kevin) follows:
%
% Startup.m - Setup initial paths when Matlab runs
%% Below should be the path to BNT source in your install
% addpath C:\usr\Source\Matlab\BNT\FullBNT\BNT
% add_BNT_to_path
%% Make output use less whitespace
% format compact
%% Go into debugger if there is problem (note that you need to use 'dbquit' to get back out)
% dbstop if error
% dbstop if warning
%% default working directory
% cd C:\Temp
%
% If you prefer emacs over the Matlab editor, you may want to get the 'matlab-mode' package
% from:
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=104&objectType=file
%
%
% This is an example of using BNT to solve a simplified burglar alarm network, a well-known
% example from Judea Pearl's classic text "Probabilistic Inference in Intelligent Systems",
% Morgan Kaufmann 1988.
%
% All variables are boolean and represent:
% B - A burglar has broken into the house of (Sherlock) Holmes
% E - An earthquake (tremor) has occurred near Holme's house
% A - The burglar alarm has gone off. This could happen due to burglary (B)
%     or earthquake (E)
% R - The radio announcer indicates that an earthquake has occurred.
%
% Define graph structure for burglar alarm network (all of the edges are arrows
% directed downwards):
%
%    B   E
%     \ / \
%      A   R
%
N = 4;
dag = zeros(N,N);
B = 1;  % Burglary
E = 2;  % Earthquake
A = 3;  % Alarm
R = 4;  % Radio
% Note that nodes must be numbered in topological order in this method.
dag(E,[A R]) = 1;
dag(B,A) = 1;

% Any DAG can be displayed via draw_graph. Note that you have to keep track of the numbering of
% the nodes to interpret the output.
draw_graph(dag);

% The values of a discrete RV will always be the set {1,...,L} where L is the number of discrete
% values it can take on. Since BNT often uses the value of a variable as the index into its
% probability table, it is important to follow this numbering convention. Note that array indices
% in Matlab start at 1, not 0.  In this example we map true and false values for binary RV's as
% follows:
FALSE = 1;  % Different from the built-in Matlab true/false variables
TRUE = 2;
% This means that if v is a vector containing the probability distribution for binary node x,
% then P(x=0) = v(FALSE) and P(x=1) = v(TRUE). Of course this choice is arbitrary and could be
% reversed, as long as the probability distributions match correctly. Note that we are careful
% to use these defined keywords when entering evidence below.

% Create Bayes Net shell
number_states = 2*ones(1,N); % vector of #values/state, all nodes are binary
bnet = mk_bnet(dag, number_states);
% Above is equivalent to:
%   onodes = [];
%   bnet = mk_bnet(dag, node_sizes, 'discrete', discrete_nodes, 'observed', onodes);
% Specify model parameters manually:
% P(E) - Prior probability of earthquake
bnet.CPD{E} = tabular_CPD(bnet, E, [0.99 0.01]);
% P(B) - Prior probability of burglary
bnet.CPD{B} = tabular_CPD(bnet, B, [0.95 0.05]);
% P(R|E) - Conditional probability of radio given earthquake
bnet.CPD{R} = tabular_CPD(bnet, R, [0.99999 0.35 0.00001 0.65]);
% P(A|B,E)
bnet.CPD{A} = tabular_CPD(bnet, A, [0.9999 0.2 0.3 0.01 0.0001 0.8 0.7 0.99]);

% Any CPT's in a bnet can be printed (and extracted if needed):
fprintf('P(A|B,E):\n');
CPT = showCPT(bnet, A);  % Get the CPT for node A

% Each CPT is a multidimensional array, with parents proceding child in node order and
% following the Matlab convention that the first index toggles fastest.
% So the array for A|B,E (in memory) corresponds to the table
% B E A  P(A|B,E)
% ---------------
% F F F  0.9999
% T F F  0.2
% F T F  0.3
% T T F  0.01
% F F T  0.0001
% T F T  0.8
% F T T  0.7
% T T T  0.99
% The vector of probabilities above can be recognized as the input to tabular_CPD. It is converted
% into a table automatically via (in this case) reshape(..., [2 2 2]). Compare the table above
% to the "raw" output that follows:
fprintf('Layout of CPT for P(A|B,E):\n');
celldisp(CPT);

% Before we can answer queries on our network we have to choose an inference engine
% (we have chosen Variable Elimination, which is the first method we will cover in class)
engine = var_elim_inf_engine(bnet)
% and enter any evidence
ev = cell(1,N); % Empty cell array = no evidence
[engine, loglik] = enter_evidence(engine, ev);
% For the first query there is no evidence and we compute marginal probabilities
marg = marginal_nodes(engine, B);
fprintf('\n\nMarginal distribution P(B):\n');
marg.T
% Note that this agrees with initial specification for B.
% T is a CPT for all of the nodes that remain (in order) after marginalization
% This type of marginalization is a "forward" use of the model, the answer depends only on
% prior model parameters, without any additional evidence

% Now look at B and E together
marg = marginal_nodes(engine, [B,E]);
fprintf('\n\nMarginal joint distribution P(B,E):\n');
dispCPT(marg.T);
sum(marg.T,2)
% Note that the above sum is identical to P(B), as expected.

% Now try something more interesting. Holmes gets a phone call informing him that his alarm has
% gone off (A = 1). The same marginalization procedure can now be used to compute the impact of
% that information on the distribution of B and E.
ev{A} = TRUE;
[engine, loglik] = enter_evidence(engine, ev);
marg = marginal_nodes(engine, [B,E]);
fprintf('\n\nP(B,E|A=1):\n');
dispCPT(marg.T);
% This type of marginalization is a "backward" use of the model, where Bayes rule is applied to
% invert the causal order and take the new information into account.

% =============================== BEGIN HOMEWORK QUESTIONS ========================
%
% PS 2, Problem 5
% Part (a)
% Answer any 2 out of the next 3 questions.
%
% i) Explain why P(B=0, E=0|A=1) is much less than P(B=0, E=0).
% ii) Explain why P(B=1, E=0|A=1) is much greater than P(B=0, E=1|A=1).
% iii) Compare P(B=1) and P(E=1) with P(B=1|A=1) and P(E=1|A=1). Do the differences make sense?
%
% Part (b)
% On his way home, suppose that Holmes turns on his radio and hears a report of a mild
% earthquake in his neighborhood (R=1).
%
% i) Compute P(B|A=1, R=1) - Be sure to show your Matlab code with your solution
% ii) Should he continue to drive home or turn around?
% iii) Compare P(B=1|A=1) and P(E=1|A=1) with P(B=1|A=1,R=1) and P(E=1|A=1,R=1). Do the
%        differences make sense?

ev{R} = TRUE;
[engine, loglik] = enter_evidence(engine, ev);
marg = marginal_nodes(engine, [B,E]);
fprintf('\n\nP(B,E|A=1,R=1):\n');
dispCPT(marg.T);
fprintf('\n\nP(B|A=1,R=1)\n');
sum(marg.T,2)
fprintf('\n\nP(E|A=1,R=1)\n');
sum(marg.T)

% EXTRA CREDIT (Each of the following questions is worth 1/2 of a homework)
%
% i) How would you modify the code to use two other inference methods, Enumerative and Junction
%    Tree, instead of variable elimination? Roughly how do these methods differ?
%
% ii) Change the table for P(A|B,E) in the previous exercise to:
%     bnet.CPD{A} = tabular_CPD(bnet, A, [0.95 0.72 0.8 0.001 0.05 0.28 0.2 0.999]);
%     Examine P(B,E|A=1) and explain why P(B=0,E=0|A=1) is much larger than P(B=1,E=1|A=1).