wolffd@0: % Lawn sprinker example from Russell and Norvig p454
wolffd@0: % For a picture, see http://www.cs.berkeley.edu/~murphyk/Bayes/usage.html#basics
wolffd@0: 
wolffd@0: N = 4; 
wolffd@0: dag = zeros(N,N);
wolffd@0: C = 1; S = 2; R = 3; W = 4;
wolffd@0: dag(C,[R S]) = 1;
wolffd@0: dag(R,W) = 1;
wolffd@0: dag(S,W)=1;
wolffd@0: 
wolffd@0: false = 1; true = 2;
wolffd@0: ns = 2*ones(1,N); % binary nodes
wolffd@0: 
wolffd@0: %bnet = mk_bnet(dag, ns);
wolffd@0: bnet = mk_bnet(dag, ns, 'names', {'cloudy','S','R','W'}, 'discrete', 1:4);
wolffd@0: names = bnet.names;
wolffd@0: %C = names{'cloudy'};
wolffd@0: bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
wolffd@0: bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
wolffd@0: bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
wolffd@0: bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
wolffd@0: 
wolffd@0: 
wolffd@0: CPD{C} = reshape([0.5 0.5], 2, 1);
wolffd@0: CPD{R} = reshape([0.8 0.2 0.2 0.8], 2, 2);
wolffd@0: CPD{S} = reshape([0.5 0.9 0.5 0.1], 2, 2);
wolffd@0: CPD{W} = reshape([1 0.1 0.1 0.01 0 0.9 0.9 0.99], 2, 2, 2);
wolffd@0: joint = zeros(2,2,2,2);
wolffd@0: for c=1:2
wolffd@0:   for r=1:2
wolffd@0:     for s=1:2
wolffd@0:       for w=1:2
wolffd@0: 	joint(c,s,r,w) = CPD{C}(c) * CPD{S}(c,s) * CPD{R}(c,r) * ...
wolffd@0: 	    CPD{W}(s,r,w);
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: joint2 = repmat(reshape(CPD{C}, [2 1 1 1]), [1 2 2 2]) .* ...
wolffd@0: 	 repmat(reshape(CPD{S}, [2 2 1 1]), [1 1 2 2]) .* ...
wolffd@0: 	 repmat(reshape(CPD{R}, [2 1 2 1]), [1 2 1 2]) .* ...
wolffd@0: 	 repmat(reshape(CPD{W}, [1 2 2 2]), [2 1 1 1]);
wolffd@0: 
wolffd@0: assert(approxeq(joint, joint2));
wolffd@0: 
wolffd@0: 
wolffd@0: engine = jtree_inf_engine(bnet);
wolffd@0: 
wolffd@0: evidence = cell(1,N);
wolffd@0: evidence{W} = true;
wolffd@0: 
wolffd@0: [engine, ll] = enter_evidence(engine, evidence);
wolffd@0: 
wolffd@0: m = marginal_nodes(engine, S);
wolffd@0: p1 = m.T(true) % P(S=true|W=true) = 0.4298
wolffd@0: lik1 = exp(ll); % P(W=true) = 0.6471
wolffd@0: assert(approxeq(p1, 0.4298));
wolffd@0: assert(approxeq(lik1, 0.6471));
wolffd@0: 
wolffd@0: pSandW = sumv(joint(:,true,:,true), [C R]); % P(S,W) = sum_cr P(CSRW)
wolffd@0: pW = sumv(joint(:,:,:,true), [C S R]);
wolffd@0: pSgivenW = pSandW / pW; % P(S=t|W=t) = P(S=t,W=t)/P(W=t)
wolffd@0: assert(approxeq(pW, lik1))
wolffd@0: assert(approxeq(pSgivenW, p1))
wolffd@0: 
wolffd@0: 
wolffd@0: m = marginal_nodes(engine, R);
wolffd@0: p2 = m.T(true)  % P(R=true|W=true) =  0.7079     
wolffd@0: 
wolffd@0: pRandW = sumv(joint(:,:,true,true), [C S]); % P(R,W) = sum_cr P(CSRW)
wolffd@0: pRgivenW = pRandW / pW; % P(R=t|W=t) = P(R=t,W=t)/P(W=t)
wolffd@0: assert(approxeq(pRgivenW, p2))
wolffd@0: 
wolffd@0: 
wolffd@0: % Add extra evidence that R=true
wolffd@0: evidence{R} = true;
wolffd@0: 
wolffd@0: [engine, ll] = enter_evidence(engine, evidence);
wolffd@0: 
wolffd@0: m = marginal_nodes(engine, S);
wolffd@0: p3 = m.T(true) % P(S=true|W=true,R=true) = 0.1945 
wolffd@0: assert(approxeq(p3, 0.1945))
wolffd@0: 
wolffd@0: 
wolffd@0: pSandRandW = sumv(joint(:,true,true,true), [C]); % P(S,R,W) = sum_c P(cSRW)
wolffd@0: pRandW = sumv(joint(:,:,true,true), [C S]); % P(R,W) = sum_cs P(cSRW)
wolffd@0: pSgivenWR = pSandRandW / pRandW; % P(S=t|W=t,R=t) = P(S=t,R=t,W=t)/P(W=t,R=t)
wolffd@0: assert(approxeq(pSgivenWR, p3))
wolffd@0: 
wolffd@0: % So the sprinkler is less likely to be on if we know that
wolffd@0: % it is raining, since the rain can "explain away" the fact
wolffd@0: % that the grass is wet.
wolffd@0: 
wolffd@0: lik3 = exp(ll); % P(W=true, R=true) = 0.4581
wolffd@0: % So the combined evidence is less likely (of course)
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: % Joint distributions
wolffd@0: 
wolffd@0: evidence = cell(1,N);
wolffd@0: [engine, ll] = enter_evidence(engine, evidence);
wolffd@0: m = marginal_nodes(engine, [S R W]);
wolffd@0: 
wolffd@0: evidence{R} = 2;
wolffd@0: [engine, ll] = enter_evidence(engine, evidence);
wolffd@0: m = marginal_nodes(engine, [S R W]);
wolffd@0: 
wolffd@0: 
wolffd@0: