Daniel@0: % Do the example in Cooper and Yoo, "Causal discovery from a mixture of experimental and
Daniel@0: % observational data", UAI 99, p120
Daniel@0: 
Daniel@0: N = 2;
Daniel@0: dag = zeros(N);
Daniel@0: A = 1; B = 2;
Daniel@0: dag(A,B) = 1;
Daniel@0: ns = 2*ones(1,N);
Daniel@0: 
Daniel@0: bnet0 = mk_bnet(dag, ns);
Daniel@0: %bnet0.CPD{A} = tabular_CPD(bnet0, A, 'unif', 1);
Daniel@0: bnet0.CPD{A} = tabular_CPD(bnet0, A, 'CPT', 'unif', 'prior_type', 'dirichlet');
Daniel@0: bnet0.CPD{B} = tabular_CPD(bnet0, B, 'CPT', 'unif', 'prior_type', 'dirichlet');
Daniel@0: 
Daniel@0: samples = [2 2;
Daniel@0: 	   2 1; 
Daniel@0: 	   2 2;
Daniel@0: 	   1 1;
Daniel@0: 	   1 2;
Daniel@0: 	   2 2;
Daniel@0: 	   1 1;
Daniel@0: 	   2 2;
Daniel@0: 	   1 2;
Daniel@0: 	   2 1;
Daniel@0: 	   1 1];
Daniel@0: 
Daniel@0: clamped = [0 0;
Daniel@0: 	   0 0;
Daniel@0: 	   0 0;
Daniel@0: 	   0 0;
Daniel@0: 	   0 0;
Daniel@0: 	   1 0;
Daniel@0: 	   1 0;
Daniel@0: 	   0 1;
Daniel@0: 	   0 1;
Daniel@0: 	   0 1;
Daniel@0: 	   0 1];
Daniel@0: 
Daniel@0: nsamples = size(samples, 1);
Daniel@0: 
Daniel@0: % sequential version
Daniel@0: LL = 0;
Daniel@0: bnet = bnet0;
Daniel@0: for l=1:nsamples
Daniel@0:   ev = num2cell(samples(l,:)');
Daniel@0:   manip = find(clamped(l,:)');
Daniel@0:   LL = LL + log_marg_lik_complete(bnet, ev, manip);
Daniel@0:   bnet = bayes_update_params(bnet, ev, manip);
Daniel@0: end
Daniel@0: assert(approxeq(exp(LL), 5.97e-7)) % compare with result from UAI paper
Daniel@0: 
Daniel@0: 
Daniel@0: % batch version
Daniel@0: cases = num2cell(samples');
Daniel@0: LL2 = log_marg_lik_complete(bnet0, cases, clamped');
Daniel@0: bnet2 = bayes_update_params(bnet0, cases, clamped');
Daniel@0: 
Daniel@0: assert(approxeq(LL, LL2))
Daniel@0: 
Daniel@0: for j=1:N
Daniel@0:   s1 = struct(bnet.CPD{j}); % violate object privacy
Daniel@0:   s2 = struct(bnet2.CPD{j});
Daniel@0:   assert(approxeq(s1.CPT, s2.CPT))
Daniel@0: end
Daniel@0: