wolffd@0: % Lawn sprinker example from Russell and Norvig p454
wolffd@0: % See www.cs.berkeley.edu/~murphyk/Bayes/usage.html for details.
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.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: CPT = cell(1,N);
wolffd@0: for i=1:N
wolffd@0:   s=struct(bnet.CPD{i});  % violate object privacy
wolffd@0:   CPT{i}=s.CPT;
wolffd@0: end
wolffd@0: 
wolffd@0: % Generate training data
wolffd@0: nsamples = 50;
wolffd@0: samples = cell(N, nsamples);
wolffd@0: for i=1:nsamples
wolffd@0:   samples(:,i) = sample_bnet(bnet);
wolffd@0: end
wolffd@0: data = cell2num(samples);
wolffd@0: 
wolffd@0: % Make a tabula rasa
wolffd@0: bnet2 = mk_bnet(dag, ns);
wolffd@0: seed = 0;
wolffd@0: rand('state', seed);
wolffd@0: bnet2.CPD{C} = tabular_CPD(bnet2, C, 'clamped', 1, 'CPT', [0.5 0.5], ...
wolffd@0: 			   'prior_type', 'dirichlet', 'dirichlet_weight', 0);
wolffd@0: bnet2.CPD{R} = tabular_CPD(bnet2, R, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
wolffd@0: bnet2.CPD{S} = tabular_CPD(bnet2, S, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
wolffd@0: bnet2.CPD{W} = tabular_CPD(bnet2, W, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
wolffd@0: 
wolffd@0: 
wolffd@0: % Find MLEs from fully observed data
wolffd@0: bnet4 = learn_params(bnet2, samples);
wolffd@0: 
wolffd@0: % Bayesian updating with 0 prior is equivalent to ML estimation
wolffd@0: bnet5 = bayes_update_params(bnet2, samples);
wolffd@0: 
wolffd@0: CPT4 = cell(1,N);
wolffd@0: for i=1:N
wolffd@0:   s=struct(bnet4.CPD{i});  % violate object privacy
wolffd@0:   CPT4{i}=s.CPT;
wolffd@0: end
wolffd@0: 
wolffd@0: CPT5 = cell(1,N);
wolffd@0: for i=1:N
wolffd@0:   s=struct(bnet5.CPD{i});  % violate object privacy
wolffd@0:   CPT5{i}=s.CPT;
wolffd@0:   assert(approxeq(CPT5{i}, CPT4{i}))
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: if 1
wolffd@0: % Find MLEs from partially observed data
wolffd@0: 
wolffd@0: % hide 50% of the nodes
wolffd@0: samplesH = samples;
wolffd@0: hide = rand(N, nsamples) > 0.5;
wolffd@0: [I,J]=find(hide);
wolffd@0: for k=1:length(I)
wolffd@0:   samplesH{I(k), J(k)} = [];
wolffd@0: end
wolffd@0: 
wolffd@0: engine = jtree_inf_engine(bnet2);
wolffd@0: max_iter = 5;
wolffd@0: [bnet6, LL] = learn_params_em(engine, samplesH, max_iter);
wolffd@0: 
wolffd@0: CPT6 = cell(1,N);
wolffd@0: for i=1:N
wolffd@0:   s=struct(bnet6.CPD{i});  % violate object privacy
wolffd@0:   CPT6{i}=s.CPT;
wolffd@0: end
wolffd@0: 
wolffd@0: end