|
wolffd@0
|
1 % 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 % Using the "difficult" params, the true model only gets a higher posterior after 2000 trials!
|
|
wolffd@0
|
9 % However, using the noisy NOT gate, the true model wins after 12 trials.
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 % ground truth
|
|
wolffd@0
|
12 N = 2;
|
|
wolffd@0
|
13 dag = zeros(N);
|
|
wolffd@0
|
14 A = 1; B = 2;
|
|
wolffd@0
|
15 dag(A,B) = 1;
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 difficult = 0;
|
|
wolffd@0
|
18 if difficult
|
|
wolffd@0
|
19 ntrials = 2000;
|
|
wolffd@0
|
20 ns = 3*ones(1,N);
|
|
wolffd@0
|
21 true_bnet = mk_bnet(dag, ns);
|
|
wolffd@0
|
22 rand('state', 0);
|
|
wolffd@0
|
23 temp = 5;
|
|
wolffd@0
|
24 for i=1:N
|
|
wolffd@0
|
25 %true_bnet.CPD{i} = tabular_CPD(true_bnet, i, temp);
|
|
wolffd@0
|
26 true_bnet.CPD{i} = tabular_CPD(true_bnet, i);
|
|
wolffd@0
|
27 end
|
|
wolffd@0
|
28 else
|
|
wolffd@0
|
29 ntrials = 25;
|
|
wolffd@0
|
30 ns = 2*ones(1,N);
|
|
wolffd@0
|
31 true_bnet = mk_bnet(dag, ns);
|
|
wolffd@0
|
32 true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]);
|
|
wolffd@0
|
33 pfail = 0.1;
|
|
wolffd@0
|
34 psucc = 1-pfail;
|
|
wolffd@0
|
35 true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [pfail psucc; psucc pfail]); % NOT gate
|
|
wolffd@0
|
36 end
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 G = mk_all_dags(N);
|
|
wolffd@0
|
39 nhyp = length(G);
|
|
wolffd@0
|
40 hyp_bnet = cell(1, nhyp);
|
|
wolffd@0
|
41 for h=1:nhyp
|
|
wolffd@0
|
42 hyp_bnet{h} = mk_bnet(G{h}, ns);
|
|
wolffd@0
|
43 for i=1:N
|
|
wolffd@0
|
44 % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct
|
|
wolffd@0
|
45 % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable.
|
|
wolffd@0
|
46 % The uniform Dirichlet prior is not score equivalent...
|
|
wolffd@0
|
47 fam = family(G{h}, i);
|
|
wolffd@0
|
48 hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ...
|
|
wolffd@0
|
49 'CPT', 'unif');
|
|
wolffd@0
|
50 end
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52 prior = normalise(ones(1, nhyp));
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54 % save results before doing sequential updating
|
|
wolffd@0
|
55 init_hyp_bnet = hyp_bnet;
|
|
wolffd@0
|
56 init_prior = prior;
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58
|
|
wolffd@0
|
59 rand('state', 0);
|
|
wolffd@0
|
60 hyp_w = zeros(ntrials+1, nhyp);
|
|
wolffd@0
|
61 hyp_w(1,:) = prior(:)';
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 data = zeros(N, ntrials);
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 % First we compute the posteriors sequentially
|
|
wolffd@0
|
66
|
|
wolffd@0
|
67 LL = zeros(1, nhyp);
|
|
wolffd@0
|
68 ll = zeros(1, nhyp);
|
|
wolffd@0
|
69 for t=1:ntrials
|
|
wolffd@0
|
70 ev = cell2num(sample_bnet(true_bnet));
|
|
wolffd@0
|
71 data(:,t) = ev;
|
|
wolffd@0
|
72 for i=1:nhyp
|
|
wolffd@0
|
73 ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev);
|
|
wolffd@0
|
74 hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev);
|
|
wolffd@0
|
75 end
|
|
wolffd@0
|
76 prior = normalise(prior .* exp(ll));
|
|
wolffd@0
|
77 LL = LL + ll;
|
|
wolffd@0
|
78 hyp_w(t+1,:) = prior;
|
|
wolffd@0
|
79 end
|
|
wolffd@0
|
80
|
|
wolffd@0
|
81 % Plot posterior model probabilities
|
|
wolffd@0
|
82 % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)
|
|
wolffd@0
|
83 % Blue = model 2 (2->1)
|
|
wolffd@0
|
84 % Green = model 3 (1->2, "ground truth")
|
|
wolffd@0
|
85
|
|
wolffd@0
|
86 if 1
|
|
wolffd@0
|
87 figure;
|
|
wolffd@0
|
88 m = size(hyp_w, 1);
|
|
wolffd@0
|
89 h=plot(1:m, hyp_w(:,1), 'r-', 1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:');
|
|
wolffd@0
|
90 axis([0 m 0 1])
|
|
wolffd@0
|
91 title('model posterior vs. time')
|
|
wolffd@0
|
92 %previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb')
|
|
wolffd@0
|
93 %exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/model_select.png',...
|
|
wolffd@0
|
94 %'format', 'png', 'height', 2, 'color', 'rgb')
|
|
wolffd@0
|
95 drawnow
|
|
wolffd@0
|
96 end
|
|
wolffd@0
|
97
|
|
wolffd@0
|
98
|
|
wolffd@0
|
99 % Now check that batch updating gives same result
|
|
wolffd@0
|
100 hyp_bnet2 = init_hyp_bnet;
|
|
wolffd@0
|
101 prior2 = init_prior;
|
|
wolffd@0
|
102
|
|
wolffd@0
|
103 cases = num2cell(data);
|
|
wolffd@0
|
104 LL2 = zeros(1, nhyp);
|
|
wolffd@0
|
105 for i=1:nhyp
|
|
wolffd@0
|
106 LL2(i) = log_marg_lik_complete(hyp_bnet2{i}, cases);
|
|
wolffd@0
|
107 hyp_bnet2{i} = bayes_update_params(hyp_bnet2{i}, cases);
|
|
wolffd@0
|
108 end
|
|
wolffd@0
|
109
|
|
wolffd@0
|
110
|
|
wolffd@0
|
111 assert(approxeq(LL, LL2))
|
|
wolffd@0
|
112 LL
|
|
wolffd@0
|
113
|
|
wolffd@0
|
114 for i=1:nhyp
|
|
wolffd@0
|
115 for j=1:N
|
|
wolffd@0
|
116 s1 = struct(hyp_bnet{i}.CPD{j});
|
|
wolffd@0
|
117 s2 = struct(hyp_bnet2{i}.CPD{j});
|
|
wolffd@0
|
118 assert(approxeq(s1.CPT, s2.CPT))
|
|
wolffd@0
|
119 end
|
|
wolffd@0
|
120 end
|
|
wolffd@0
|
121
|
