Daniel@0: % Online Bayesian model selection demo.
Daniel@0: 
Daniel@0: % We generate data from the model A->B
Daniel@0: % and compute the posterior prob of all 3 dags on 2 nodes:
Daniel@0: %  (1) A B,  (2) A <- B , (3) A -> B
Daniel@0: % Models 2 and 3 are Markov equivalent, and therefore indistinguishable from 
Daniel@0: % observational data alone.
Daniel@0: 
Daniel@0: % We control the dependence of B on A by setting
Daniel@0: % P(B|A) = 0.5 - epislon and vary epsilon
Daniel@0: % as in Koller & Friedman book p512
Daniel@0: 
Daniel@0: % ground truth
Daniel@0: N = 2;
Daniel@0: dag = zeros(N);
Daniel@0: A = 1; B = 2; 
Daniel@0: dag(A,B) = 1;
Daniel@0: 
Daniel@0: ntrials = 100;
Daniel@0: ns = 2*ones(1,N);
Daniel@0: true_bnet = mk_bnet(dag, ns);
Daniel@0: true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]);
Daniel@0: 
Daniel@0: % hypothesis space
Daniel@0: G = mk_all_dags(N);
Daniel@0: nhyp = length(G);
Daniel@0: hyp_bnet = cell(1, nhyp);
Daniel@0: for h=1:nhyp
Daniel@0:   hyp_bnet{h} = mk_bnet(G{h}, ns);
Daniel@0:   for i=1:N
Daniel@0:     % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct
Daniel@0:     % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable.
Daniel@0:     % The uniform Dirichlet prior is not score equivalent...
Daniel@0:     fam = family(G{h}, i);
Daniel@0:     hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ...
Daniel@0: 				    'CPT', 'unif');
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: clf
Daniel@0: seeds = 1:3;
Daniel@0: expt = 1;
Daniel@0: for seedi=1:length(seeds)
Daniel@0:   seed = seeds(seedi);
Daniel@0:   rand('state', seed);
Daniel@0:   randn('state', seed);
Daniel@0:     
Daniel@0:   es = [0.05 0.1 0.15 0.2];
Daniel@0:   for ei=1:length(es)
Daniel@0:     e = es(ei);
Daniel@0:     true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [0.5+e 0.5-e; 0.5-e 0.5+e]);
Daniel@0: 
Daniel@0:     prior = normalise(ones(1, nhyp));
Daniel@0:     hyp_w = zeros(ntrials+1, nhyp);
Daniel@0:     hyp_w(1,:) = prior(:)';
Daniel@0:     LL = zeros(1, nhyp);
Daniel@0:     ll = zeros(1, nhyp);
Daniel@0:     for t=1:ntrials
Daniel@0:       ev = cell2num(sample_bnet(true_bnet));
Daniel@0:       for i=1:nhyp
Daniel@0: 	ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev);
Daniel@0: 	hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev);
Daniel@0:       end
Daniel@0:       prior = normalise(prior .* exp(ll));
Daniel@0:       LL = LL + ll;
Daniel@0:       hyp_w(t+1,:) = prior;
Daniel@0:     end
Daniel@0: 
Daniel@0:     % Plot posterior model probabilities
Daniel@0:     % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)
Daniel@0:     % Blue = model 2 (2->1)
Daniel@0:     % Green = model 3 (1->2, "ground truth")
Daniel@0:     
Daniel@0:     subplot2(length(seeds), length(es), seedi, ei);
Daniel@0:     m = size(hyp_w,1);
Daniel@0:     h=plot(1:m, hyp_w(:,1), 'r-',  1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:');
Daniel@0:     axis([0 m   0 1])
Daniel@0:     %title('model posterior vs. time')
Daniel@0:     title(sprintf('e=%3.2f, seed=%d', e, seed));
Daniel@0:     drawnow
Daniel@0:     expt = expt + 1;
Daniel@0:   end
Daniel@0: end