|
wolffd@0
|
1 % Online Bayesian model selection demo.
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 % We generate data from the model A->B
|
|
wolffd@0
|
4 % and compute the posterior prob of all 3 dags on 2 nodes:
|
|
wolffd@0
|
5 % (1) A B, (2) A <- B , (3) A -> B
|
|
wolffd@0
|
6 % Models 2 and 3 are Markov equivalent, and therefore indistinguishable from
|
|
wolffd@0
|
7 % observational data alone.
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 % We control the dependence of B on A by setting
|
|
wolffd@0
|
10 % P(B|A) = 0.5 - epislon and vary epsilon
|
|
wolffd@0
|
11 % as in Koller & Friedman book p512
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % ground truth
|
|
wolffd@0
|
14 N = 2;
|
|
wolffd@0
|
15 dag = zeros(N);
|
|
wolffd@0
|
16 A = 1; B = 2;
|
|
wolffd@0
|
17 dag(A,B) = 1;
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 ntrials = 100;
|
|
wolffd@0
|
20 ns = 2*ones(1,N);
|
|
wolffd@0
|
21 true_bnet = mk_bnet(dag, ns);
|
|
wolffd@0
|
22 true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]);
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 % hypothesis space
|
|
wolffd@0
|
25 G = mk_all_dags(N);
|
|
wolffd@0
|
26 nhyp = length(G);
|
|
wolffd@0
|
27 hyp_bnet = cell(1, nhyp);
|
|
wolffd@0
|
28 for h=1:nhyp
|
|
wolffd@0
|
29 hyp_bnet{h} = mk_bnet(G{h}, ns);
|
|
wolffd@0
|
30 for i=1:N
|
|
wolffd@0
|
31 % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct
|
|
wolffd@0
|
32 % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable.
|
|
wolffd@0
|
33 % The uniform Dirichlet prior is not score equivalent...
|
|
wolffd@0
|
34 fam = family(G{h}, i);
|
|
wolffd@0
|
35 hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ...
|
|
wolffd@0
|
36 'CPT', 'unif');
|
|
wolffd@0
|
37 end
|
|
wolffd@0
|
38 end
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 clf
|
|
wolffd@0
|
41 seeds = 1:3;
|
|
wolffd@0
|
42 expt = 1;
|
|
wolffd@0
|
43 for seedi=1:length(seeds)
|
|
wolffd@0
|
44 seed = seeds(seedi);
|
|
wolffd@0
|
45 rand('state', seed);
|
|
wolffd@0
|
46 randn('state', seed);
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48 es = [0.05 0.1 0.15 0.2];
|
|
wolffd@0
|
49 for ei=1:length(es)
|
|
wolffd@0
|
50 e = es(ei);
|
|
wolffd@0
|
51 true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [0.5+e 0.5-e; 0.5-e 0.5+e]);
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 prior = normalise(ones(1, nhyp));
|
|
wolffd@0
|
54 hyp_w = zeros(ntrials+1, nhyp);
|
|
wolffd@0
|
55 hyp_w(1,:) = prior(:)';
|
|
wolffd@0
|
56 LL = zeros(1, nhyp);
|
|
wolffd@0
|
57 ll = zeros(1, nhyp);
|
|
wolffd@0
|
58 for t=1:ntrials
|
|
wolffd@0
|
59 ev = cell2num(sample_bnet(true_bnet));
|
|
wolffd@0
|
60 for i=1:nhyp
|
|
wolffd@0
|
61 ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev);
|
|
wolffd@0
|
62 hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev);
|
|
wolffd@0
|
63 end
|
|
wolffd@0
|
64 prior = normalise(prior .* exp(ll));
|
|
wolffd@0
|
65 LL = LL + ll;
|
|
wolffd@0
|
66 hyp_w(t+1,:) = prior;
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 % Plot posterior model probabilities
|
|
wolffd@0
|
70 % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)
|
|
wolffd@0
|
71 % Blue = model 2 (2->1)
|
|
wolffd@0
|
72 % Green = model 3 (1->2, "ground truth")
|
|
wolffd@0
|
73
|
|
wolffd@0
|
74 subplot2(length(seeds), length(es), seedi, ei);
|
|
wolffd@0
|
75 m = size(hyp_w,1);
|
|
wolffd@0
|
76 h=plot(1:m, hyp_w(:,1), 'r-', 1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:');
|
|
wolffd@0
|
77 axis([0 m 0 1])
|
|
wolffd@0
|
78 %title('model posterior vs. time')
|
|
wolffd@0
|
79 title(sprintf('e=%3.2f, seed=%d', e, seed));
|
|
wolffd@0
|
80 drawnow
|
|
wolffd@0
|
81 expt = expt + 1;
|
|
wolffd@0
|
82 end
|
|
wolffd@0
|
83 end
|
