Merge pull request #1 from ziejd2/octave_php_separate

Octave php separate

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)