wolffd@0: % compare BIC and Bayesian score 
wolffd@0: 
wolffd@0: N = 4;
wolffd@0: dag = zeros(N,N);
wolffd@0: %C = 1; S = 2; R = 3; W = 4; % topological order
wolffd@0: C = 4; S = 2; R = 3; W = 1; % arbitrary order
wolffd@0: dag(C,[R S]) = 1;
wolffd@0: dag(R,W) = 1;
wolffd@0: dag(S,W)=1;
wolffd@0: 
wolffd@0: 
wolffd@0: false = 1; true = 2;
wolffd@0: ns = 2*ones(1,N); % binary nodes
wolffd@0: bnet = mk_bnet(dag, ns);
wolffd@0: bnet.CPD{C} = tabular_CPD(bnet, C, 'CPT', [0.5 0.5]);
wolffd@0: bnet.CPD{R} = tabular_CPD(bnet, R, 'CPT', [0.8 0.2 0.2 0.8]);
wolffd@0: bnet.CPD{S} = tabular_CPD(bnet, S, 'CPT', [0.5 0.9 0.5 0.1]);
wolffd@0: bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
wolffd@0: 
wolffd@0: 
wolffd@0: seed = 0;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: ncases = 1000;
wolffd@0: data = cell(N, ncases);
wolffd@0: for m=1:ncases
wolffd@0:   data(:,m) = sample_bnet(bnet);
wolffd@0: end
wolffd@0: 
wolffd@0: priors = [0.1 1 10];
wolffd@0: P = length(priors);
wolffd@0: params = cell(1,P);
wolffd@0: for p=1:P
wolffd@0:   params{p} = cell(1,N);
wolffd@0:   for i=1:N
wolffd@0:     %params{p}{i} = {'prior', priors(p)};
wolffd@0:     params{p}{i} = {'prior_type', 'dirichlet', 'dirichlet_weight', priors(p)};
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: %sz = 1000:1000:10000;
wolffd@0: sz = 10:10:100;
wolffd@0: S = length(sz);
wolffd@0: bic_score = zeros(S, 1);
wolffd@0: bayes_score = zeros(S, P);
wolffd@0: for i=1:S
wolffd@0:   bic_score(i) = score_dags(data(:,1:sz(i)), ns, {dag}, 'scoring_fn', 'bic', 'params', []);
wolffd@0: end
wolffd@0: diff = zeros(S,P);
wolffd@0: for p=1:P
wolffd@0:   for i=1:S
wolffd@0:     bayes_score(i,p) = score_dags(data(:,1:sz(i)), ns, {dag}, 'params', params{p});
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: for p=1:P
wolffd@0:   for i=1:S
wolffd@0:     diff(i,p) = bayes_score(i,p)/ bic_score(i);
wolffd@0:     %diff(i,p) = abs(bayes_score(i,p) - bic_score(i));
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: if 0
wolffd@0: plot(sz, diff(:,1), 'g--*', sz, diff(:,2), 'b-.+', sz, diff(:,3), 'k:s');
wolffd@0: title('Relative BIC error vs. size of data set')
wolffd@0: legend('BDeu 0.1', 'BDeu 1', 'Bdeu 10', 2)
wolffd@0: end
wolffd@0: 
wolffd@0: if 0
wolffd@0: plot(sz, bic_score, 'r-o',  sz, bayes_score(:,1), 'g--*', sz, bayes_score(:,2), 'b-.+', sz, bayes_score(:,3), 'k:s');
wolffd@0: legend('bic', 'BDeu 0.01', 'BDeu 1', 'Bdeu 100')
wolffd@0: ylabel('score')
wolffd@0: title('score vs. size of data set')
wolffd@0: end
wolffd@0: 
wolffd@0: %xlabel('num. data cases')
wolffd@0: 
wolffd@0: %previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb')
wolffd@0: %exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/bic.png', 'format', 'png', 'height', 2, 'color', 'rgb')