|
wolffd@0
|
1 % compare BIC and Bayesian score
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 N = 4;
|
|
wolffd@0
|
4 dag = zeros(N,N);
|
|
wolffd@0
|
5 %C = 1; S = 2; R = 3; W = 4; % topological order
|
|
wolffd@0
|
6 C = 4; S = 2; R = 3; W = 1; % arbitrary order
|
|
wolffd@0
|
7 dag(C,[R S]) = 1;
|
|
wolffd@0
|
8 dag(R,W) = 1;
|
|
wolffd@0
|
9 dag(S,W)=1;
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11
|
|
wolffd@0
|
12 false = 1; true = 2;
|
|
wolffd@0
|
13 ns = 2*ones(1,N); % binary nodes
|
|
wolffd@0
|
14 bnet = mk_bnet(dag, ns);
|
|
wolffd@0
|
15 bnet.CPD{C} = tabular_CPD(bnet, C, 'CPT', [0.5 0.5]);
|
|
wolffd@0
|
16 bnet.CPD{R} = tabular_CPD(bnet, R, 'CPT', [0.8 0.2 0.2 0.8]);
|
|
wolffd@0
|
17 bnet.CPD{S} = tabular_CPD(bnet, S, 'CPT', [0.5 0.9 0.5 0.1]);
|
|
wolffd@0
|
18 bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20
|
|
wolffd@0
|
21 seed = 0;
|
|
wolffd@0
|
22 rand('state', seed);
|
|
wolffd@0
|
23 randn('state', seed);
|
|
wolffd@0
|
24 ncases = 1000;
|
|
wolffd@0
|
25 data = cell(N, ncases);
|
|
wolffd@0
|
26 for m=1:ncases
|
|
wolffd@0
|
27 data(:,m) = sample_bnet(bnet);
|
|
wolffd@0
|
28 end
|
|
wolffd@0
|
29
|
|
wolffd@0
|
30 priors = [0.1 1 10];
|
|
wolffd@0
|
31 P = length(priors);
|
|
wolffd@0
|
32 params = cell(1,P);
|
|
wolffd@0
|
33 for p=1:P
|
|
wolffd@0
|
34 params{p} = cell(1,N);
|
|
wolffd@0
|
35 for i=1:N
|
|
wolffd@0
|
36 %params{p}{i} = {'prior', priors(p)};
|
|
wolffd@0
|
37 params{p}{i} = {'prior_type', 'dirichlet', 'dirichlet_weight', priors(p)};
|
|
wolffd@0
|
38 end
|
|
wolffd@0
|
39 end
|
|
wolffd@0
|
40
|
|
wolffd@0
|
41 %sz = 1000:1000:10000;
|
|
wolffd@0
|
42 sz = 10:10:100;
|
|
wolffd@0
|
43 S = length(sz);
|
|
wolffd@0
|
44 bic_score = zeros(S, 1);
|
|
wolffd@0
|
45 bayes_score = zeros(S, P);
|
|
wolffd@0
|
46 for i=1:S
|
|
wolffd@0
|
47 bic_score(i) = score_dags(data(:,1:sz(i)), ns, {dag}, 'scoring_fn', 'bic', 'params', []);
|
|
wolffd@0
|
48 end
|
|
wolffd@0
|
49 diff = zeros(S,P);
|
|
wolffd@0
|
50 for p=1:P
|
|
wolffd@0
|
51 for i=1:S
|
|
wolffd@0
|
52 bayes_score(i,p) = score_dags(data(:,1:sz(i)), ns, {dag}, 'params', params{p});
|
|
wolffd@0
|
53 end
|
|
wolffd@0
|
54 end
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 for p=1:P
|
|
wolffd@0
|
57 for i=1:S
|
|
wolffd@0
|
58 diff(i,p) = bayes_score(i,p)/ bic_score(i);
|
|
wolffd@0
|
59 %diff(i,p) = abs(bayes_score(i,p) - bic_score(i));
|
|
wolffd@0
|
60 end
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 if 0
|
|
wolffd@0
|
64 plot(sz, diff(:,1), 'g--*', sz, diff(:,2), 'b-.+', sz, diff(:,3), 'k:s');
|
|
wolffd@0
|
65 title('Relative BIC error vs. size of data set')
|
|
wolffd@0
|
66 legend('BDeu 0.1', 'BDeu 1', 'Bdeu 10', 2)
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 if 0
|
|
wolffd@0
|
70 plot(sz, bic_score, 'r-o', sz, bayes_score(:,1), 'g--*', sz, bayes_score(:,2), 'b-.+', sz, bayes_score(:,3), 'k:s');
|
|
wolffd@0
|
71 legend('bic', 'BDeu 0.01', 'BDeu 1', 'Bdeu 100')
|
|
wolffd@0
|
72 ylabel('score')
|
|
wolffd@0
|
73 title('score vs. size of data set')
|
|
wolffd@0
|
74 end
|
|
wolffd@0
|
75
|
|
wolffd@0
|
76 %xlabel('num. data cases')
|
|
wolffd@0
|
77
|
|
wolffd@0
|
78 %previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb')
|
|
wolffd@0
|
79 %exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/bic.png', 'format', 'png', 'height', 2, 'color', 'rgb')
|
