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
|
function [LL, prior, transmat, mu, Sigma, mixmat] = ...
mhmm_em(data, prior, transmat, mu, Sigma, mixmat, varargin);
% LEARN_MHMM Compute the ML parameters of an HMM with (mixtures of) Gaussians output using EM.
% [ll_trace, prior, transmat, mu, sigma, mixmat] = learn_mhmm(data, ...
% prior0, transmat0, mu0, sigma0, mixmat0, ...)
%
% Notation: Q(t) = hidden state, Y(t) = observation, M(t) = mixture variable
%
% INPUTS:
% data{ex}(:,t) or data(:,t,ex) if all sequences have the same length
% prior(i) = Pr(Q(1) = i),
% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
% mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k ]
% Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k]
% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to [] or ones(Q,1) if only one mixture component
%
% Optional parameters may be passed as 'param_name', param_value pairs.
% Parameter names are shown below; default values in [] - if none, argument is mandatory.
%
% 'max_iter' - max number of EM iterations [10]
% 'thresh' - convergence threshold [1e-4]
% 'verbose' - if 1, print out loglik at every iteration [1]
% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
%
% To clamp some of the parameters, so learning does not change them:
% 'adj_prior' - if 0, do not change prior [1]
% 'adj_trans' - if 0, do not change transmat [1]
% 'adj_mix' - if 0, do not change mixmat [1]
% 'adj_mu' - if 0, do not change mu [1]
% 'adj_Sigma' - if 0, do not change Sigma [1]
%
% If the number of mixture components differs depending on Q, just set the trailing
% entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2,
% then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0.
if ~isstr(varargin{1}) % catch old syntax
error('optional arguments should be passed as string/value pairs')
end
[max_iter, thresh, verbose, cov_type, adj_prior, adj_trans, adj_mix, adj_mu, adj_Sigma] = ...
process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ...
'cov_type', 'full', 'adj_prior', 1, 'adj_trans', 1, 'adj_mix', 1, ...
'adj_mu', 1, 'adj_Sigma', 1);
previous_loglik = -inf;
loglik = 0;
converged = 0;
num_iter = 1;
LL = [];
if ~iscell(data)
data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell
end
numex = length(data);
O = size(data{1},1);
Q = length(prior);
if isempty(mixmat)
mixmat = ones(Q,1);
end
M = size(mixmat,2);
if M == 1
adj_mix = 0;
end
while (num_iter <= max_iter) & ~converged
% E step
[loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
ess_mhmm(prior, transmat, mixmat, mu, Sigma, data);
% M step
if adj_prior
prior = normalise(exp_num_visits1);
end
if adj_trans
transmat = mk_stochastic(exp_num_trans);
end
if adj_mix
mixmat = mk_stochastic(postmix);
end
if adj_mu | adj_Sigma
[mu2, Sigma2] = mixgauss_Mstep(postmix, m, op, ip, 'cov_type', cov_type);
if adj_mu
mu = reshape(mu2, [O Q M]);
end
if adj_Sigma
Sigma = reshape(Sigma2, [O O Q M]);
end
end
if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
num_iter = num_iter + 1;
converged = em_converged(loglik, previous_loglik, thresh);
previous_loglik = loglik;
LL = [LL loglik];
end
%%%%%%%%%
function [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
ess_mhmm(prior, transmat, mixmat, mu, Sigma, data)
% ESS_MHMM Compute the Expected Sufficient Statistics for a MOG Hidden Markov Model.
%
% Outputs:
% exp_num_trans(i,j) = sum_l sum_{t=2}^T Pr(Q(t-1) = i, Q(t) = j| Obs(l))
% exp_num_visits1(i) = sum_l Pr(Q(1)=i | Obs(l))
%
% Let w(i,k,t,l) = P(Q(t)=i, M(t)=k | Obs(l))
% where Obs(l) = Obs(:,:,l) = O_1 .. O_T for sequence l
% Then
% postmix(i,k) = sum_l sum_t w(i,k,t,l) (posterior mixing weights/ responsibilities)
% m(:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)
% ip(i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)' * Obs(:,t,l)
% op(:,:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l) * Obs(:,t,l)'
verbose = 0;
%[O T numex] = size(data);
numex = length(data);
O = size(data{1},1);
Q = length(prior);
M = size(mixmat,2);
exp_num_trans = zeros(Q,Q);
exp_num_visits1 = zeros(Q,1);
postmix = zeros(Q,M);
m = zeros(O,Q,M);
op = zeros(O,O,Q,M);
ip = zeros(Q,M);
mix = (M>1);
loglik = 0;
if verbose, fprintf(1, 'forwards-backwards example # '); end
for ex=1:numex
if verbose, fprintf(1, '%d ', ex); end
%obs = data(:,:,ex);
obs = data{ex};
T = size(obs,2);
if mix
[B, B2] = mixgauss_prob(obs, mu, Sigma, mixmat);
[alpha, beta, gamma, current_loglik, xi, gamma2] = ...
fwdback(prior, transmat, B, 'obslik2', B2, 'mixmat', mixmat);
else
B = mixgauss_prob(obs, mu, Sigma);
[alpha, beta, gamma, current_loglik, xi] = fwdback(prior, transmat, B);
end
loglik = loglik + current_loglik;
if verbose, fprintf(1, 'll at ex %d = %f\n', ex, loglik); end
exp_num_trans = exp_num_trans + sum(xi,3);
exp_num_visits1 = exp_num_visits1 + gamma(:,1);
if mix
postmix = postmix + sum(gamma2,3);
else
postmix = postmix + sum(gamma,2);
gamma2 = reshape(gamma, [Q 1 T]); % gamma2(i,m,t) = gamma(i,t)
end
for i=1:Q
for k=1:M
w = reshape(gamma2(i,k,:), [1 T]); % w(t) = w(i,k,t,l)
wobs = obs .* repmat(w, [O 1]); % wobs(:,t) = w(t) * obs(:,t)
m(:,i,k) = m(:,i,k) + sum(wobs, 2); % m(:) = sum_t w(t) obs(:,t)
op(:,:,i,k) = op(:,:,i,k) + wobs * obs'; % op(:,:) = sum_t w(t) * obs(:,t) * obs(:,t)'
ip(i,k) = ip(i,k) + sum(sum(wobs .* obs, 2)); % ip = sum_t w(t) * obs(:,t)' * obs(:,t)
end
end
end
if verbose, fprintf(1, '\n'); end
|
