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
|
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);
|
