BNW using Octave instead of Matlab.

This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning.

I am calling this BNW_1.02. It can be accessed at:
compbio.uthsc.edu/BNW_1.02

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)