1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
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);
|
