set(gca,'xticklabel',{'With Bronchitis','Without Bronchitis'},'yticklabel', {'Smoking','Lung Cancer'});
title('Conditional probabilities with evidence of abnormal x-ray results')
view(50,35);

%%
% When bronchitis is ruled out (|B = f|), the probability of smoking
% history decreases with respect to the case in which the
% bronchitis is confirmed (|B = t|). The effect is propagated across the
% network and affects the probability of lung cancer in a similar manner.


%% Expanding the Network 
% Among various symptoms related to lung cancer and bronchitis is shortness
% of breath (dyspnea). We want to model the relationship of this condition
% within the considered Bayesian Network. We introduce a node |D| and modify
% the adjacency matrix accordingly.

%=== add node D to the network
D = 5;
CPT{D}(:,:,t) = [.75 .1; .5 .05];
CPT{D}(:,:,f) = 1 - CPT{D}(:,:,t);
values{D} =  [1 2];
adj(end+1,:) = [0 0 0 0];
adj(:,end+1) = [0 1 1 0 0];

%=== draw the updated network 
nodeLabels = {'Smoking', 'Bronchitis', 'Lung Cancer', 'Xrays', 'Dyspnea'};
nodeSymbols = {'S', 'B', 'L', 'X', 'D'};
bg = biograph(adj, nodeLabels, 'arrowsize', 4);
nodeHandles= bg.Nodes;
set(nodeHandles, 'shape', 'ellipse');
view(bg)
   
%%
% With the introduction of node |D|, the network is not singly-connected
% anymore. In fact, there are more than one chain between any two nodes
% (i.e., |S| and |D|). We can check this property by considering the undirected
% graph associated with the network and veryfying that it is not acyclic.

isAcyclic = graphisdag(sparse(adj | adj'))

%%
% In order to use the algorithm for exact inference described above, we
% must transform the new, multiply-connected network into a singly-connected
% network. Several approaches can be used, including clustering of parent
% nodes (in our case |B| and |L|) into a single node as follows. 
% 
% First, we combine the adjacency matrix entries corresponding to the nodes
% |B| and |L| into one entry associated to node |BL|. The node |BL| can take up to
% four values, corresponding to all possible combinations of values of the
% original nodes |B| and |L|. Then, we update the conditional probability
% distribution considering that |B| and |L| are conditionally independent given
% the node |S|, that is P(BL|S) = P(B,L|S) = P(B|S) * P(L|S).

%=== combine B and L
adj(B,:) = adj(B,:) | adj(L,:); 
adj(:,B) = adj(:,B) | adj(:,L);
adj(L,:) = []; adj(:,L) = [];   

%=== update the probability distribution accordingly
b1 = kron(CPT{B}(1,:), CPT{L}(1,:));  
b2 = kron(CPT{B}(2,:), CPT{L}(2,:));
x = [CPT{X}(1,:) CPT{X}(1,:)];
d =  reshape((CPT{D}(:,:,1))', 1, 4);

%=== update the node values
S = 1; BL = 2; X = 3; D = 4; 
nodeNames = {'S', 'BL', 'X', 'D'}; 
tt = 1; tf = 2; ft = 3; ff = 4;
values{BL} = 1:4; 
values(L) = [];

%=== create a clustered Conditional Probability Table
cCPT{S} = CPT{S}; 
cCPT{BL}(t,:) = b1; cCPT{BL}(f,:) = b2;
cCPT{D}(:,t) = d; cCPT{D}(:,f) = 1 - d;
cCPT{X}(:,t) = x; cCPT{X}(:,f) = 1 - x;

%=== create and intiate the net
root = find(sum(adj,1)==0); % root (node with no parent)
[cNodes, cEdges] = bnMsgPassCreate(adj, values, cCPT);
[cNodes, cEdges] = bnMsgPassInitiate(cNodes, cEdges, root);

for i = 1:n
    disp(['P(' nodeNames{i}, '|[]) = ' num2str(cNodes(i).P(1))]);
end

%% Drawing the Expanded Network

%=== draw the network
nodeLabels = {'Smoking', 'Bronchitis or Lung Cancer', 'Abnormal X-rays', 'Dyspnea'};
cbg = biograph(adj, nodeNames, 'arrowsize', 4);
set(cbg.Nodes, 'shape', 'ellipse');
cbgInViewer = view(cbg);

%=== assign relevant info to each node handle
cnodeHandles = cbgInViewer.Nodes;
for i = 1:n
    cnodeHandles(i).UserData.Distribution = [cNodes(i).P];
end

%=== draw customized nodes
set(cnodeHandles, 'shape','circle')
colormap(summer)
cbgInViewer.ShowTextInNodes = 'none';
cbgInViewer.CustomNodeDrawFcn = @(node) customnodedraw(node);
cbgInViewer.Scale = .7
cbgInViewer.dolayout



%% Performing Exact Inference on Clustered Trees
% Suppose a patient complains of dyspnea (|D=t|). We would
% like to evaluate the likelihood that this symptom is related to either
% lung cancer or bronchitis.

[n5, e5, A5, a5] = bnMsgPassUpdate(cNodes, cEdges, [], [], D, t);

%%
% Because node |B| and node |L| are clustered into the node |BL|, we have to
% calculate their individual conditional probabilities by considering the
% appropriate value combinations. The conditional probabilities
% in |BL| correspond to the following |B| and |L| value combinations: |BL = tt| if
% |B = t| and |L = t|; |BL = tf| if |B = t| and |L = f|; |BL = ft| if |B = f| and |L = t|; |BL
% = ff| if |B = f| and |L = f|. Therefore P(B|evidence) is equal to the sum of
% the first two elements of P(BL|evidence), and similarly, P(L|evidence) is
% equal to the sum of the first and third elements in P(BL|evidence).

p(1,:) = n5(S).P; 
p(2,:) = [sum(n5(BL).P([tt,tf])), 1-sum(n5(BL).P([tt,tf]))];     % P(B|evidence)
p(3,:) = [sum(n5(BL).P([tt,ft])), 1-sum(n5(BL).P([tt,ft]))]; % P(L|evidence)
p(4,:) = n5(X).P;
p(5,:) = n5(D).P;

for i = 1:5
    disp(['P(' nodeSymbols{i}, '|D=t) = ' num2str(p(i,1))]);
end
%%
% When dyspnea is present, both the likelihood of bronchitis and lung cancer
% increase. This makes sense, since both illnesses have dyspnea as
% symptom and the patient is  indeed exhibiting this symptom.

%% Explaining Away the Lung Cancer
% As we can see in the graph, the
% dyspnea symptom has dependency both on bronchitis and lung cancer.
% Consider the effect of a bronchitis diagnosis on the likelihood of lung
% cancer.

%=== adjust the CPT to reflect B = 1 before clustering into BL node
B = 2; L = 3; 
CPT{B}(:,1) = [1 1] ; CPT{B}(:,2) = 1 - CPT{B}(:,1);
b1 = kron(CPT{B}(1,:), CPT{L}(1,:));  
b2 = kron(CPT{B}(2,:), CPT{L}(2,:));

%=== create a clustered Conditional Probability Table
BL = 2;  
cCPT{BL}(1,:) = b1; cCPT{BL}(2,:) = b2;

%=== create and intiate the net
root = find(sum(adj,1)==0); % root (node with no parent)
[cNodes, cEdges] = bnMsgPassCreate(adj, values, cCPT);
[cNodes, cEdges] = bnMsgPassInitiate(cNodes, cEdges, root);

%=== instantiate for F = 1
[n7, e7, A7, a7] = bnMsgPassUpdate(cNodes, cEdges, [], [], D, t);
w(1,:) = n7(S).P;
w(2,:) = [sum(n7(BL).P([tt,tf])), 1-sum(n7(BL).P([tt,tf]))];     % P(B|evidence)
w(3,:) = [sum(n7(BL).P([tt,ft])), 1-sum(n7(BL).P([tt,ft]))]; % P(L|evidence)
w(4,:) = n7(X).P;
w(5,:) = n7(D).P;

for i = 1:5
    disp(['P(' nodeSymbols{i}, '|B=t,D=t) = ' num2str(w(i,1))]);
end

%%
% When a patient complains of dyspnea and is diagnosed with bronchitis, the
% conditional probability of lung cancer is lower.
% 
% Consider now the effect of a lung cancer diagnosis on the likellihood of
% bronchitis.

%=== adjust the CPT to reflect L = 1 before clustering into BL node
B = 2; L = 3; 
CPT{B}(:,t) = [.25 .05] ; CPT{B}(:,f) = 1 - CPT{B}(:,t);
CPT{L}(:,t) = [1 1]; CPT{L}(:,f) = 1 - CPT{L}(:,t);

b1 = kron(CPT{B}(t,:), CPT{L}(t,:));  
b2 = kron(CPT{B}(f,:), CPT{L}(f,:));

BL = 2;  
cCPT{BL}(t,:) = b1; cCPT{BL}(f,:) = b2;

%=== create and intiate the net
root = find(sum(adj,1)==0); % root (node with no parent)
[cNodes, cEdges] = bnMsgPassCreate(adj, values, cCPT);
[cNodes, cEdges] = bnMsgPassInitiate(cNodes, cEdges, root);

%=== instantiate for D = 1
[n8, e8, A8, a8] = bnMsgPassUpdate(cNodes, cEdges, [], [], D, t);
v(1,:) = n8(S).P;
v(2,:) = [sum(n8(BL).P([tt,tf])), 1-sum(n8(BL).P([tt,tf]))];     % P(B|evidence)
v(3,:) = [sum(n8(BL).P([tt,ft])), 1-sum(n8(BL).P([tt,ft]))]; % P(L|evidence)
v(4,:) = n8(X).P;
v(5,:) = n8(D).P;

for i = 1:5
    disp(['P(' nodeSymbols{i}, '|L=t,D=t) = ' num2str(v(i,1))]);
end

%%
% If a patient is diagnosed with lung cancer in presence of
% dyspnea, the likelihood of bronchitis decreases significantly.
% This phenomenon is called "explaining away" and refers to the situations
% in which the chances of one cause decrease significantly when the chances
% of the competing cause increase. 

%%
% We can observe the "explaining away" phenomenon in the two situations
% described above by comparing the conditional probabilities of node |L| and
% |B| in the two cases. When the evidence for |B| is high, the likelihood
% of |L| is relatively low, and viceversa, when the evidence for |L| is high, the
% likelihood of |B| is low. 

y = [w(2:3,1)  v(2:3,1)];
figure();
bar(y);
set(gca, 'xticklabel', {'Bronchitis', 'Lung Cancer'});
ylabel('Probability'); title('Explaining away with evidence of dyspnea')
legend('B = t', 'L = t', 'location', 'SouthEastOutside');
colormap(summer)

%% References
% [1] Pearl J., "Probabilistic Reasoning in Intelligent Systems", Morgan
% Kaufmann, San Mateo, California, 1988.
%
% [2] Neapolitan R., "Learning Bayesian Networks", Pearson Prentice Hall,
% Upper Saddle River, New Jersey, 2004.


##### SOURCE END #####-->   </body></html>