wolffd@0
|
1 function bnet = mk_higher_order_dbn(intra, inter, node_sizes, varargin)
|
wolffd@0
|
2 % MK_DBN Make a Dynamic Bayesian Network.
|
wolffd@0
|
3 %
|
wolffd@0
|
4 % BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
|
wolffd@0
|
5 % from i in slice t to j in slice t iff intra(i,j) = 1, and
|
wolffd@0
|
6 % from i in slice t to j in slice t+1 iff inter(i,j) = 1,
|
wolffd@0
|
7 % for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
|
wolffd@0
|
8 % node_sizes(i) is the number of values node i can take on.
|
wolffd@0
|
9 % The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
|
wolffd@0
|
10 % See also mk_bnet.
|
wolffd@0
|
11 %
|
wolffd@0
|
12 % Optional arguments [default in brackets]
|
wolffd@0
|
13 % 'discrete' - list of discrete nodes [1:n]
|
wolffd@0
|
14 % 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
|
wolffd@0
|
15 % 'eclass1' - equiv class for slice 1 [1:n]
|
wolffd@0
|
16 % 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
|
wolffd@0
|
17 % equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
|
wolffd@0
|
18 % i.e., nodes i and j have tied parameters.
|
wolffd@0
|
19 % 'intra1' - topology of first slice, if different from others
|
wolffd@0
|
20 % 'names' - a cell array of strings to be associated with nodes 1:n [{}]
|
wolffd@0
|
21 % This creates an associative array, so you write e.g.
|
wolffd@0
|
22 % 'evidence(bnet.names{'bar'}) = 42' instead of 'evidence(2} = 42'
|
wolffd@0
|
23 % assuming names = { 'foo', 'bar', ...}.
|
wolffd@0
|
24 %
|
wolffd@0
|
25 % For backwards compatibility with BNT2, arguments can also be specified as follows
|
wolffd@0
|
26 % bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
|
wolffd@0
|
27 %
|
wolffd@0
|
28 % After calling this function, you must specify the parameters (conditional probability
|
wolffd@0
|
29 % distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.
|
wolffd@0
|
30
|
wolffd@0
|
31
|
wolffd@0
|
32 n = length(intra);
|
wolffd@0
|
33 ss = n;
|
wolffd@0
|
34 bnet.nnodes_per_slice = ss;
|
wolffd@0
|
35 bnet.intra = intra;
|
wolffd@0
|
36 bnet.inter = inter;
|
wolffd@0
|
37 bnet.intra1 = intra;
|
wolffd@0
|
38
|
wolffd@0
|
39 % As this method is used to generate a higher order Markov Model
|
wolffd@0
|
40 % also connect from time slice t - i -> t with i > 1 has to be
|
wolffd@0
|
41 % taken into account.
|
wolffd@0
|
42
|
wolffd@0
|
43 %inter should be a three dimensional array where inter(:,:,i)
|
wolffd@0
|
44 %describes the connections from time-slice t - i to t.
|
wolffd@0
|
45 [rows,columns,order] = size(inter);
|
wolffd@0
|
46 assert(rows == n);
|
wolffd@0
|
47 assert(columns == n);
|
wolffd@0
|
48 dag = zeros((order + 1)*n);
|
wolffd@0
|
49
|
wolffd@0
|
50 i = 0;
|
wolffd@0
|
51 while i <= order
|
wolffd@0
|
52 j = i;
|
wolffd@0
|
53 while j <= order
|
wolffd@0
|
54 if j == i
|
wolffd@0
|
55 dag(1 + i*n:(i+1)*n,1+i*n:(i+1)*n) = intra;
|
wolffd@0
|
56 else
|
wolffd@0
|
57 dag(1+i*n:(i+1)*n,1+j*n:(j+1)*n) = inter(:,:,j - i);
|
wolffd@0
|
58 end
|
wolffd@0
|
59 j = j + 1;
|
wolffd@0
|
60 end;
|
wolffd@0
|
61 i = i + 1;
|
wolffd@0
|
62 end;
|
wolffd@0
|
63
|
wolffd@0
|
64 bnet.dag = dag;
|
wolffd@0
|
65 bnet.names = {};
|
wolffd@0
|
66
|
wolffd@0
|
67 directed = 1;
|
wolffd@0
|
68 if ~acyclic(dag,directed)
|
wolffd@0
|
69 error('graph must be acyclic')
|
wolffd@0
|
70 end
|
wolffd@0
|
71
|
wolffd@0
|
72 % Calculation of the equivalence classes
|
wolffd@0
|
73 bnet.eclass1 = 1:n;
|
wolffd@0
|
74 bnet.eclass = zeros(order + 1,ss);
|
wolffd@0
|
75 bnet.eclass(1,:) = 1:n;
|
wolffd@0
|
76 for i = 1:order
|
wolffd@0
|
77 bnet.eclass(i+1,:) = bnet.eclass(i,:);
|
wolffd@0
|
78 for j = 1:ss
|
wolffd@0
|
79 if(isequal(parents(dag,(i-1)*n+j)+ss,parents(dag,(i*n + j))))
|
wolffd@0
|
80 %fprintf('%d has isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))
|
wolffd@0
|
81 else
|
wolffd@0
|
82 bnet.eclass(i + 1,j) = max(bnet.eclass(i+1,:))+1;
|
wolffd@0
|
83 %fprintf('%d has non isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))
|
wolffd@0
|
84 end;
|
wolffd@0
|
85 end;
|
wolffd@0
|
86 end;
|
wolffd@0
|
87 bnet.eclass1 = 1:n;
|
wolffd@0
|
88
|
wolffd@0
|
89 % To be compatible with whe rest of the code
|
wolffd@0
|
90 bnet.eclass2 = bnet.eclass(2,:);
|
wolffd@0
|
91
|
wolffd@0
|
92 dnodes = 1:n;
|
wolffd@0
|
93 bnet.observed = [];
|
wolffd@0
|
94
|
wolffd@0
|
95 if nargin >= 4
|
wolffd@0
|
96 args = varargin;
|
wolffd@0
|
97 nargs = length(args);
|
wolffd@0
|
98 if ~isstr(args{1})
|
wolffd@0
|
99 if nargs >= 1 dnodes = args{1}; end
|
wolffd@0
|
100 if nargs >= 2 bnet.eclass1 = args{2}; bnet.eclass(1,:) = args{2}; end
|
wolffd@0
|
101 if nargs >= 3 bnet.eclass2 = args{3}; bnet.eclass(2,:) = args{2}; end
|
wolffd@0
|
102 if nargs >= 4 bnet.intra1 = args{4}; end
|
wolffd@0
|
103 else
|
wolffd@0
|
104 for i=1:2:nargs
|
wolffd@0
|
105 switch args{i},
|
wolffd@0
|
106 case 'discrete', dnodes = args{i+1};
|
wolffd@0
|
107 case 'observed', bnet.observed = args{i+1};
|
wolffd@0
|
108 case 'eclass1', bnet.eclass1 = args{i+1}; bnet.eclass(1,:) = args{i+1};
|
wolffd@0
|
109 case 'eclass2', bnet.eclass2 = args{i+1}; bnet.eclass(2,:) = args{i+1};
|
wolffd@0
|
110 case 'eclass', bnet.eclass = args{i+1};
|
wolffd@0
|
111 case 'intra1', bnet.intra1 = args{i+1};
|
wolffd@0
|
112 %case 'ar_hmm', bnet.ar_hmm = args{i+1}; % should check topology
|
wolffd@0
|
113 case 'names', bnet.names = assocarray(args{i+1}, num2cell(1:n));
|
wolffd@0
|
114 otherwise,
|
wolffd@0
|
115 error(['invalid argument name ' args{i}]);
|
wolffd@0
|
116 end
|
wolffd@0
|
117 end
|
wolffd@0
|
118 end
|
wolffd@0
|
119 end
|
wolffd@0
|
120
|
wolffd@0
|
121 bnet.observed = sort(bnet.observed); % for comparing sets
|
wolffd@0
|
122 ns = node_sizes;
|
wolffd@0
|
123 bnet.node_sizes_slice = ns(:)';
|
wolffd@0
|
124 bnet.node_sizes = repmat(ns(:),1,order + 1);
|
wolffd@0
|
125
|
wolffd@0
|
126 cnodes = mysetdiff(1:n, dnodes);
|
wolffd@0
|
127 bnet.dnodes_slice = dnodes;
|
wolffd@0
|
128 bnet.cnodes_slice = cnodes;
|
wolffd@0
|
129 bnet.dnodes = dnodes;
|
wolffd@0
|
130 bnet.cnodes = cnodes;
|
wolffd@0
|
131 % To adapt the function to higher order Markov models include dnodes for more
|
wolffd@0
|
132 % time slices
|
wolffd@0
|
133 for i = 1:order
|
wolffd@0
|
134 bnet.dnodes = [bnet.dnodes dnodes+i*n];
|
wolffd@0
|
135 bnet.cnodes = [bnet.cnodes cnodes+i*n];
|
wolffd@0
|
136 end
|
wolffd@0
|
137
|
wolffd@0
|
138 % Generieren einer Matrix, deren i-te Spalte die Aequivalenzklassen
|
wolffd@0
|
139 % der i-ten Zeitscheibe enthaelt.
|
wolffd@0
|
140 bnet.equiv_class = [bnet.eclass(1,:)]';
|
wolffd@0
|
141 for i = 2:(order + 1)
|
wolffd@0
|
142 bnet.equiv_class = [bnet.equiv_class bnet.eclass(i,:)'];
|
wolffd@0
|
143 end
|
wolffd@0
|
144
|
wolffd@0
|
145 bnet.CPD = cell(1,max(bnet.equiv_class(:)));
|
wolffd@0
|
146
|
wolffd@0
|
147 ss = n;
|
wolffd@0
|
148 onodes = bnet.observed;
|
wolffd@0
|
149 hnodes = mysetdiff(1:ss, onodes);
|
wolffd@0
|
150 bnet.hidden_bitv = zeros(1,(order + 1)*ss);
|
wolffd@0
|
151 for i = 0:order
|
wolffd@0
|
152 bnet.hidden_bitv(hnodes +i*ss) = 1;
|
wolffd@0
|
153 end;
|
wolffd@0
|
154
|
wolffd@0
|
155 bnet.parents = cell(1, (order + 1)*ss);
|
wolffd@0
|
156 for i=1:(order + 1)*ss
|
wolffd@0
|
157 bnet.parents{i} = parents(bnet.dag, i);
|
wolffd@0
|
158 end
|
wolffd@0
|
159
|
wolffd@0
|
160 bnet.auto_regressive = zeros(1,ss);
|
wolffd@0
|
161 % ar(i)=1 means (observed) node i depends on i in the previous slice
|
wolffd@0
|
162 for o=bnet.observed(:)'
|
wolffd@0
|
163 if any(bnet.parents{o+ss} <= ss)
|
wolffd@0
|
164 bnet.auto_regressive(o) = 1;
|
wolffd@0
|
165 end
|
wolffd@0
|
166 end
|
wolffd@0
|
167
|
wolffd@0
|
168
|
wolffd@0
|
169
|
wolffd@0
|
170
|
wolffd@0
|
171
|
wolffd@0
|
172
|
wolffd@0
|
173
|
wolffd@0
|
174
|
wolffd@0
|
175
|
wolffd@0
|
176
|
wolffd@0
|
177
|
wolffd@0
|
178
|
wolffd@0
|
179
|
wolffd@0
|
180
|
wolffd@0
|
181
|