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