
% Exercise 2-4.

% conditional probability tables
% note that Matlab indices start from 1 so 0,1 need to be changed to 1,2
PA = [0.5; 0.5];
PB = [0.8 0.2];
PC_given_AB = reshape([0.8 0.6 0.7 0.3 0.2 0.4 0.3 0.7],[2,2,2]);
PD_given_C = reshape([0.9 0.2 0.1 0.8],[1,1,2,2]);

% Use sampling to estimate P(B=1|D=1) 
nsamples = 1000;
data = zeros(nsamples,4);
for i=1:nsamples,
  data(i,1) = rand < PA(2);
  data(i,2) = rand < PB(1,2); 
  data(i,3) = rand < PC_given_AB(data(i,1)+1,data(i,2)+1,2);
  data(i,4) = rand < PD_given_C(1,1,data(i,3)+1,2); 
end;
count_D1 = sum(data(:,4));
count_B1D1 = sum(and(data(:,2)==1,data(:,4)==1));
estimate = count_B1D1/count_D1;

% Find the accurate P(B=1|D=1)
% joint probability of P(A), P(A,B), P(A,B,C), P(A,B,C,D)
P1 = PA;
P2 = repmat(PA,1,2) .* repmat(PB,2,1);
P3 = repmat(P2,[1,1,2]).*PC_given_AB;
P4 = repmat(P3,[1,1,1,2]).*repmat(PD_given_C,[2,2,1,1]);
accurate = sum(sum(P4(:,2,:,2))) / sum(sum(sum(P4(:,:,:,2))));

fprintf('P(B=1|D=1) is estimated as %6.4f and the accurate value is %6.4f', ...
        estimate,accurate);

% Plot a figure of the evolution of the estimate during sampling
x = zeros(nsamples,1);
conf = zeros(nsamples,1);
for i=1:nsamples,
  x(i) = sum(and(data(1:i,2)==1,data(1:i,4)==1))) / ...
         sum(data(1:i,4));
  conf(i) = sum(data(1:i,4))^-0.5;
end;
plot([x, x-conf, x+conf, repmat(accurate,nsamples,1)]);
axis([0 1000 0 1]);


