From 8070dc963753142bb86c4ed698d91fd623ed28e7 Mon Sep 17 00:00:00 2001 From: ziejd2 Date: Thu, 28 Sep 2017 15:04:40 -0500 Subject: 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 --- sourcecodes/bnt-master/BNT/learning/CVS/Entries | 23 ++ sourcecodes/bnt-master/BNT/learning/CVS/Repository | 1 + sourcecodes/bnt-master/BNT/learning/CVS/Root | 1 + sourcecodes/bnt-master/BNT/learning/CovMat.m | 39 +++ .../bnt-master/BNT/learning/bayes_update_params.m | 38 +++ .../bnt-master/BNT/learning/bic_score_family.m | 23 ++ .../BNT/learning/compute_cooling_schedule.m | 13 + .../BNT/learning/dirichlet_score_family.m | 42 +++ .../BNT/learning/kpm_learn_struct_mcmc.m | 265 +++++++++++++++++++ sourcecodes/bnt-master/BNT/learning/learn_params.m | 26 ++ .../bnt-master/BNT/learning/learn_params_dbn.m | 36 +++ .../bnt-master/BNT/learning/learn_params_dbn_em.m | 179 +++++++++++++ .../bnt-master/BNT/learning/learn_params_em.m | 88 +++++++ .../bnt-master/BNT/learning/learn_struct_K2.m | 104 ++++++++ .../BNT/learning/learn_struct_dbn_reveal.m | 101 ++++++++ .../bnt-master/BNT/learning/learn_struct_mcmc.m | 284 +++++++++++++++++++++ .../BNT/learning/learn_struct_pdag_ic_star.m | 155 +++++++++++ .../bnt-master/BNT/learning/learn_struct_pdag_pc.m | 130 ++++++++++ .../BNT/learning/learn_struct_pdag_pc_constrain.m | 177 +++++++++++++ .../bnt-master/BNT/learning/mcmc_sample_to_hist.m | 36 +++ sourcecodes/bnt-master/BNT/learning/mk_schedule.m | 10 + .../bnt-master/BNT/learning/mk_tetrad_data_file.m | 25 ++ .../bnt-master/BNT/learning/score_dags_old.m | 62 +++++ .../bnt-master/BNT/learning/score_dags_wocache.m | 67 +++++ .../bnt-master/BNT/learning/score_family_wocache.m | 57 +++++ 25 files changed, 1982 insertions(+) create mode 100644 sourcecodes/bnt-master/BNT/learning/CVS/Entries create mode 100644 sourcecodes/bnt-master/BNT/learning/CVS/Repository create mode 100644 sourcecodes/bnt-master/BNT/learning/CVS/Root create mode 100644 sourcecodes/bnt-master/BNT/learning/CovMat.m create mode 100644 sourcecodes/bnt-master/BNT/learning/bayes_update_params.m create mode 100644 sourcecodes/bnt-master/BNT/learning/bic_score_family.m create mode 100644 sourcecodes/bnt-master/BNT/learning/compute_cooling_schedule.m create mode 100644 sourcecodes/bnt-master/BNT/learning/dirichlet_score_family.m create mode 100644 sourcecodes/bnt-master/BNT/learning/kpm_learn_struct_mcmc.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_params.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_params_dbn.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_params_dbn_em.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_params_em.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_K2.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_dbn_reveal.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_mcmc.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_ic_star.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc.m create mode 100644 sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc_constrain.m create mode 100644 sourcecodes/bnt-master/BNT/learning/mcmc_sample_to_hist.m create mode 100644 sourcecodes/bnt-master/BNT/learning/mk_schedule.m create mode 100644 sourcecodes/bnt-master/BNT/learning/mk_tetrad_data_file.m create mode 100644 sourcecodes/bnt-master/BNT/learning/score_dags_old.m create mode 100644 sourcecodes/bnt-master/BNT/learning/score_dags_wocache.m create mode 100644 sourcecodes/bnt-master/BNT/learning/score_family_wocache.m (limited to 'sourcecodes/bnt-master/BNT/learning') diff --git a/sourcecodes/bnt-master/BNT/learning/CVS/Entries b/sourcecodes/bnt-master/BNT/learning/CVS/Entries new file mode 100644 index 00000000..02e90b7f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/CVS/Entries @@ -0,0 +1,23 @@ +/CovMat.m/1.1.1.1/Sun Jul 28 23:09:42 2002// +/bayes_update_params.m/1.1.1.1/Wed May 29 15:59:56 2002// +/bic_score_family.m/1.1.1.1/Wed May 29 15:59:56 2002// +/compute_cooling_schedule.m/1.1.1.1/Wed May 29 15:59:56 2002// +/dirichlet_score_family.m/1.1.1.1/Wed May 29 15:59:56 2002// +/kpm_learn_struct_mcmc.m/1.1.1.1/Wed May 29 15:59:56 2002// +/learn_params.m/1.1.1.1/Thu Jun 10 01:28:08 2004// +/learn_params_dbn.m/1.1.1.1/Sun Feb 2 00:23:38 2003// +/learn_params_dbn_em.m/1.1.1.1/Mon Aug 18 21:50:34 2003// +/learn_params_em.m/1.1.1.1/Wed May 29 15:59:56 2002// +/learn_struct_K2.m/1.1.1.1/Thu Sep 26 08:39:16 2002// +/learn_struct_dbn_reveal.m/1.1.1.1/Wed May 29 15:59:56 2002// +/learn_struct_mcmc.m/1.1.1.1/Fri Sep 5 14:06:40 2003// +/learn_struct_pdag_ic_star.m/1.1.1.1/Wed May 29 15:59:56 2002// +/learn_struct_pdag_pc.m/1.2/Mon Feb 12 01:15:13 2007// +/learn_struct_pdag_pc_constrain.m/1.1.1.1/Sun Jul 28 23:09:38 2002// +/mcmc_sample_to_hist.m/1.1.1.1/Wed May 29 15:59:56 2002// +/mk_schedule.m/1.1.1.1/Wed May 29 15:59:56 2002// +/mk_tetrad_data_file.m/1.1.1.1/Wed May 29 15:59:56 2002// +/score_dags.m/1.1.1.1/Mon Apr 4 01:39:46 2005// +/score_dags_old.m/1.1.1.1/Wed May 29 15:59:56 2002// +/score_family.m/1.1.1.1/Thu Jun 10 01:33:14 2004// +D diff --git a/sourcecodes/bnt-master/BNT/learning/CVS/Repository b/sourcecodes/bnt-master/BNT/learning/CVS/Repository new file mode 100644 index 00000000..ed7d464f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/learning diff --git a/sourcecodes/bnt-master/BNT/learning/CVS/Root b/sourcecodes/bnt-master/BNT/learning/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/learning/CovMat.m b/sourcecodes/bnt-master/BNT/learning/CovMat.m new file mode 100644 index 00000000..9b788172 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/CovMat.m @@ -0,0 +1,39 @@ +function [CovMatrix, obs, varfields] = CovMat(filename,row_cols) +%[CovMatrix, obs, varfields] = CovMat(filename,row_cols) +%% generates a Covariance Matrix from a file of data consisting of N columns of M data rows +%% filename string name (with path and extension) of file to open +%% row_cols Number_of_converstions_per_row (turns into [3 inf]) +%% Return +%% CovMatrix Covariance matrix +%% obs Number of observations read in +%% varfields Labels of the variables see filename structure below +%% +%% Filename structure: +%% Comma separated, starting with the variable labels, then the data in rows. +%% filename test.txt consists of: +%% +%% Earthquake,Burglar,Radio,Alarm,Call +%% 1,2,3,4,5 +%% 11,22,33,44,55 +%% . . . +%% +%% Example call: +%% [cvmat numdat lables] = CovMat('test.txt',5); +%% +%% Returns Covariance matrix, number of date rows and variable field names +%% Gary R. Bradski 7/2002 + +fmtstr = '%f'; +for i = 2:row_cols + fmtstr = strcat(fmtstr,',%f'); +end + +%% load data +fidCov = fopen(filename,'r'); + +varfields = fgetl(fidCov); +Corx = fscanf(fidCov,fmtstr,[row_cols inf]); +Corx= Corx'; +[obs bla] = size(Corx); +CovMatrix = cov(Corx); +fclose(fidCov); \ No newline at end of file diff --git a/sourcecodes/bnt-master/BNT/learning/bayes_update_params.m b/sourcecodes/bnt-master/BNT/learning/bayes_update_params.m new file mode 100644 index 00000000..4a0a28f4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/bayes_update_params.m @@ -0,0 +1,38 @@ +function bnet = bayes_update_params(bnet, cases, clamped) +% BAYES_UPDATE_PARAMS Bayesian parameter updating given completely observed data +% bnet = bayes_update_params(bnet, cases, clamped) +% +% If there is a missing data, you must use EM. +% cases(i,m) is the value assigned to node i in case m (this can also be a cell array). +% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros). +% Clamped nodes are not updated. +% If there is a single case, clamped is a list of the clamped nodes, not a bit vector. + + +%if iscell(cases), usecell = 1; else usecell = 0; end + +n = length(bnet.dag); +ncases = size(cases, 2); +if n ~= size(cases, 1) + error('data must be of size nnodes * ncases'); +end + +if ncases == 1 % clamped is a list of nodes + if nargin < 3, clamped = []; end + clamp_set = clamped; + clamped = zeros(n,1); + clamped(clamp_set) = 1; +else % each row of clamped is a bit vector + if nargin < 3, clamped = zeros(n,ncases); end +end + +for i=1:n + e = bnet.equiv_class(i); + if adjustable_CPD(bnet.CPD{e}) + u = find(clamped(i,:)==0); + ps = parents(bnet.dag, i); + bnet.CPD{e} = bayes_update_params(bnet.CPD{e}, cases(i,u), cases(ps,u)); + end +end + + diff --git a/sourcecodes/bnt-master/BNT/learning/bic_score_family.m b/sourcecodes/bnt-master/BNT/learning/bic_score_family.m new file mode 100644 index 00000000..f3725ecb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/bic_score_family.m @@ -0,0 +1,23 @@ +function [S, LL] = bic_score(counts, CPT, ncases) +% BIC_SCORE Bayesian Information Criterion score for a single family +% [S, LL] = bic_score(counts, CPT, ncases) +% +% S is a large sample approximation to the log marginal likelihood, +% which can be computed using dirichlet_score. +% +% S = \log [ prod_j _prod_k theta_ijk ^ N_ijk ] - 0.5*d*log(ncases) +% where counts encode N_ijk, theta_ijk is the MLE comptued from counts, +% and d is the num of free parameters. + +%CPT = mk_stochastic(counts); +tiny = exp(-700); +LL = sum(log(CPT(:) + tiny) .* counts(:)); +% CPT(i) = 0 iff counts(i) = 0 so it is okay to add tiny + +ns = mysize(counts); +ns_ps = ns(1:end-1); +ns_self = ns(end); +nparams = prod([ns_ps (ns_self-1)]); +% sum-to-1 constraint reduces the effective num. vals of the node by 1 + +S = LL - 0.5*nparams*log(ncases); diff --git a/sourcecodes/bnt-master/BNT/learning/compute_cooling_schedule.m b/sourcecodes/bnt-master/BNT/learning/compute_cooling_schedule.m new file mode 100644 index 00000000..397c2f8c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/compute_cooling_schedule.m @@ -0,0 +1,13 @@ +function temp_schedule = compute_cooling_schedule(init_temp, final_temp, anneal_rate) + +temp_schedule = []; +i = 1; +temp_schedule(i)=init_temp; +while temp_schedule(i) > final_temp + i = i + 1; + temp_schedule(i)=temp_schedule(i-1)*anneal_rate; +end + + + + diff --git a/sourcecodes/bnt-master/BNT/learning/dirichlet_score_family.m b/sourcecodes/bnt-master/BNT/learning/dirichlet_score_family.m new file mode 100644 index 00000000..abfb046b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/dirichlet_score_family.m @@ -0,0 +1,42 @@ +function LL = dirichlet_score_family(counts, prior) +% DIRICHLET_SCORE Compute the log marginal likelihood of a single family +% LL = dirichlet_score(counts, prior) +% +% counts(a, b, ..., z) is the number of times parent 1 = a, parent 2 = b, ..., child = z +% prior is an optional multidimensional array of the same shape as counts. +% It defaults to a uniform prior. +% +% We marginalize out the parameters: +% LL = log \int \prod_m P(x(i,m) | x(Pa_i,m), theta_i) P(theta_i) d(theta_i) + + +% LL = log[ prod_j gamma(alpha_ij)/gamma(alpha_ij + N_ij) * +% prod_k gamma(alpha_ijk + N_ijk)/gamma(alpha_ijk) ] +% Call the prod_k term U and the prod_j term V. +% We reshape all quantities into (j,k) matrices +% This formula was first derived by Cooper and Herskovits, 1992. +% See also "Learning Bayesian Networks", Heckerman, Geiger and Chickering, MLJ 95. + +ns = mysize(counts); +ns_ps = ns(1:end-1); +ns_self = ns(end); + +if nargin < 2, prior = normalise(myones(ns)); end + + +if 1 + prior = reshape(prior(:), [prod(ns_ps) ns_self]); + counts = reshape(counts, [prod(ns_ps) ns_self]); + %U = prod(gamma(prior + counts) ./ gamma(prior), 2); % mult over k + LU = sum(gammaln(prior + counts) - gammaln(prior), 2); + alpha_ij = sum(prior, 2); % sum over k + N_ij = sum(counts, 2); + %V = gamma(alpha_ij) ./ gamma(alpha_ij + N_ij); + LV = gammaln(alpha_ij) - gammaln(alpha_ij + N_ij); + %L = prod(U .* V); + LL = sum(LU + LV); +else + CPT = mk_stochastic(prior + counts); + LL = sum(log(CPT(:) .* counts(:))); +end + diff --git a/sourcecodes/bnt-master/BNT/learning/kpm_learn_struct_mcmc.m b/sourcecodes/bnt-master/BNT/learning/kpm_learn_struct_mcmc.m new file mode 100644 index 00000000..b73910d2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/kpm_learn_struct_mcmc.m @@ -0,0 +1,265 @@ +function [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, varargin) +% LEARN_STRUCT_MCMC Monte Carla Markov Chain search over DAGs assuming fully observed data +% [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, ...) +% +% data(i,m) is the value of node i in case m. +% ns(i) is the number of discrete values node i can take on. +% +% sampled_graphs{m} is the m'th sampled graph. +% accept_ratio(t) = acceptance ratio at iteration t +% num_edges(t) = number of edges in model at iteration t +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% nsamples - number of samples to draw from the chain after burn-in [ 100*N ] +% burnin - number of steps to take before drawing samples [ 5*N ] +% init_dag - starting point for the search [ zeros(N,N) ] +% +% e.g., samples = learn_struct_mcmc(data, ns, 'nsamples', 1000); +% +% This interface is not backwards compatible with BNT2, +% but is designed to be compatible with the other learn_struct_xxx routines. +% +% Note: We currently assume a uniform structural prior. + +[n ncases] = size(data); + + +% set default params +type = cell(1,n); +params = cell(1,n); +for i=1:n + type{i} = 'tabular'; + %params{i} = { 'prior', 1 }; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; +clamped = zeros(n, ncases); +nsamples = 100*n; +burnin = 5*n; +dag = zeros(n); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'nsamples', nsamples = args{i+1}; + case 'burnin', burnin = args{i+1}; + case 'init_dag', dag = args{i+1}; + case 'scoring_fn', scoring_fn = args{i+1}; + case 'type', type = args{i+1}; + case 'discrete', discrete = args{i+1}; + case 'clamped', clamped = args{i+1}; + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + end +end + +% We implement the fast acyclicity check described by P. Giudici and R. Castelo, +% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001. +use_giudici = 1; +if use_giudici + [nbrs, ops, nodes] = mk_nbrs_of_digraph(dag); + A = init_ancestor_matrix(dag); +else + [nbrs, ops, nodes] = mk_nbrs_of_dag(dag); + A = []; +end + +num_accepts = 1; +num_rejects = 1; +T = burnin + nsamples; +accept_ratio = zeros(1, T); +num_edges = zeros(1, T); +sampled_graphs = cell(1, nsamples); +%sampled_bitv = zeros(nsamples, n^2); + +for t=1:T + [dag, nbrs, ops, nodes, A, accept] = take_step(dag, nbrs, ops, nodes, ns, data, clamped, A, ... + scoring_fn, discrete, type, params); + num_edges(t) = sum(dag(:)); + num_accepts = num_accepts + accept; + num_rejects = num_rejects + (1-accept); + accept_ratio(t) = num_accepts/num_rejects; + if t > burnin + sampled_graphs{t-burnin} = dag; + %sampled_bitv(t-burnin, :) = dag(:)'; + end +end + + +%%%%%%%%% + + +function [new_dag, new_nbrs, new_ops, new_nodes, A, accept] = ... + take_step(dag, nbrs, ops, nodes, ns, data, clamped, A, ... + scoring_fn, discrete, type, params) + + +use_giudici = ~isempty(A); +if use_giudici + [new_dag, op, i, j] = pick_digraph_nbr(dag, nbrs, ops, nodes, A); + %assert(acyclic(new_dag)); + [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_digraph(new_dag); +else + d = sample_discrete(normalise(ones(1, length(nbrs)))); + new_dag = nbrs{d}; + op = ops{d}; + i = nodes(d, 1); j = nodes(d, 2); + [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_dag(new_dag); +end + +bf = bayes_factor(dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params); + +%R = bf * (new_prior / prior) * (length(nbrs) / length(new_nbrs)); +R = bf * (length(nbrs) / length(new_nbrs)); +u = rand(1,1); +if u > min(1,R) % reject the move + accept = 0; + new_dag = dag; + new_nbrs = nbrs; + new_ops = ops; + new_nodes = nodes; +else + accept = 1; + if use_giudici + A = update_ancestor_matrix(A, op, i, j, new_dag); + end +end + + +%%%%%%%%% + +function bfactor = bayes_factor(old_dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params) + +u = find(clamped(j,:)==0); +LLnew = score_family(j, parents(new_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); +LLold = score_family(j, parents(old_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); +bf1 = exp(LLnew - LLold); + +if strcmp(op, 'rev') % must also multiply in the changes to i's family + u = find(clamped(i,:)==0); + LLnew = score_family(i, parents(new_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i}); + LLold = score_family(i, parents(old_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i}); + bf2 = exp(LLnew - LLold); +else + bf2 = 1; +end +bfactor = bf1 * bf2; + + +%%%%%%%% Giudici stuff follows %%%%%%%%%% + + +function [new_dag, op, i, j] = pick_digraph_nbr(dag, digraph_nbrs, ops, nodes, A) + +legal = 0; +while ~legal + d = sample_discrete(normalise(ones(1, length(digraph_nbrs)))); + i = nodes(d, 1); j = nodes(d, 2); + switch ops{d} + case 'add', + if A(i,j)==0 + legal = 1; + end + case 'del', + legal = 1; + case 'rev', + ps = mysetdiff(parents(dag, j), i); + % if any(A(ps,i)) then there is a path i -> parent of j -> j + % so reversing i->j would create a cycle + legal = ~any(A(ps, i)); + end +end +%new_dag = digraph_nbrs{d}; +new_dag = digraph_nbrs(:,:,d); +op = ops{d}; +i = nodes(d, 1); j = nodes(d, 2); + + +%%%%%%%%%%%%%% + + +function A = update_ancestor_matrix(A, op, i, j, dag) + +switch op + case 'add', + A = do_addition(A, op, i, j, dag); + case 'del', + A = do_removal(A, op, i, j, dag); + case 'rev', + A = do_removal(A, op, i, j, dag); + A = do_addition(A, op, j, i, dag); +end + + +%%%%%%%%%%%% + +function A = do_addition(A, op, i, j, dag) + +A(j,i) = 1; % i is an ancestor of j +anci = find(A(i,:)); +if ~isempty(anci) + A(j,anci) = 1; % all of i's ancestors are added to Anc(j) +end +ancj = find(A(j,:)); +descj = find(A(:,j)); +if ~isempty(ancj) + for k=descj(:)' + A(k,ancj) = 1; % all of j's ancestors are added to each descendant of j + end +end + +%%%%%%%%%%% + +function A = do_removal(A, op, i, j, dag) + +% find all the descendants of j, and put them in topological order +%descj = find(A(:,j)); +R = reachability_graph(dag); +descj = find(R(j,:)); +order = topological_sort(dag); +descj_topnum = order(descj); +[junk, perm] = sort(descj_topnum); +descj = descj(perm); +% Update j and all its descendants +A = update_row(A, j, dag); +for k = descj(:)' + A = update_row(A, k, dag); +end + +%%%%%%%%% + +function A = update_row(A, j, dag) + +% We compute row j of A +A(j, :) = 0; +ps = parents(dag, j); +if ~isempty(ps) + A(j, ps) = 1; +end +for k=ps(:)' + anck = find(A(k,:)); + if ~isempty(anck) + A(j, anck) = 1; + end +end + +%%%%%%%% + +function A = init_ancestor_matrix(dag) + +order = topological_sort(dag); +A = zeros(length(dag)); +for j=order(:)' + A = update_row(A, j, dag); +end diff --git a/sourcecodes/bnt-master/BNT/learning/learn_params.m b/sourcecodes/bnt-master/BNT/learning/learn_params.m new file mode 100644 index 00000000..1bf9c843 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_params.m @@ -0,0 +1,26 @@ +function bnet = learn_params(bnet, data) +% LEARN_PARAMS Find the maximum likelihood params for a fully observed model +% bnet = learn_params(bnet, data) +% +% data(i,m) is the value of node i in case m (can be a cell array) +% +% We set bnet.CPD{i} to its ML/MAP estimate. +% +% Currently we assume no param tying + +% AND THAT EACH DATA POINT IS A SCALAR - no longer assumed + +%if iscell(data) +% data=cell2num(data); +%end +[n ncases] = size(data); +for j=1:n + e = bnet.equiv_class(j); + assert(e==j); + if adjustable_CPD(bnet.CPD{e}) + fam = family(bnet.dag,j); + %bnet.CPD{j} = learn_params(bnet.CPD{j}, data(fam,:)); + bnet.CPD{j} = learn_params(bnet.CPD{j}, fam, data, bnet.node_sizes, bnet.cnodes); + end +end + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_params_dbn.m b/sourcecodes/bnt-master/BNT/learning/learn_params_dbn.m new file mode 100644 index 00000000..4d0dc501 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_params_dbn.m @@ -0,0 +1,36 @@ +function bnet = learn_params_dbn(bnet, data) +% LEARN_PARAM_DBN Estimate params of a DBN for a fully observed model +% bnet = learn_params_dbn(bnet, data) +% +% data(i,t) is the value of node i in slice t (can be a cell array) +% We currently assume there is a single time series +% +% We set bnet.CPD{i} to its ML/MAP estimate. +% +% Currently we assume each node in the first 2 slices has its own CPD (no param tying); +% all nodes in slices >2 share their params with slice 2 as usual. + +[ss T] = size(data); + +% slice 1 +for j=1:ss + if adjustable_CPD(bnet.CPD{j}) + fam = family(bnet.dag,j); + bnet.CPD{j} = learn_params(bnet.CPD{j}, data(fam,1)); + end +end + + +% slices 2:T +% data2(:,t) contains [data(:,t-1); data(:,t)]. +% Then we extract out the rows corresponding to the parents in the current and previous slice. +data2 = [data(:,1:T-1); + data(:,2:T)]; +for j=1:ss + j2 = j+ss; + if adjustable_CPD(bnet.CPD{j2}) + fam = family(bnet.dag,j2); + bnet.CPD{j2} = learn_params(bnet.CPD{j2}, data2(fam,:)); + end +end + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_params_dbn_em.m b/sourcecodes/bnt-master/BNT/learning/learn_params_dbn_em.m new file mode 100644 index 00000000..e6f7f9f7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_params_dbn_em.m @@ -0,0 +1,179 @@ +function [bnet, LL, engine] = learn_params_dbn_em(engine, evidence, varargin) +% LEARN_PARAMS_DBN Set the parameters in a DBN to their ML/MAP values using batch EM. +% [bnet, LLtrace, engine] = learn_params_dbn_em(engine, data, ...) +% +% data{l}{i,t} = value of node i in slice t of time-series l, or [] if hidden. +% Suppose you have L time series, each of length T, in an O*T*L array D, +% where O is the num of observed scalar nodes, and N is the total num nodes per slice. +% Then you can create data as follows, where onodes is the index of the observable nodes: +% data = cell(1,L); +% for l=1:L +% data{l} = cell(N, T); +% data{l}(onodes,:) = num2cell(D(:,:,l)); +% end +% Of course it is possible for different sets of nodes to be observed in +% each slice/ sequence, and for each sequence to be a different length. +% +% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration. +% +% Optional arguments [default] +% +% max_iter - specifies the maximum number of iterations [100] +% thresh - specifies the thresold for stopping EM [1e-3] +% We stop when |f(t) - f(t-1)| / avg < threshold, +% where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik. +% verbose - display loglik at each iteration [1] +% anneal - 1 means do deterministic annealing (only for entropic priors) [0] +% anneal_rate - geometric cooling rate [0.8] +% init_temp - initial annealing temperature [10] +% final_temp - final annealing temperature [1e-3] +% + +max_iter = 100; +thresh = 1e-3; +anneal = 0; +anneal_rate = 0.8; +init_temp = 10; +final_temp = 1e-3; +verbose = 1; + +for i=1:2:length(varargin) + switch varargin{i} + case 'max_iter', max_iter = varargin{i+1}; + case 'thresh', thresh = varargin{i+1}; + case 'anneal', anneal = varargin{i+1}; + case 'anneal_rate', anneal_rate = varargin{i+1}; + case 'init_temp', init_temp = varargin{i+1}; + case 'final_temp', final_temp = varargin{i+1}; + otherwise, error(['unrecognized argument' varargin{i}]) + end +end + +% take 1 EM step at each temperature value, then when temp=0, run to convergence +% When using an entropic prior, Z = 1-T, so +% T=2 => Z=-1 (max entropy) +% T=1 => Z=0 (max likelihood) +% T=0 => Z=1 (min entropy / max structure) +num_iter = 1; +LL = []; +if anneal + temperature = init_temp; + while temperature > final_temp + [engine, loglik, logpost] = EM_step(engine, evidence, temperature); + if verbose + fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f, temp=%8.4f\n', ... + num_iter, loglik, logpost, temperature); + end + num_iter = num_iter + 1; + LL = [LL loglik]; + temperature = temperature * anneal_rate; + end + temperature = 0; + previous_loglik = loglik; + previous_logpost = logpost; +else + temperature = 0; + previous_loglik = -inf; + previous_logpost = -inf; +end + +converged = 0; +while ~converged & (num_iter <= max_iter) + [engine, loglik, logpost] = EM_step(engine, evidence, temperature); + if verbose + %fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f\n', ... + % num_iter, loglik, logpost); + fprintf('EM iteration %d, loglik = %8.4f\n', num_iter, loglik); + end + num_iter = num_iter + 1; + [converged, decreased] = em_converged(loglik, previous_loglik, thresh); + %[converged, decreased] = em_converged(logpost, previous_logpost, thresh); + previous_loglik = loglik; + previous_logpost = logpost; + LL = [LL loglik]; +end + +bnet = bnet_from_engine(engine); + +%%%%%%%%% + +function [engine, loglik, logpost] = EM_step(engine, cases, temp) + +bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step +ss = length(bnet.intra); +CPDs = bnet.CPD; % these are the new params that get maximized +num_CPDs = length(CPDs); + +% log P(theta|D) = (log P(D|theta) + log P(theta)) - log(P(D)) +% where log P(D|theta) = sum_cases log P(case|theta) +% and log P(theta) = sum_CPDs log P(CPD) - only count once even if tied! +% logpost = log P(theta,D) (un-normalized) +% This should be negative, and increase at every step. + +adjustable = zeros(1,num_CPDs); +logprior = zeros(1, num_CPDs); +for e=1:num_CPDs + adjustable(e) = adjustable_CPD(CPDs{e}); +end +adj = find(adjustable); + +for e=adj(:)' + logprior(e) = log_prior(CPDs{e}); + CPDs{e} = reset_ess(CPDs{e}); +end + +loglik = 0; +for l=1:length(cases) + evidence = cases{l}; + if ~iscell(evidence) + error('training data must be a cell array of cell arrays') + end + [engine, ll] = enter_evidence(engine, evidence); + assert(~isnan(ll)) + loglik = loglik + ll; + T = size(evidence, 2); + + % We unroll ns etc because in update_ess, we refer to nodes by their unrolled number + % so that they extract evidence from the right place. + % (The CPD should really store its own version of ns and cnodes...) + ns = repmat(bnet.node_sizes_slice(:), [1 T]); + cnodes = unroll_set(bnet.cnodes_slice, ss, T); + + %hidden_bitv = repmat(bnet.hidden_bitv(1:ss), [1 T]); + hidden_bitv = zeros(ss, T); + hidden_bitv(isemptycell(evidence))=1; + % hidden_bitv(i) = 1 means node i is hidden. + % We pass this in, rather than using isemptycell(evidence(dom)), because + % isemptycell is very slow. + + t = 1; + for i=1:ss + e = bnet.equiv_class(i,1); + if adjustable(e) + fmarg = marginal_family(engine, i, t); + CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:)); + end + end + + for i=1:ss + e = bnet.equiv_class(i,2); + if adjustable(e) + for t=2:T + fmarg = marginal_family(engine, i, t); + CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:)); + end + end + end +end + +logpost = loglik + sum(logprior(:)); + +for e=adj(:)' + CPDs{e} = maximize_params(CPDs{e}, temp); +end + +engine = update_engine(engine, CPDs); + + + + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_params_em.m b/sourcecodes/bnt-master/BNT/learning/learn_params_em.m new file mode 100644 index 00000000..f673e712 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_params_em.m @@ -0,0 +1,88 @@ +function [bnet, LL, engine] = learn_params_em(engine, evidence, max_iter, thresh) +% LEARN_PARAMS_EM Set the parameters of each adjustable node to their ML/MAP values using batch EM. +% [bnet, LLtrace, engine] = learn_params_em(engine, data, max_iter, thresh) +% +% data{i,l} is the value of node i in case l, or [] if hidden. +% Suppose you have L training cases in an O*L array, D, where O is the num observed +% scalar nodes, and N is the total num nodes. +% Then you can create 'data' as follows, where onodes is the index of the observable nodes: +% data = cell(N, L); +% data(onodes,:) = num2cell(D); +% Of course it is possible for different sets of nodes to be observed in each case. +% +% We return the modified bnet and engine. +% To see the learned parameters for node i, use the construct +% s = struct(bnet.CPD{i}); % violate object privacy +% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration. +% +% max_iter specifies the maximum number of iterations. Default: 10. +% +% thresh specifies the thresold for stopping EM. Default: 1e-3. +% We stop when |f(t) - f(t-1)| / avg < threshold, +% where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik. + +if nargin < 3, max_iter = 10; end +if nargin < 4, thresh = 1e-3; end + +verbose = 1; + +loglik = 0; +previous_loglik = -inf; +converged = 0; +num_iter = 1; +LL = []; + +while ~converged & (num_iter <= max_iter) + [engine, loglik] = EM_step(engine, evidence); + if verbose, fprintf('EM iteration %d, ll = %8.4f\n', num_iter, loglik); end + num_iter = num_iter + 1; + converged = em_converged(loglik, previous_loglik, thresh); + previous_loglik = loglik; + LL = [LL loglik]; +end +if verbose, fprintf('\n'); end + +bnet = bnet_from_engine(engine); + +%%%%%%%%% + +function [engine, loglik] = EM_step(engine, cases) + +bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step +CPDs = bnet.CPD; % these are the new params that get maximized +num_CPDs = length(CPDs); +adjustable = zeros(1,num_CPDs); +for e=1:num_CPDs + adjustable(e) = adjustable_CPD(CPDs{e}); +end +adj = find(adjustable); +n = length(bnet.dag); + +for e=adj(:)' + CPDs{e} = reset_ess(CPDs{e}); +end + +loglik = 0; +ncases = size(cases, 2); +for l=1:ncases + evidence = cases(:,l); + [engine, ll] = enter_evidence(engine, evidence); + loglik = loglik + ll; + hidden_bitv = zeros(1,n); + hidden_bitv(isemptycell(evidence))=1; + for i=1:n + e = bnet.equiv_class(i); + if adjustable(e) + fmarg = marginal_family(engine, i); + CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, bnet.node_sizes, bnet.cnodes, hidden_bitv); + end + end +end + +for e=adj(:)' + CPDs{e} = maximize_params(CPDs{e}); +end + +engine = update_engine(engine, CPDs); + + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_K2.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_K2.m new file mode 100644 index 00000000..8a01a2e9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_K2.m @@ -0,0 +1,104 @@ +function dag = learn_struct_K2(data, ns, order, varargin) +% LEARN_STRUCT_K2 Greedily learn the best structure compatible with a fixed node ordering +% best_dag = learn_struct_K2(data, node_sizes, order, ...) +% +% data(i,m) = value of node i in case m (can be a cell array). +% node_sizes(i) is the size of node i. +% order(i) is the i'th node in the topological ordering. +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% max_fan_in - this the largest number of parents we allow per node [N] +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% verbose - 'yes' means display output while running [ 'no' ] +% +% e.g., dag = learn_struct_K2(data, ns, order, 'scoring_fn', 'bic', 'params', []) +% +% To be backwards compatible with BNT2, you can also specify arguments as follows +% dag = learn_struct_K2(data, node_sizes, order, max_fan_in) +% +% This algorithm is described in +% - Cooper and Herskovits, "A Bayesian method for the induction of probabilistic +% networks from data", Machine Learning Journal 9:308--347, 1992 + +[n ncases] = size(data); + +% set default params +type = cell(1,n); +params = cell(1,n); +for i=1:n + type{i} = 'tabular'; + %params{i} = { 'prior', 1 }; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; +clamped = zeros(n, ncases); + +max_fan_in = n; +verbose = 0; + +args = varargin; +nargs = length(args); +if length(args) > 0 + if isstr(args{1}) + for i=1:2:nargs + switch args{i}, + case 'verbose', verbose = strcmp(args{i+1}, 'yes'); + case 'max_fan_in', max_fan_in = args{i+1}; + case 'scoring_fn', scoring_fn = args{i+1}; + case 'type', type = args{i+1}; + case 'discrete', discrete = args{i+1}; + case 'clamped', clamped = args{i+1}; + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + end + end + else + max_fan_in = args{1}; + end +end + +dag = zeros(n,n); + +for i=1:n + ps = []; + j = order(i); + u = find(clamped(j,:)==0); + score = score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); + if verbose, fprintf('\nnode %d, empty score %6.4f\n', j, score); end + done = 0; + while ~done & (length(ps) <= max_fan_in) + pps = mysetdiff(order(1:i-1), ps); % potential parents + nps = length(pps); + pscore = zeros(1, nps); + for pi=1:nps + p = pps(pi); + pscore(pi) = score_family(j, [ps p], type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); + if verbose, fprintf('considering adding %d to %d, score %6.4f\n', p, j, pscore(pi)); end + end + [best_pscore, best_p] = max(pscore); + best_p = pps(best_p); + if best_pscore > score + score = best_pscore; + ps = [ps best_p]; + if verbose, fprintf('* adding %d to %d, score %6.4f\n', best_p, j, best_pscore); end + else + done = 1; + end + end + if ~isempty(ps) % need this check for matlab 5.2 + dag(ps, j) = 1; + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_dbn_reveal.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_dbn_reveal.m new file mode 100644 index 00000000..7be4e5da --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_dbn_reveal.m @@ -0,0 +1,101 @@ +function inter = learn_struct_dbn_reveal(seqs, ns, max_fan_in, penalty) +% LEARN_STRUCT_DBN_REVEAL Learn inter-slice adjacency matrix given fully observable discrete time series +% inter = learn_struct_dbn_reveal(seqs, node_sizes, max_fan_in, penalty) +% +% seqs{l}{i,t} = value of node i in slice t of time-series l. +% If you have a single time series in an N*T array D, use +% seqs = { num2cell(D) }. +% If you have L time series, each of length T, in an N*T*L array D, use +% seqs= cell(1,L); for l=1:L, seqs{l} = num2cell(D(:,:,l)); end +% or, in vectorized form, +% seqs = squeeze(num2cell(num2cell(D),[1 2])); +% Currently the data is assumed to be discrete (1,2,...) +% +% node_sizes(i) is the number of possible values for node i +% max_fan_in is the largest number of parents we allow per node (default: N) +% penalty is weight given to the complexity penalty (default: 0.5) +% A penalty of 0.5 gives the BIC score. +% A penalty of 0 gives the ML score. +% Maximizing likelihood is equivalent to maximizing mutual information between parents and child. +% +% inter(i,j) = 1 iff node in slice t connects to node j in slice t+1 +% +% The parent set for each node in slice 2 is computed by evaluating all subsets of nodes in slice 1, +% and picking the largest scoring one. This takes O(n^k) time per node, where n is the num. nodes +% per slice, and k <= n is the max fan in. +% Since all the nodes are observed, we do not need to use an inference engine. +% And since we are only learning the inter-slice matrix, we do not need to check for cycles. +% +% This algorithm is described in +% - "REVEAL: A general reverse engineering algorithm for inference of genetic network +% architectures", Liang et al. PSB 1998 +% - "Extended dependency analysis of large systems", +% Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141 +% - "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998. + +n = length(ns); + +if nargin < 3, max_fan_in = n; end +if nargin < 4, penalty = 0.5; end + +inter = zeros(n,n); + +if ~iscell(seqs) + data{1} = seqs; +end + +nseq = length(seqs); +nslices = 0; +data = cell(1, nseq); +for l=1:nseq + nslices = nslices + size(seqs{l}, 2); + data{l} = cell2num(seqs{l})'; % each row is a case +end +ndata = nslices - nseq; % subtract off the initial slice of each sequence + +% We concatenate the sequences as in the following example. +% Let there be 2 sequences of lengths 4 and 5, with n nodes per slice, +% and let i be the target node. +% Then we construct following matrix D +% +% s{1}{1,1} ... s{1}{1,3} s{2}{1,1} ... s{2}{1,4} +% .... +% s{1}{n,1} ... s{1}{n,3} s{2}{n,1} ... s{2}{n,4} +% s{1}{i,2} ... s{1}{i,4} s{2}{i,2} ... s{2}{i,5} +% +% D(1:n, i) is the i'th input and D(n+1, i) is the i'th output. +% +% We concatenate each sequence separately to avoid treating the transition +% from the end of one sequence to the beginning of another as a "normal" transition. + + +for i=1:n + D = []; + for l=1:nseq + T = size(seqs{l}, 2); + A = cell2num(seqs{l}(:, 1:T-1)); + B = cell2num(seqs{l}(i, 2:T)); + C = [A;B]; + D = [D C]; + end + SS = subsets(1:n, max_fan_in, 1); % skip the empty set + nSS = length(SS); + bic_score = zeros(1, nSS); + ll_score = zeros(1, nSS); + target = n+1; + ns2 = [ns ns(i)]; + for h=1:nSS + ps = SS{h}; + dom = [ps target]; + counts = compute_counts(D(dom, :), ns2(dom)); + CPT = mk_stochastic(counts); + [bic_score(h), ll_score(h)] = bic_score_family(counts, CPT, ndata); + end + if penalty == 0 + h = argmax(ll_score); + else + h = argmax(bic_score); + end + ps = SS{h}; + inter(ps, i) = 1; +end diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_mcmc.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_mcmc.m new file mode 100644 index 00000000..d39af881 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_mcmc.m @@ -0,0 +1,284 @@ +function [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, varargin) +% MY_LEARN_STRUCT_MCMC Monte Carlo Markov Chain search over DAGs assuming fully observed data +% [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, ...) +% +% data(i,m) is the value of node i in case m. +% ns(i) is the number of discrete values node i can take on. +% +% sampled_graphs{m} is the m'th sampled graph. +% accept_ratio(t) = acceptance ratio at iteration t +% num_edges(t) = number of edges in model at iteration t +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% nsamples - number of samples to draw from the chain after burn-in [ 100*N ] +% burnin - number of steps to take before drawing samples [ 5*N ] +% init_dag - starting point for the search [ zeros(N,N) ] +% +% e.g., samples = my_learn_struct_mcmc(data, ns, 'nsamples', 1000); +% +% Modified by Sonia Leach (SML) 2/4/02, 9/5/03 + + + +[n ncases] = size(data); + + +% set default params +type = cell(1,n); +params = cell(1,n); +for i=1:n + type{i} = 'tabular'; + %params{i} = { 'prior', 1}; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; +clamped = zeros(n, ncases); +nsamples = 100*n; +burnin = 5*n; +dag = zeros(n); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'nsamples', nsamples = args{i+1}; + case 'burnin', burnin = args{i+1}; + case 'init_dag', dag = args{i+1}; + case 'scoring_fn', scoring_fn = args{i+1}; + case 'type', type = args{i+1}; + case 'discrete', discrete = args{i+1}; + case 'clamped', clamped = args{i+1}; + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + end +end + +% We implement the fast acyclicity check described by P. Giudici and R. Castelo, +% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001. + +% SML: also keep descendant matrix C +use_giudici = 1; +if use_giudici + [nbrs, ops, nodes, A] = mk_nbrs_of_digraph(dag); +else + [nbrs, ops, nodes] = mk_nbrs_of_dag(dag); + A = []; +end + +num_accepts = 1; +num_rejects = 1; +T = burnin + nsamples; +accept_ratio = zeros(1, T); +num_edges = zeros(1, T); +sampled_graphs = cell(1, nsamples); +%sampled_bitv = zeros(nsamples, n^2); + +for t=1:T + [dag, nbrs, ops, nodes, A, accept] = take_step(dag, nbrs, ops, ... + nodes, ns, data, clamped, A, ... + scoring_fn, discrete, type, params); + num_edges(t) = sum(dag(:)); + num_accepts = num_accepts + accept; + num_rejects = num_rejects + (1-accept); + accept_ratio(t) = num_accepts/num_rejects; + if t > burnin + sampled_graphs{t-burnin} = dag; + %sampled_bitv(t-burnin, :) = dag(:)'; + end +end + + +%%%%%%%%% + + +function [new_dag, new_nbrs, new_ops, new_nodes, A, accept] = ... + take_step(dag, nbrs, ops, nodes, ns, data, clamped, A, ... + scoring_fn, discrete, type, params, prior_w) + + +use_giudici = ~isempty(A); +if use_giudici + [new_dag, op, i, j, new_A] = pick_digraph_nbr(dag, nbrs, ops, nodes,A); % updates A + [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_digraph(new_dag, new_A); +else + d = sample_discrete(normalise(ones(1, length(nbrs)))); + new_dag = nbrs{d}; + op = ops{d}; + i = nodes(d, 1); j = nodes(d, 2); + [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_dag(new_dag); +end + +bf = bayes_factor(dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params); + +%R = bf * (new_prior / prior) * (length(nbrs) / length(new_nbrs)); +R = bf * (length(nbrs) / length(new_nbrs)); +u = rand(1,1); +if u > min(1,R) % reject the move + accept = 0; + new_dag = dag; + new_nbrs = nbrs; + new_ops = ops; + new_nodes = nodes; +else + accept = 1; + if use_giudici +A = new_A; % new_A already updated in pick_digraph_nbr + end +end + + +%%%%%%%%% + +function bfactor = bayes_factor(old_dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params) + +u = find(clamped(j,:)==0); +LLnew = score_family(j, parents(new_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); +LLold = score_family(j, parents(old_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); +bf1 = exp(LLnew - LLold); + +if strcmp(op, 'rev') % must also multiply in the changes to i's family + u = find(clamped(i,:)==0); + LLnew = score_family(i, parents(new_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i}); + LLold = score_family(i, parents(old_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i}); + bf2 = exp(LLnew - LLold); +else + bf2 = 1; +end +bfactor = bf1 * bf2; + + +%%%%%%%% Giudici stuff follows %%%%%%%%%% + + +% SML: This now updates A as it goes from digraph it choses +function [new_dag, op, i, j, new_A] = pick_digraph_nbr(dag, digraph_nbrs, ops, nodes, A) + +d = sample_discrete(normalise(ones(1, length(digraph_nbrs)))); +%d = myunidrnd(length(digraph_nbrs),1,1); +i = nodes(d, 1); j = nodes(d, 2); +new_dag = digraph_nbrs(:,:,d); +op = ops{d}; +new_A = update_ancestor_matrix(A, op, i, j, new_dag); + + +%%%%%%%%%%%%%% + + +function A = update_ancestor_matrix(A, op, i, j, dag) + +switch op +case 'add', + A = do_addition(A, op, i, j, dag); +case 'del', + A = do_removal(A, op, i, j, dag); +case 'rev', + A = do_removal(A, op, i, j, dag); + A = do_addition(A, op, j, i, dag); +end + + +%%%%%%%%%%%% + +function A = do_addition(A, op, i, j, dag) + +A(j,i) = 1; % i is an ancestor of j +anci = find(A(i,:)); +if ~isempty(anci) + A(j,anci) = 1; % all of i's ancestors are added to Anc(j) +end +ancj = find(A(j,:)); +descj = find(A(:,j)); +if ~isempty(ancj) + for k=descj(:)' + A(k,ancj) = 1; % all of j's ancestors are added to each descendant of j + end +end + +%%%%%%%%%%% +function A = do_removal(A, op, i, j, dag) + +% find all the descendants of j, and put them in topological order + +% SML: originally Kevin had the next line commented and the %* lines +% being used but I think this is equivalent and much less expensive +% I assume he put it there for debugging and never changed it back...? +descj = find(A(:,j)); +%* R = reachability_graph(dag); +%* descj = find(R(j,:)); + +order = topological_sort(dag); + +% SML: originally Kevin used the %* line but this was extracting the +% wrong things to sort +%* descj_topnum = order(descj); +[junk, perm] = sort(order); %SML:node i is perm(i)-TH in order +descj_topnum = perm(descj); %SML:descj(i) is descj_topnum(i)-th in order + +% SML: now re-sort descj by rank in descj_topnum +[junk, perm] = sort(descj_topnum); +descj = descj(perm); + +% Update j and all its descendants +A = update_row(A, j, dag); +for k = descj(:)' + A = update_row(A, k, dag); +end + +%%%%%%%%%%% + +function A = old_do_removal(A, op, i, j, dag) + +% find all the descendants of j, and put them in topological order +% SML: originally Kevin had the next line commented and the %* lines +% being used but I think this is equivalent and much less expensive +% I assume he put it there for debugging and never changed it back...? +descj = find(A(:,j)); +%* R = reachability_graph(dag); +%* descj = find(R(j,:)); + +order = topological_sort(dag); +descj_topnum = order(descj); +[junk, perm] = sort(descj_topnum); +descj = descj(perm); +% Update j and all its descendants +A = update_row(A, j, dag); +for k = descj(:)' + A = update_row(A, k, dag); +end + +%%%%%%%%% + +function A = update_row(A, j, dag) + +% We compute row j of A +A(j, :) = 0; +ps = parents(dag, j); +if ~isempty(ps) + A(j, ps) = 1; +end +for k=ps(:)' + anck = find(A(k,:)); + if ~isempty(anck) + A(j, anck) = 1; + end +end + +%%%%%%%% + +function A = init_ancestor_matrix(dag) + +order = topological_sort(dag); +A = zeros(length(dag)); +for j=order(:)' + A = update_row(A, j, dag); +end diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_ic_star.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_ic_star.m new file mode 100644 index 00000000..4bbba20e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_ic_star.m @@ -0,0 +1,155 @@ +function [pdag, G] = learn_struct_pdag_ic_star(cond_indep, n, k, varargin) +% LEARN_STRUCT_PDAG_IC_STAR Learn a partially oriented DAG (pattern) with latent +% variables using the IC* algorithm +% P = learn_struct_pdag_ic_star(cond_indep, n, k, ...) +% +% n is the number of nodes. +% k is an optional upper bound on the fan-in (default: n) +% cond_indep is a boolean function that will be called as follows: +% feval(cond_indep, x, y, S, ...) +% where x and y are nodes, and S is a set of nodes (positive integers), +% and ... are any optional parameters passed to this function. +% +% The output P is an adjacency matrix, in which +% P(i,j) = -1 if there is either a latent variable L such that i <-L-> j +% OR there is a directed edge from i->j. +% P(i,j) = -2 if there is a marked directed i-*>j edge. +% P(i,j) = P(j,i) = 1 if there is and undirected edge i--j +% P(i,j) = P(j,i) = 2 if there is a latent variable L such that i<-L->j. +% +% The IC* algorithm learns a latent structure associated with a set of observed +% variables. +% The latent structure revealed is the projection in which every latent variable is +% 1) a root node +% 2) linked to exactly two observed variables. +% Latent variables in the projection are represented using a bidirectional graph, +% and thus remain implicit. +% +% See Pearl, "Causality: Models, Reasoning, and Inference", 2000, p52 for more details. +% Written by Tamar Kushnir, 2000 + +sep = cell(n,n); +ord = 0; +done = 0; +G = ones(n,n); +G = setdiag(G,0); +while ~done + done = 1; + [X,Y] = find(G); + for i=1:length(X) + x = X(i); y = Y(i); + nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]); + if length(nbrs) >= ord & G(x,y) ~= 0 + done = 0; + SS = subsets(nbrs, ord, ord); % all subsets of size ord + for si=1:length(SS) + S = SS{si}; + if feval(cond_indep, x, y, S, varargin{:}) + G(x,y) = 0; + G(y,x) = 0; + sep{x,y} = myunion(sep{x,y}, S); + sep{y,x} = myunion(sep{y,x}, S); + break; % no need to check any more subsets + end + end + end + end + ord = ord + 1; +end + +% Create the minimal pattern, +% i.e., the only directed edges are V structures. +pdag = G; +[X, Y] = find(G); +% We want to generate all unique triples x,y,z +% where y is a common neighbor to x and z +for i=1:length(X) + x = X(i); + y = Y(i); + Z = find(G(y,:)); + Z = mysetdiff(Z, x); + for z=Z(:)' + if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x}) + pdag(x,y) = -1; pdag(y,x) = 0; + pdag(z,y) = -1; pdag(y,z) = 0; + end + end +end + +% Convert the minimal pattern to a complete one using the following rules: +% Rule 1: +% if a and b are non-adjacent nodes with a common neighbor c, +% if a->c and not b->c then c-*>b (marked arrow). +% Rule 2: +% if a and b are adjacent and there is a directed path (marked links) from a to b +% then a->b (add arrowhead). +%Pearl (2000) + +arrowin = [-1 -2 2]; +old_pdag = zeros(n); +iter = 0; +while ~isequal(pdag, old_pdag) + iter = iter + 1; + old_pdag = pdag; + % rule 1 + [X, Y] = find(pdag); + for i=1:length(X) + x = X(i); + y = Y(i); + Z = find(pdag(y,:)); + Z = mysetdiff(Z, x); + for z=Z(:)' + if G(x,z)==0 & ismember(pdag(x,y),arrowin) & ~ismember(pdag(z,y),arrowin) + pdag(y,z) = -2; pdag(z,y) = 0; + end + end + end + % rule 2 + [X, Y] = find(G); + %check all adjacent nodes because if pdag(x,y) = -1 + %and pdag(y,x) = 0 there could still be an bidirected edge between x & y. + for i=1:length(X) + x = X(i); + y = Y(i); + if ~ismember(pdag(x,y), arrowin) %x->y doesn't exist yet + %find marked path from x to y + add_arrow = marked_path(x,y,pdag); + if add_arrow + if pdag(y,x)==-1 %bidirected edge + pdag(x,y) = 2; pdag(y,x) = 2; + else + pdag(x,y) = -1;pdag(y,x) = 0; + end + end + end + end +end + + +%%%%%%%%%%%%% + +function t = marked_path(x,y,L) +% MARKED_PATH is a boolean function which returns 1 if a marked path +% between nodes x and y exists in the partially directed latent structure L. +% +% t = marked_path(x,y,L) +% +% x and y are the starting and ending nodes in the path, respectively. +% L is a latent structure (partially directed graph with possible latent variables). +% +% Rule 2 of IC* algorithm (see Pearl, 2000) + +t=0; + +%find set of marked links from x +marked = find(L(x,:)==-2); +if ismember(y,marked) + t=1; %marked path found +else + for m=marked(:)' + t = marked_path(m,y,L); + if t==1 + break; %stop when marked path found + end + end +end diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc.m new file mode 100644 index 00000000..d734a0bd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc.m @@ -0,0 +1,130 @@ +function [pdag, G] = learn_struct_pdag_pc(cond_indep, n, k, varargin) +% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm +% P = learn_struct_pdag_pc(cond_indep, n, k, ...) +% +% n is the number of nodes. +% k is an optional upper bound on the fan-in (default: n) +% cond_indep is a boolean function that will be called as follows: +% feval(cond_indep, x, y, S, ...) +% where x and y are nodes, and S is a set of nodes (positive integers), +% and ... are any optional parameters passed to this function. +% +% The output P is an adjacency matrix, in which +% P(i,j) = -1 if there is an i->j edge. +% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j +% +% The PC algorithm does structure learning assuming all variables are observed. +% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117. +% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in, +% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n). + + +sep = cell(n,n); +ord = 0; +done = 0; +G = ones(n,n); +G=setdiag(G,0); +while ~done + done = 1; + [X,Y] = find(G); + for i=1:length(X) + x = X(i); y = Y(i); + %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]); + nbrs = mysetdiff(neighbors(G, y), x); % bug fix by Raanan Yehezkel 6/27/04 + nbrs = unique(nbrs); % bug fix Immer Ebert ebert@tree.com 3/278/11 + if length(nbrs) >= ord & G(x,y) ~= 0 + done = 0; + %SS = subsets(nbrs, ord, ord); % all subsets of size ord + SS = subsets1(nbrs, ord); + for si=1:length(SS) + S = SS{si}; + if feval(cond_indep, x, y, S, varargin{:}) + %if isempty(S) + % fprintf('%d indep of %d ', x, y); + %else + % fprintf('%d indep of %d given ', x, y); fprintf('%d ', S); + %end + %fprintf('\n'); + + % diagnostic + %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:}); + %fprintf(': r = %6.4f\n', r); + + G(x,y) = 0; + G(y,x) = 0; + sep{x,y} = myunion(sep{x,y}, S); + sep{y,x} = myunion(sep{y,x}, S); + break; % no need to check any more subsets + end + end + end + end + ord = ord + 1; +end + + +% Create the minimal pattern, +% i.e., the only directed edges are V structures. +pdag = G; +[X, Y] = find(G); +% We want to generate all unique triples x,y,z +% This code generates x,y,z and z,y,x. +for i=1:length(X) + x = X(i); + y = Y(i); + Z = find(G(y,:)); + Z = mysetdiff(Z, x); + for z=Z(:)' + if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x}) + %fprintf('%d -> %d <- %d\n', x, y, z); + pdag(x,y) = -1; pdag(y,x) = 0; + pdag(z,y) = -1; pdag(y,z) = 0; + end + end +end + +% Convert the minimal pattern to a complete one, +% i.e., every directed edge in P is compelled +% (must be directed in all Markov equivalent models), +% and every undirected edge in P is reversible. +% We use the rules of Pearl (2000) p51 (derived in Meek (1995)) + +old_pdag = zeros(n); +iter = 0; +while ~isequal(pdag, old_pdag) + iter = iter + 1; + old_pdag = pdag; + % rule 1 + [A,B] = find(pdag==-1); % a -> b + for i=1:length(A) + a = A(i); b = B(i); + C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a + if ~isempty(C) + pdag(b,C) = -1; pdag(C,b) = 0; + %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C); + end + end + % rule 2 + [A,B] = find(pdag==1); % unoriented a-b edge + for i=1:length(A) + a = A(i); b = B(i); + if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' ); + pdag(a,b) = -1; pdag(b,a) = 0; + %fprintf('rule 2: %d -> %d\n', a, b); + end + end + % rule 3 + [A,B] = find(pdag==1); % a-b + for i=1:length(A) + a = A(i); b = B(i); + C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' ); + % C contains nodes c s.t. a-c->ba + G2 = setdiag(G(C, C), 1); + if any(G2(:)==0) % there are 2 different non adjacent elements of C + pdag(a,b) = -1; pdag(b,a) = 0; + %fprintf('rule 3: %d -> %d\n', a, b); + end + end +end + + diff --git a/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc_constrain.m b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc_constrain.m new file mode 100644 index 00000000..f226b87e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/learn_struct_pdag_pc_constrain.m @@ -0,0 +1,177 @@ +function [pdag, G] = dn_learn_struct_pdag_pc_constrain(adj, cond_indep, n, k, varargin) +% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm +% Pdag = learn_struct_pdag_pc_constrain(adj, cond_indep, n, k, ...) +% +% adj = adjacency matrix learned from dependency network P(i,j) = 1 => i--j; 0 => i j +% n is the number of nodes. +% k is an optional upper bound on the fan-in (default: n) +% cond_indep is a boolean function that will be called as follows: +% feval(cond_indep, x, y, S, ...) +% where x and y are nodes, and S is a set of nodes (positive integers), +% and ... are any optional parameters passed to this function. +% +%Output +% pdag Partially directed graph +% G Resulting adjacency graph prior to setting direction arrows +% +% The output P is an adjacency matrix, in which +% P(i,j) = -1 if there is an i->j edge. +% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j +% +% The PC algorithm does structure learning assuming all variables are observed. +% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117. +% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in, +% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n). +% +%% Example +%% Given data in a comma separated, filename starting with the variable labels, then the data in rows. +%% filename test.txt consists of: +%% +%% Earthquake,Burglar,Radio,Alarm,Call +%% 1,2,2,2,1 +%% 1,1,2,1,2 +%% . . . +%[CovMatrix, obs, varfields] = CovMat('test.txt',5); +% +%dn = zeros(5,5); +%dn(1,2) = 1; % This was the known Markov blanket of the system that generated test.txt +%dn(2,1) = 1; +%dn(2,4) = 1; +%dn(4,2) = 1; +%dn(1,3) = 1; +%dn(3,1) = 1; +%dn(1,4) = 1; +%dn(4,1) = 1; +%dn(4,5) = 1; +%dn(5,4) = 1; +%dn(3,5) = 1; %loop r->c +%dn(5,3) = 1; %loop c-r +%dn(3,4) = 1; +%dn(4,3) = 1; +% +%max_fan_in = 4; +%alpha = 0.05; +% +%[pdag G] = learn_struct_pdag_pc_constrain(dn,'cond_indep_fisher_z', 5, max_fan_in, CovMatrix, obs, alpha); +%% +%% +%% Gary Bradski, 7/2002 Modified this to take an adjacency matrix from a dependency network. + + +sep = cell(n,n); +ord = 0; +done = 0; +G = ones(n,n); +G=setdiag(G,0); + +while ~done + done = 1; + [X,Y] = find(G); + for i=1:length(X) + x = X(i); y = Y(i); +% nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);%parents, children, but not self + nbrs = mysetdiff(myunion(neighbors(adj, x), neighbors(adj,y)), [x y]);%parents, children, but not self + + if length(nbrs) >= ord & G(x,y) ~= 0 + done = 0; + SS = subsets(nbrs, ord, ord); % all subsets of size ord + for si=1:length(SS) + S = SS{si}; + %if (feval(dsep,x,y,S,adj)) | (feval(cond_indep, x, y, S, varargin{:})) + if feval(cond_indep, x, y, S, varargin{:}) + %if isempty(S) + % fprintf('%d indep of %d ', x, y); + %else + % fprintf('%d indep of %d given ', x, y); fprintf('%d ', S); + %end + %fprintf('\n'); + + % diagnostic + %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:}); + %fprintf(': r = %6.4f\n', r); + + G(x,y) = 0; + G(y,x) = 0; + adj(x,y) = 0; %make sure found cond. independencies are marked out + adj(y,x) = 0; + sep{x,y} = myunion(sep{x,y}, S); + sep{y,x} = myunion(sep{y,x}, S); + break; % no need to check any more subsets + end + end + end + end + ord = ord + 1; +end + + + + +% Create the minimal pattern, +% i.e., the only directed edges are V structures. + +pdag = G; +[X, Y] = find(G); +% We want to generate all unique triples x,y,z +% This code generates x,y,z and z,y,x. +for i=1:length(X) + x = X(i); + y = Y(i); + Z = find(G(y,:)); + Z = mysetdiff(Z, x); + for z=Z(:)' + if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x}) + %fprintf('%d -> %d <- %d\n', x, y, z); + pdag(x,y) = -1; pdag(y,x) = 0; + pdag(z,y) = -1; pdag(y,z) = 0; + end + end +end + +% Convert the minimal pattern to a complete one, +% i.e., every directed edge in P is compelled +% (must be directed in all Markov equivalent models), +% and every undirected edge in P is reversible. +% We use the rules of Pearl (2000) p51 (derived in Meek (1995)) + +old_pdag = zeros(n); +iter = 0; +while ~isequal(pdag, old_pdag) + iter = iter + 1; + old_pdag = pdag; + % rule 1 + [A,B] = find(pdag==-1); % a -> b + for i=1:length(A) + a = A(i); b = B(i); + C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a + if ~isempty(C) + pdag(b,C) = -1; pdag(C,b) = 0; + %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C); + end + end + % rule 2 + [A,B] = find(pdag==1); % unoriented a-b edge + for i=1:length(A) + a = A(i); b = B(i); + if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' ); + pdag(a,b) = -1; pdag(b,a) = 0; + %fprintf('rule 2: %d -> %d\n', a, b); + end + end + % rule 3 + [A,B] = find(pdag==1); % a-b + for i=1:length(A) + a = A(i); b = B(i); + C = find( (G(a,:)==1) & (pdag(:,b)==-1)' ); + % C contains nodes c s.t. a-c->ba + G2 = setdiag(G(C, C), 1); + if any(G2(:)==0) % there are 2 different non adjacent elements of C + pdag(a,b) = -1; pdag(b,a) = 0; + %fprintf('rule 3: %d -> %d\n', a, b); + end + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/learning/mcmc_sample_to_hist.m b/sourcecodes/bnt-master/BNT/learning/mcmc_sample_to_hist.m new file mode 100644 index 00000000..8b277c8a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/mcmc_sample_to_hist.m @@ -0,0 +1,36 @@ +function mcmc_post = mcmc_sample_to_hist(sampled_graphs, dags) +% MCMC_SAMPLE_TO_HIST Convert a set of sampled dags into a histogram over dags +% hist = mcmc_sample_to_hist(sampled_graphs, dags) +% +% sampled_graphs{m} is the m'th sampled dag +% dags{i} is the i'th dag in the hypothesis space +% hist(i) = Pr(model i | data) + +ndags = length(dags); +nsamples = length(sampled_graphs); +nnodes = length(dags{1}); +% sampled_bitv(m, :) is the m'th sampled graph represented as a vector of n^2 bits, computed +% by stacking the columns of the adjacency matrix vertically. +sampled_bitvs = zeros(nsamples, nnodes*nnodes); +for m=1:nsamples + sampled_bitvs(m, :) = sampled_graphs{m}(:)'; +end + +[ugraphs, I, J] = unique(sampled_bitvs, 'rows'); % each row of ugraphs is a unique bit vector +sampled_indices = subv2ind(2*ones(1,nnodes*nnodes), ugraphs+1); +counts = hist(J, 1:size(ugraphs,1)); % counts(i) = number of times graphs(i,:) occurs in the sample + +mcmc_post = zeros(1, ndags); +for i=1:ndags + bitv = dags{i}(:)'; + % Find the samples that corresponds to this graph by converting the graphs to bitvectors and + % then to integers. + ndx = subv2ind(2*ones(1,nnodes*nnodes), bitv+1); + locn = find(ndx == sampled_indices); + if ~isempty(locn) + mcmc_post(i) = counts(locn); + end +end +mcmc_post = normalise(mcmc_post); + + diff --git a/sourcecodes/bnt-master/BNT/learning/mk_schedule.m b/sourcecodes/bnt-master/BNT/learning/mk_schedule.m new file mode 100644 index 00000000..018f3d16 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/mk_schedule.m @@ -0,0 +1,10 @@ +function schedule = mk_schedule(init_temp, final_temp, anneal_rate) + +init_temp = 10; final_temp = 1e-2; anneal_rate = 0.8; +schedule = []; +temp=init_temp; +schedule = [schedule temp]; +while temp > final_temp + temp = temp * anneal_rate; + schedule = [schedule temp]; +end diff --git a/sourcecodes/bnt-master/BNT/learning/mk_tetrad_data_file.m b/sourcecodes/bnt-master/BNT/learning/mk_tetrad_data_file.m new file mode 100644 index 00000000..9c556556 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/mk_tetrad_data_file.m @@ -0,0 +1,25 @@ +function mk_tetrad_data_file(filename, samples, sig) +% MK_TETRAD_DATA_FILE Make a file containing raw discrete data for input to TETRAD +% mk_tetrad_data_file(filename, samples, sig) +% +% samples(i,j) is the value for case i, variable j +% The resulting file can be used for the 'build' part of Tetrad. +% For details on tetrad, see hss.cmu.edu/html/departments/philosophy/TETRAD/tetrad.html + +[nsamples N] = size(samples); + +fid = fopen(filename, 'w'); +fprintf(fid, '/Raw\n'); +fprintf(fid, '%d\n', nsamples); +for i=1:N + fprintf(fid, 'x%d ', i); +end +fprintf(fid, '\n'); +for i=1:nsamples + fprintf(fid, '%d ', samples(i,:)-1); % tetrad counts from 0 + fprintf(fid, '\n'); +end +%fprintf(fid, '/Knowledge\n'); +%fprintf(fid, 'Significance %4.2f\n', sig); +fclose(fid); + diff --git a/sourcecodes/bnt-master/BNT/learning/score_dags_old.m b/sourcecodes/bnt-master/BNT/learning/score_dags_old.m new file mode 100644 index 00000000..507721f6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/score_dags_old.m @@ -0,0 +1,62 @@ +function score = score_dags(data, ns, dags, varargin) +% SCORE_DAGS Compute the score of one or more DAGs +% score = score_dags(data, ns, dags, varargin) +% +% data{i,m} = value of node i in case m (can be a cell array). +% node_sizes(i) is the number of size of node i. +% dags{g} is the g'th dag +% score(g) is the score of the i'th dag +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% +% e.g., score = score_dags(data, ns, mk_all_dags(n), 'scoring_fn', 'bic', 'params', []); +% +% If the DAGs have a lot of families in common, we can cache the sufficient statistics, +% making this potentially more efficient than scoring the DAGs one at a time. +% (Caching is not currently implemented, however.) + +[n ncases] = size(data); + +% set default params +type = cell(1,n); +params = cell(1,n); +for i=1:n + type{i} = 'tabular'; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; +clamped = zeros(n, ncases); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'scoring_fn', scoring_fn = args{i+1}; + case 'type', type = args{i+1}; + case 'discrete', discrete = args{i+1}; + case 'clamped', clamped = args{i+1}; + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + end +end + +NG = length(dags); +score = zeros(1, NG); +for g=1:NG + dag = dags{g}; + for j=1:n + u = find(clamped(j,:)==0); + ps = parents(dag, j); + score(g) = score(g) + score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); + end +end diff --git a/sourcecodes/bnt-master/BNT/learning/score_dags_wocache.m b/sourcecodes/bnt-master/BNT/learning/score_dags_wocache.m new file mode 100644 index 00000000..79684080 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/score_dags_wocache.m @@ -0,0 +1,67 @@ +function score = score_dags(data, ns, dags, varargin) +% SCORE_DAGS Compute the score of one or more DAGs +% score = score_dags(data, ns, dags, varargin) +% +% data{i,m} = value of node i in case m (can be a cell array). +% node_sizes(i) is the number of size of node i. +% dags{g} is the g'th dag +% score(g) is the score of the i'th dag +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% +% e.g., score = score_dags(data, ns, mk_all_dags(n), 'scoring_fn', 'bic', 'params', []); +% +% If the DAGs have a lot of families in common, we can cache the sufficient statistics, +% making this potentially more efficient than scoring the DAGs one at a time. +% (Caching is not currently implemented, however.) + +[n ncases] = size(data); + +% set default params +type = cell(1,n); +params = cell(1,n); +for i=1:n + type{i} = 'tabular'; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; + +u = [1:ncases]'; % DWH +isclamped = 0; %DWH +clamped = zeros(n, ncases); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'scoring_fn', scoring_fn = args{i+1}; + case 'type', type = args{i+1}; + case 'discrete', discrete = args{i+1}; + case 'clamped', clamped = args{i+1}, isclamped = 1; %DWH + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + end +end + +NG = length(dags); +score = zeros(1, NG); +for g=1:NG + dag = dags{g}; + for j=1:n + if isclamped %DWH + u = find(clamped(j,:)==0); + end + ps = parents(dag, j); + score(g) = score(g) + score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j}); + end +end diff --git a/sourcecodes/bnt-master/BNT/learning/score_family_wocache.m b/sourcecodes/bnt-master/BNT/learning/score_family_wocache.m new file mode 100644 index 00000000..a0710ba5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/learning/score_family_wocache.m @@ -0,0 +1,57 @@ +function score = score_family(j, ps, node_type, scoring_fn, ns, discrete, data, args) +% SCORE_FAMILY_COMPLETE Compute the score of a node and its parents given completely observed data +% score = score_family(j, ps, node_type, scoring_fn, ns, discrete, data, args) +% +% data(i,m) is the value of node i in case m (can be a cell array) +% args is a cell array containing optional arguments passed to the constructor, +% or is [] if none +% +% We create a whole Bayes net which only connects parents to node, +% where node has a CPD of the specified type (with default parameters). +% We then evaluate its score ('bic' or 'bayesian') + +% We should use a cache to avoid unnecessary computation. +% In particular, log_marginal_prob_node for tabular CPDs calls gammaln +% and compute_counts, both of which are slow. + +[n ncases] = size(data); +dag = zeros(n,n); +% SML added to sort ps b/c mk_bnet, learn_params use sorted ps to make +% CPTs +% Kevin had: if ~isempty(ps), dag(ps, j) = 1; end +if ~isempty(ps), dag(ps, j) = 1;, ps = sort(ps);, end + +bnet = mk_bnet(dag, ns, 'discrete', discrete); +%bnet.CPD{j} = xxx_CPD(bnet, j); +%eval(sprintf('bnet.CPD{j} = %s_CPD(bnet, j);', node_type)); +fname = sprintf('%s_CPD', node_type); +%fprintf('score CPD %d\n', j); +if isempty(args) + bnet.CPD{j} = feval(fname, bnet, j); +else + bnet.CPD{j} = feval(fname, bnet, j, args{:}); +end +switch scoring_fn + case 'bic', + fam = [ps j]; + %score = BIC_score_CPD(bnet.CPD{j}, fam, data, ns, bnet.cnodes); + %bnet.CPD{j} = learn_params(bnet.CPD{j}, fam, data, ns, bnet.cnodes); + + % SML 03/16/04 had to special case gaussian b/c generic_CPD/learn_params + % no longer supported because of simple interface to learn_params + % introduced by KPM for tabular nodes below: + % KPM 9 June 04 - tabular nodes have changed back! + if 1 % (isempty(find(j==discrete))) + bnet.CPD{j} = learn_params(bnet.CPD{j}, fam, data, ns, bnet.cnodes); + else + bnet.CPD{j} = learn_params(bnet.CPD{j}, data(fam, :)); + end + L = log_prob_node(bnet.CPD{j}, data(j,:), data(ps,:)); + S = struct(bnet.CPD{j}); % violate object privacy + score = L - 0.5*S.nparams*log(ncases); + case 'bayesian', + %score = bayesian_score_CPD(bnet.CPD{j}, data(fam, :)); + score = log_marg_prob_node(bnet.CPD{j}, data(j,:), data(ps,:)); + otherwise, + error(['unrecognized scoring fn ' scoring_fn]); +end -- cgit 1.4.1