wolffd@0: % Bayesian model selection demo.
wolffd@0: 
wolffd@0: % We generate data from the model A->B
wolffd@0: % and compute the posterior prob of all 3 dags on 2 nodes:
wolffd@0: %  (1) A B,  (2) A <- B , (3) A -> B
wolffd@0: % Models 2 and 3 are Markov equivalent, and therefore indistinguishable from 
wolffd@0: % observational data alone.
wolffd@0: % Using the "difficult" params, the true model only gets a higher posterior after 2000 trials!
wolffd@0: % However, using the noisy NOT gate, the true model wins after 12 trials.
wolffd@0: 
wolffd@0: % ground truth
wolffd@0: N = 2;
wolffd@0: dag = zeros(N);
wolffd@0: A = 1; B = 2; 
wolffd@0: dag(A,B) = 1;
wolffd@0: 
wolffd@0: difficult = 0;
wolffd@0: if difficult
wolffd@0:   ntrials = 2000;
wolffd@0:   ns = 3*ones(1,N);
wolffd@0:   true_bnet = mk_bnet(dag, ns);
wolffd@0:   rand('state', 0);
wolffd@0:   temp = 5;
wolffd@0:   for i=1:N
wolffd@0:     %true_bnet.CPD{i} = tabular_CPD(true_bnet, i, temp);
wolffd@0:     true_bnet.CPD{i} = tabular_CPD(true_bnet, i);
wolffd@0:   end
wolffd@0: else
wolffd@0:   ntrials = 25;
wolffd@0:   ns = 2*ones(1,N);
wolffd@0:   true_bnet = mk_bnet(dag, ns);
wolffd@0:   true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]);
wolffd@0:   pfail = 0.1;
wolffd@0:   psucc = 1-pfail;
wolffd@0:   true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [pfail psucc; psucc pfail]); % NOT gate
wolffd@0: end
wolffd@0: 
wolffd@0: G = mk_all_dags(N);
wolffd@0: nhyp = length(G);
wolffd@0: hyp_bnet = cell(1, nhyp);
wolffd@0: for h=1:nhyp
wolffd@0:   hyp_bnet{h} = mk_bnet(G{h}, ns);
wolffd@0:   for i=1:N
wolffd@0:     % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct
wolffd@0:     % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable.
wolffd@0:     % The uniform Dirichlet prior is not score equivalent...
wolffd@0:     fam = family(G{h}, i);
wolffd@0:     hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ...
wolffd@0: 				    'CPT', 'unif');
wolffd@0:   end
wolffd@0: end
wolffd@0: prior = normalise(ones(1, nhyp));
wolffd@0: 
wolffd@0: % save results before doing sequential updating
wolffd@0: init_hyp_bnet = hyp_bnet; 
wolffd@0: init_prior = prior;
wolffd@0: 
wolffd@0: 
wolffd@0: rand('state', 0);
wolffd@0: hyp_w = zeros(ntrials+1, nhyp);
wolffd@0: hyp_w(1,:) = prior(:)';
wolffd@0: 
wolffd@0: data = zeros(N, ntrials);
wolffd@0: 
wolffd@0: % First we compute the posteriors sequentially
wolffd@0: 
wolffd@0: LL = zeros(1, nhyp);
wolffd@0: ll = zeros(1, nhyp);
wolffd@0: for t=1:ntrials
wolffd@0:   ev = cell2num(sample_bnet(true_bnet));
wolffd@0:   data(:,t) = ev;
wolffd@0:   for i=1:nhyp
wolffd@0:     ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev);
wolffd@0:     hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev);
wolffd@0:   end
wolffd@0:   prior = normalise(prior .* exp(ll));
wolffd@0:   LL = LL + ll;
wolffd@0:   hyp_w(t+1,:) = prior;
wolffd@0: end
wolffd@0: 
wolffd@0: % Plot posterior model probabilities
wolffd@0: % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)
wolffd@0: % Blue = model 2 (2->1)
wolffd@0: % Green = model 3 (1->2, "ground truth")
wolffd@0: 
wolffd@0: if 1
wolffd@0:   figure;
wolffd@0: m = size(hyp_w, 1);
wolffd@0: h=plot(1:m, hyp_w(:,1), 'r-',  1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:');
wolffd@0: axis([0 m   0 1])
wolffd@0: title('model posterior vs. time')
wolffd@0: %previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb')
wolffd@0: %exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/model_select.png',...
wolffd@0: %'format', 'png', 'height', 2, 'color', 'rgb')
wolffd@0: drawnow
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Now check that batch updating gives same result
wolffd@0: hyp_bnet2 = init_hyp_bnet;
wolffd@0: prior2 = init_prior;
wolffd@0: 
wolffd@0: cases = num2cell(data);
wolffd@0: LL2 = zeros(1, nhyp);
wolffd@0: for i=1:nhyp
wolffd@0:   LL2(i) = log_marg_lik_complete(hyp_bnet2{i}, cases);
wolffd@0:   hyp_bnet2{i} = bayes_update_params(hyp_bnet2{i}, cases);
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: assert(approxeq(LL, LL2))
wolffd@0: LL
wolffd@0: 
wolffd@0: for i=1:nhyp
wolffd@0:   for j=1:N
wolffd@0:     s1 = struct(hyp_bnet{i}.CPD{j});
wolffd@0:     s2 = struct(hyp_bnet2{i}.CPD{j});
wolffd@0:     assert(approxeq(s1.CPT, s2.CPT))
wolffd@0:   end
wolffd@0: end
wolffd@0: