diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/examples/static/StructLearn')
11 files changed, 490 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Entries new file mode 100644 index 00000000..0586b211 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Entries @@ -0,0 +1,9 @@ +/bic1.m/1.1.1.1/Wed May 29 15:59:54 2002// +/cooper_yoo.m/1.1.1.1/Wed May 29 15:59:54 2002// +/k2demo1.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mcmc1.m/1.1.1.1/Wed May 29 15:59:54 2002// +/model_select1.m/1.1.1.1/Sat Nov 6 20:55:18 2004// +/model_select2.m/1.1.1.1/Sat Nov 6 21:52:42 2004// +/pc1.m/1.1.1.1/Wed May 29 15:59:54 2002// +/pc2.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Repository new file mode 100644 index 00000000..5b40c54c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/static/StructLearn diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Root b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/bic1.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/bic1.m new file mode 100644 index 00000000..22473564 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/bic1.m @@ -0,0 +1,79 @@ +% compare BIC and Bayesian score + +N = 4; +dag = zeros(N,N); +%C = 1; S = 2; R = 3; W = 4; % topological order +C = 4; S = 2; R = 3; W = 1; % arbitrary order +dag(C,[R S]) = 1; +dag(R,W) = 1; +dag(S,W)=1; + + +false = 1; true = 2; +ns = 2*ones(1,N); % binary nodes +bnet = mk_bnet(dag, ns); +bnet.CPD{C} = tabular_CPD(bnet, C, 'CPT', [0.5 0.5]); +bnet.CPD{R} = tabular_CPD(bnet, R, 'CPT', [0.8 0.2 0.2 0.8]); +bnet.CPD{S} = tabular_CPD(bnet, S, 'CPT', [0.5 0.9 0.5 0.1]); +bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [1 0.1 0.1 0.01 0 0.9 0.9 0.99]); + + +seed = 0; +rand('state', seed); +randn('state', seed); +ncases = 1000; +data = cell(N, ncases); +for m=1:ncases + data(:,m) = sample_bnet(bnet); +end + +priors = [0.1 1 10]; +P = length(priors); +params = cell(1,P); +for p=1:P + params{p} = cell(1,N); + for i=1:N + %params{p}{i} = {'prior', priors(p)}; + params{p}{i} = {'prior_type', 'dirichlet', 'dirichlet_weight', priors(p)}; + end +end + +%sz = 1000:1000:10000; +sz = 10:10:100; +S = length(sz); +bic_score = zeros(S, 1); +bayes_score = zeros(S, P); +for i=1:S + bic_score(i) = score_dags(data(:,1:sz(i)), ns, {dag}, 'scoring_fn', 'bic', 'params', []); +end +diff = zeros(S,P); +for p=1:P + for i=1:S + bayes_score(i,p) = score_dags(data(:,1:sz(i)), ns, {dag}, 'params', params{p}); + end +end + +for p=1:P + for i=1:S + diff(i,p) = bayes_score(i,p)/ bic_score(i); + %diff(i,p) = abs(bayes_score(i,p) - bic_score(i)); + end +end + +if 0 +plot(sz, diff(:,1), 'g--*', sz, diff(:,2), 'b-.+', sz, diff(:,3), 'k:s'); +title('Relative BIC error vs. size of data set') +legend('BDeu 0.1', 'BDeu 1', 'Bdeu 10', 2) +end + +if 0 +plot(sz, bic_score, 'r-o', sz, bayes_score(:,1), 'g--*', sz, bayes_score(:,2), 'b-.+', sz, bayes_score(:,3), 'k:s'); +legend('bic', 'BDeu 0.01', 'BDeu 1', 'Bdeu 100') +ylabel('score') +title('score vs. size of data set') +end + +%xlabel('num. data cases') + +%previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb') +%exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/bic.png', 'format', 'png', 'height', 2, 'color', 'rgb') diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/cooper_yoo.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/cooper_yoo.m new file mode 100644 index 00000000..97ceb44a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/cooper_yoo.m @@ -0,0 +1,65 @@ +% Do the example in Cooper and Yoo, "Causal discovery from a mixture of experimental and +% observational data", UAI 99, p120 + +N = 2; +dag = zeros(N); +A = 1; B = 2; +dag(A,B) = 1; +ns = 2*ones(1,N); + +bnet0 = mk_bnet(dag, ns); +%bnet0.CPD{A} = tabular_CPD(bnet0, A, 'unif', 1); +bnet0.CPD{A} = tabular_CPD(bnet0, A, 'CPT', 'unif', 'prior_type', 'dirichlet'); +bnet0.CPD{B} = tabular_CPD(bnet0, B, 'CPT', 'unif', 'prior_type', 'dirichlet'); + +samples = [2 2; + 2 1; + 2 2; + 1 1; + 1 2; + 2 2; + 1 1; + 2 2; + 1 2; + 2 1; + 1 1]; + +clamped = [0 0; + 0 0; + 0 0; + 0 0; + 0 0; + 1 0; + 1 0; + 0 1; + 0 1; + 0 1; + 0 1]; + +nsamples = size(samples, 1); + +% sequential version +LL = 0; +bnet = bnet0; +for l=1:nsamples + ev = num2cell(samples(l,:)'); + manip = find(clamped(l,:)'); + LL = LL + log_marg_lik_complete(bnet, ev, manip); + bnet = bayes_update_params(bnet, ev, manip); +end +assert(approxeq(exp(LL), 5.97e-7)) % compare with result from UAI paper + + +% batch version +cases = num2cell(samples'); +LL2 = log_marg_lik_complete(bnet0, cases, clamped'); +bnet2 = bayes_update_params(bnet0, cases, clamped'); + +assert(approxeq(LL, LL2)) + +for j=1:N + s1 = struct(bnet.CPD{j}); % violate object privacy + s2 = struct(bnet2.CPD{j}); + assert(approxeq(s1.CPT, s2.CPT)) +end + diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/k2demo1.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/k2demo1.m new file mode 100644 index 00000000..a6288286 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/k2demo1.m @@ -0,0 +1,45 @@ +N = 4; +dag = zeros(N,N); +%C = 1; S = 2; R = 3; W = 4; +C = 4; S = 2; R = 3; W = 1; % arbitrary order +dag(C,[R S]) = 1; +dag(R,W) = 1; +dag(S,W)=1; + +false = 1; true = 2; +ns = 2*ones(1,N); % binary nodes + +bnet = mk_bnet(dag, ns); +bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]); +bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]); +bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]); +bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]); + +seed = 0; +rand('state', seed); +randn('state', seed); +ncases = 100; +data = zeros(N, ncases); +for m=1:ncases + data(:,m) = cell2num(sample_bnet(bnet)); +end + +order = [C S R W]; +max_fan_in = 2; + +%dag2 = learn_struct_K2(data, ns, order, 'max_fan_in', max_fan_in, 'verbose', 'yes'); + +sz = 5:5:50; +for i=1:length(sz) + dag2 = learn_struct_K2(data(:,1:sz(i)), ns, order, 'max_fan_in', max_fan_in); + correct(i) = isequal(dag, dag2); +end +correct + +for i=1:length(sz) + dag3 = learn_struct_K2(data(:,1:sz(i)), ns, order, 'max_fan_in', max_fan_in, 'scoring_fn', 'bic', 'params', []); + correct(i) = isequal(dag, dag3); +end +correct + + diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/mcmc1.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/mcmc1.m new file mode 100644 index 00000000..241d0686 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/mcmc1.m @@ -0,0 +1,35 @@ +% We compare MCMC structure learning with exhaustive enumeration of all dags. + +N = 3; +%N = 4; +dag = mk_rnd_dag(N); +ns = 2*ones(1,N); +bnet = mk_bnet(dag, ns); +for i=1:N + bnet.CPD{i} = tabular_CPD(bnet, i); +end + +ncases = 100; +data = zeros(N, ncases); +for m=1:ncases + data(:,m) = cell2num(sample_bnet(bnet)); +end + +dags = mk_all_dags(N); +score = score_dags(data, ns, dags); +post = normalise(exp(score)); + +[sampled_graphs, accept_ratio] = learn_struct_mcmc(data, ns, 'nsamples', 100, 'burnin', 10); +mcmc_post = mcmc_sample_to_hist(sampled_graphs, dags); + +if 0 + subplot(2,1,1) + bar(post) + subplot(2,1,2) + bar(mcmc_post) + print(gcf, '-djpeg', '/home/cs/murphyk/public_html/Bayes/Figures/mcmc_post.jpg') + + clf + plot(accept_ratio) + print(gcf, '-djpeg', '/home/cs/murphyk/public_html/Bayes/Figures/mcmc_accept.jpg') +end diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select1.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select1.m new file mode 100644 index 00000000..c79131c3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select1.m @@ -0,0 +1,121 @@ +% Bayesian model selection demo. + +% We generate data from the model A->B +% and compute the posterior prob of all 3 dags on 2 nodes: +% (1) A B, (2) A <- B , (3) A -> B +% Models 2 and 3 are Markov equivalent, and therefore indistinguishable from +% observational data alone. +% Using the "difficult" params, the true model only gets a higher posterior after 2000 trials! +% However, using the noisy NOT gate, the true model wins after 12 trials. + +% ground truth +N = 2; +dag = zeros(N); +A = 1; B = 2; +dag(A,B) = 1; + +difficult = 0; +if difficult + ntrials = 2000; + ns = 3*ones(1,N); + true_bnet = mk_bnet(dag, ns); + rand('state', 0); + temp = 5; + for i=1:N + %true_bnet.CPD{i} = tabular_CPD(true_bnet, i, temp); + true_bnet.CPD{i} = tabular_CPD(true_bnet, i); + end +else + ntrials = 25; + ns = 2*ones(1,N); + true_bnet = mk_bnet(dag, ns); + true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]); + pfail = 0.1; + psucc = 1-pfail; + true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [pfail psucc; psucc pfail]); % NOT gate +end + +G = mk_all_dags(N); +nhyp = length(G); +hyp_bnet = cell(1, nhyp); +for h=1:nhyp + hyp_bnet{h} = mk_bnet(G{h}, ns); + for i=1:N + % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct + % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable. + % The uniform Dirichlet prior is not score equivalent... + fam = family(G{h}, i); + hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ... + 'CPT', 'unif'); + end +end +prior = normalise(ones(1, nhyp)); + +% save results before doing sequential updating +init_hyp_bnet = hyp_bnet; +init_prior = prior; + + +rand('state', 0); +hyp_w = zeros(ntrials+1, nhyp); +hyp_w(1,:) = prior(:)'; + +data = zeros(N, ntrials); + +% First we compute the posteriors sequentially + +LL = zeros(1, nhyp); +ll = zeros(1, nhyp); +for t=1:ntrials + ev = cell2num(sample_bnet(true_bnet)); + data(:,t) = ev; + for i=1:nhyp + ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev); + hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev); + end + prior = normalise(prior .* exp(ll)); + LL = LL + ll; + hyp_w(t+1,:) = prior; +end + +% Plot posterior model probabilities +% Red = model 1 (no arcs), blue/green = models 2/3 (1 arc) +% Blue = model 2 (2->1) +% Green = model 3 (1->2, "ground truth") + +if 1 + figure; +m = size(hyp_w, 1); +h=plot(1:m, hyp_w(:,1), 'r-', 1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:'); +axis([0 m 0 1]) +title('model posterior vs. time') +%previewfig(gcf, 'format', 'png', 'height', 2, 'color', 'rgb') +%exportfig(gcf, '/home/cs/murphyk/public_html/Bayes/Figures/model_select.png',... +%'format', 'png', 'height', 2, 'color', 'rgb') +drawnow +end + + +% Now check that batch updating gives same result +hyp_bnet2 = init_hyp_bnet; +prior2 = init_prior; + +cases = num2cell(data); +LL2 = zeros(1, nhyp); +for i=1:nhyp + LL2(i) = log_marg_lik_complete(hyp_bnet2{i}, cases); + hyp_bnet2{i} = bayes_update_params(hyp_bnet2{i}, cases); +end + + +assert(approxeq(LL, LL2)) +LL + +for i=1:nhyp + for j=1:N + s1 = struct(hyp_bnet{i}.CPD{j}); + s2 = struct(hyp_bnet2{i}.CPD{j}); + assert(approxeq(s1.CPT, s2.CPT)) + end +end + diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select2.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select2.m new file mode 100644 index 00000000..d34c75c4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/model_select2.m @@ -0,0 +1,83 @@ +% Online Bayesian model selection demo. + +% We generate data from the model A->B +% and compute the posterior prob of all 3 dags on 2 nodes: +% (1) A B, (2) A <- B , (3) A -> B +% Models 2 and 3 are Markov equivalent, and therefore indistinguishable from +% observational data alone. + +% We control the dependence of B on A by setting +% P(B|A) = 0.5 - epislon and vary epsilon +% as in Koller & Friedman book p512 + +% ground truth +N = 2; +dag = zeros(N); +A = 1; B = 2; +dag(A,B) = 1; + +ntrials = 100; +ns = 2*ones(1,N); +true_bnet = mk_bnet(dag, ns); +true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]); + +% hypothesis space +G = mk_all_dags(N); +nhyp = length(G); +hyp_bnet = cell(1, nhyp); +for h=1:nhyp + hyp_bnet{h} = mk_bnet(G{h}, ns); + for i=1:N + % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct + % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable. + % The uniform Dirichlet prior is not score equivalent... + fam = family(G{h}, i); + hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ... + 'CPT', 'unif'); + end +end + +clf +seeds = 1:3; +expt = 1; +for seedi=1:length(seeds) + seed = seeds(seedi); + rand('state', seed); + randn('state', seed); + + es = [0.05 0.1 0.15 0.2]; + for ei=1:length(es) + e = es(ei); + true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [0.5+e 0.5-e; 0.5-e 0.5+e]); + + prior = normalise(ones(1, nhyp)); + hyp_w = zeros(ntrials+1, nhyp); + hyp_w(1,:) = prior(:)'; + LL = zeros(1, nhyp); + ll = zeros(1, nhyp); + for t=1:ntrials + ev = cell2num(sample_bnet(true_bnet)); + for i=1:nhyp + ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev); + hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev); + end + prior = normalise(prior .* exp(ll)); + LL = LL + ll; + hyp_w(t+1,:) = prior; + end + + % Plot posterior model probabilities + % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc) + % Blue = model 2 (2->1) + % Green = model 3 (1->2, "ground truth") + + subplot2(length(seeds), length(es), seedi, ei); + m = size(hyp_w,1); + h=plot(1:m, hyp_w(:,1), 'r-', 1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:'); + axis([0 m 0 1]) + %title('model posterior vs. time') + title(sprintf('e=%3.2f, seed=%d', e, seed)); + drawnow + expt = expt + 1; + end +end diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc1.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc1.m new file mode 100644 index 00000000..a6fbb9bd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc1.m @@ -0,0 +1,30 @@ +% SGS p118 +% Try learning the structure using an oracle for the cond indep tests + +n = 5; + +A = 1; B = 2; C = 3; D = 4; E = 5; + +G = zeros(n); +G(A,B)=1; +G(B,[C D]) = 1; +G(C,E)=1; +G(D,E)=1; + +k = 2; + +pdag = learn_struct_pdag_pc('dsep', n, k, G) + + + + +if 0 +N = 4; +dag = zeros(N,N); +C = 1; S = 2; R = 3; W = 4; +dag(C,[R S]) = 1; +dag(R,W) = 1; +dag(S,W)=1; + +pdag = learn_struct_pdag_pc('dsep', N, 2, dag) +end diff --git a/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc2.m b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc2.m new file mode 100644 index 00000000..7b7e2066 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/static/StructLearn/pc2.m @@ -0,0 +1,21 @@ +% SGS p141 (female orgasm data set) + +C = eye(7,7); +C(2,1:1) = [-0.132]; +C(3,1:2) = [0.009 -0.136]; +C(4,1:3) = [0.22 -0.166 0.403]; +C(5,1:4) = [-0.008 0.008 0.598 0.282]; +C(6,1:5) = [0.119 -0.076 0.264 0.514 0.176]; +C(7,1:6) = [0.118 -0.137 0.368 0.414 0.336 0.338]; + +n = 7; +for i=1:n + for j=i+1:n + C(i,j)=C(j,i); + end +end + +max_fan_in = 4; +nsamples = 281; +alpha = 0.05; +pdag = learn_struct_pdag_pc('cond_indep_fisher_z', n, max_fan_in, C, nsamples, alpha) |