|
wolffd@0
|
1 function dag = learn_struct_K2(data, ns, order, varargin)
|
|
wolffd@0
|
2 % LEARN_STRUCT_K2 Greedily learn the best structure compatible with a fixed node ordering
|
|
wolffd@0
|
3 % best_dag = learn_struct_K2(data, node_sizes, order, ...)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % data(i,m) = value of node i in case m (can be a cell array).
|
|
wolffd@0
|
6 % node_sizes(i) is the size of node i.
|
|
wolffd@0
|
7 % order(i) is the i'th node in the topological ordering.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % The following optional arguments can be specified in the form of name/value pairs:
|
|
wolffd@0
|
10 % [default value in brackets]
|
|
wolffd@0
|
11 %
|
|
wolffd@0
|
12 % max_fan_in - this the largest number of parents we allow per node [N]
|
|
wolffd@0
|
13 % scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
|
|
wolffd@0
|
14 % Currently, only networks with all tabular nodes support Bayesian scoring.
|
|
wolffd@0
|
15 % type - type{i} is the type of CPD to use for node i, where the type is a string
|
|
wolffd@0
|
16 % of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
|
|
wolffd@0
|
17 % params - params{i} contains optional arguments passed to the CPD constructor for node i,
|
|
wolffd@0
|
18 % or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
|
|
wolffd@0
|
19 % discrete - the list of discrete nodes [ 1:N ]
|
|
wolffd@0
|
20 % clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
|
|
wolffd@0
|
21 % verbose - 'yes' means display output while running [ 'no' ]
|
|
wolffd@0
|
22 %
|
|
wolffd@0
|
23 % e.g., dag = learn_struct_K2(data, ns, order, 'scoring_fn', 'bic', 'params', [])
|
|
wolffd@0
|
24 %
|
|
wolffd@0
|
25 % To be backwards compatible with BNT2, you can also specify arguments as follows
|
|
wolffd@0
|
26 % dag = learn_struct_K2(data, node_sizes, order, max_fan_in)
|
|
wolffd@0
|
27 %
|
|
wolffd@0
|
28 % This algorithm is described in
|
|
wolffd@0
|
29 % - Cooper and Herskovits, "A Bayesian method for the induction of probabilistic
|
|
wolffd@0
|
30 % networks from data", Machine Learning Journal 9:308--347, 1992
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32 [n ncases] = size(data);
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % set default params
|
|
wolffd@0
|
35 type = cell(1,n);
|
|
wolffd@0
|
36 params = cell(1,n);
|
|
wolffd@0
|
37 for i=1:n
|
|
wolffd@0
|
38 type{i} = 'tabular';
|
|
wolffd@0
|
39 %params{i} = { 'prior', 1 };
|
|
wolffd@0
|
40 params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
|
|
wolffd@0
|
41 end
|
|
wolffd@0
|
42 scoring_fn = 'bayesian';
|
|
wolffd@0
|
43 discrete = 1:n;
|
|
wolffd@0
|
44 clamped = zeros(n, ncases);
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46 max_fan_in = n;
|
|
wolffd@0
|
47 verbose = 0;
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 args = varargin;
|
|
wolffd@0
|
50 nargs = length(args);
|
|
wolffd@0
|
51 if length(args) > 0
|
|
wolffd@0
|
52 if isstr(args{1})
|
|
wolffd@0
|
53 for i=1:2:nargs
|
|
wolffd@0
|
54 switch args{i},
|
|
wolffd@0
|
55 case 'verbose', verbose = strcmp(args{i+1}, 'yes');
|
|
wolffd@0
|
56 case 'max_fan_in', max_fan_in = args{i+1};
|
|
wolffd@0
|
57 case 'scoring_fn', scoring_fn = args{i+1};
|
|
wolffd@0
|
58 case 'type', type = args{i+1};
|
|
wolffd@0
|
59 case 'discrete', discrete = args{i+1};
|
|
wolffd@0
|
60 case 'clamped', clamped = args{i+1};
|
|
wolffd@0
|
61 case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end
|
|
wolffd@0
|
62 end
|
|
wolffd@0
|
63 end
|
|
wolffd@0
|
64 else
|
|
wolffd@0
|
65 max_fan_in = args{1};
|
|
wolffd@0
|
66 end
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 dag = zeros(n,n);
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 for i=1:n
|
|
wolffd@0
|
72 ps = [];
|
|
wolffd@0
|
73 j = order(i);
|
|
wolffd@0
|
74 u = find(clamped(j,:)==0);
|
|
wolffd@0
|
75 score = score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
|
|
wolffd@0
|
76 if verbose, fprintf('\nnode %d, empty score %6.4f\n', j, score); end
|
|
wolffd@0
|
77 done = 0;
|
|
wolffd@0
|
78 while ~done & (length(ps) <= max_fan_in)
|
|
wolffd@0
|
79 pps = mysetdiff(order(1:i-1), ps); % potential parents
|
|
wolffd@0
|
80 nps = length(pps);
|
|
wolffd@0
|
81 pscore = zeros(1, nps);
|
|
wolffd@0
|
82 for pi=1:nps
|
|
wolffd@0
|
83 p = pps(pi);
|
|
wolffd@0
|
84 pscore(pi) = score_family(j, [ps p], type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
|
|
wolffd@0
|
85 if verbose, fprintf('considering adding %d to %d, score %6.4f\n', p, j, pscore(pi)); end
|
|
wolffd@0
|
86 end
|
|
wolffd@0
|
87 [best_pscore, best_p] = max(pscore);
|
|
wolffd@0
|
88 best_p = pps(best_p);
|
|
wolffd@0
|
89 if best_pscore > score
|
|
wolffd@0
|
90 score = best_pscore;
|
|
wolffd@0
|
91 ps = [ps best_p];
|
|
wolffd@0
|
92 if verbose, fprintf('* adding %d to %d, score %6.4f\n', best_p, j, best_pscore); end
|
|
wolffd@0
|
93 else
|
|
wolffd@0
|
94 done = 1;
|
|
wolffd@0
|
95 end
|
|
wolffd@0
|
96 end
|
|
wolffd@0
|
97 if ~isempty(ps) % need this check for matlab 5.2
|
|
wolffd@0
|
98 dag(ps, j) = 1;
|
|
wolffd@0
|
99 end
|
|
wolffd@0
|
100 end
|
|
wolffd@0
|
101
|
|
wolffd@0
|
102
|
|
wolffd@0
|
103
|
|
wolffd@0
|
104
|
