diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/learning')
25 files changed, 1982 insertions, 0 deletions
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 <raanany@ee.bgu.ac.il> 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 |