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
180
181
|
function [mix, options, errlog] = gmmem(mix, x, options)
%GMMEM EM algorithm for Gaussian mixture model.
%
% Description
% [MIX, OPTIONS, ERRLOG] = GMMEM(MIX, X, OPTIONS) uses the Expectation
% Maximization algorithm of Dempster et al. to estimate the parameters
% of a Gaussian mixture model defined by a data structure MIX. The
% matrix X represents the data whose expectation is maximized, with
% each row corresponding to a vector. The optional parameters have
% the following interpretations.
%
% OPTIONS(1) is set to 1 to display error values; also logs error
% values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then
% only warning messages are displayed. If OPTIONS(1) is -1, then
% nothing is displayed.
%
% OPTIONS(3) is a measure of the absolute precision required of the
% error function at the solution. If the change in log likelihood
% between two steps of the EM algorithm is less than this value, then
% the function terminates.
%
% OPTIONS(5) is set to 1 if a covariance matrix is reset to its
% original value when any of its singular values are too small (less
% than MIN_COVAR which has the value eps). With the default value of
% 0 no action is taken.
%
% OPTIONS(14) is the maximum number of iterations; default 100.
%
% The optional return value OPTIONS contains the final error value
% (i.e. data log likelihood) in OPTIONS(8).
%
% See also
% GMM, GMMINIT
%
% Copyright (c) Ian T Nabney (1996-2001)
% Check that inputs are consistent
errstring = consist(mix, 'gmm', x);
if ~isempty(errstring)
error(errstring);
end
[ndata, xdim] = size(x);
% Sort out the options
if (options(14))
niters = options(14);
else
niters = 100;
end
display = options(1);
store = 0;
if (nargout > 2)
store = 1; % Store the error values to return them
errlog = zeros(1, niters);
end
test = 0;
if options(3) > 0.0
test = 1; % Test log likelihood for termination
end
check_covars = 0;
if options(5) >= 1
if display >= 0
disp('check_covars is on');
end
check_covars = 1; % Ensure that covariances don't collapse
MIN_COVAR = eps; % Minimum singular value of covariance matrix
init_covars = mix.covars;
end
% Main loop of algorithm
for n = 1:niters
% Calculate posteriors based on old parameters
[post, act] = gmmpost(mix, x);
% Calculate error value if needed
if (display || store || test)
prob = act*(mix.priors)';
% Error value is negative log likelihood of data
e = - sum(log(prob));
if store
errlog(n) = e;
end
if display > 0
fprintf(1, 'Cycle %4d Error %11.6f\n', n, e);
end
if test
if (n > 1 && abs(e - eold) < options(3))
options(8) = e;
return;
else
eold = e;
end
end
end
% Adjust the new estimates for the parameters
new_pr = sum(post, 1);
new_c = post' * x;
% Now move new estimates to old parameter vectors
mix.priors = new_pr ./ ndata;
mix.centres = new_c ./ (new_pr' * ones(1, mix.nin));
switch mix.covar_type
case 'spherical'
n2 = dist2(x, mix.centres);
for j = 1:mix.ncentres
v(j) = (post(:,j)'*n2(:,j));
end
mix.covars = ((v./new_pr))./mix.nin;
if check_covars
% Ensure that no covariance is too small
for j = 1:mix.ncentres
if mix.covars(j) < MIN_COVAR
mix.covars(j) = init_covars(j);
end
end
end
case 'diag'
for j = 1:mix.ncentres
diffs = x - (ones(ndata, 1) * mix.centres(j,:));
mix.covars(j,:) = sum((diffs.*diffs).*(post(:,j)*ones(1, ...
mix.nin)), 1)./new_pr(j);
end
if check_covars
% Ensure that no covariance is too small
for j = 1:mix.ncentres
if min(mix.covars(j,:)) < MIN_COVAR
mix.covars(j,:) = init_covars(j,:);
end
end
end
case 'full'
for j = 1:mix.ncentres
diffs = x - (ones(ndata, 1) * mix.centres(j,:));
diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
mix.covars(:,:,j) = (diffs'*diffs)/new_pr(j);
end
if check_covars
% Ensure that no covariance is too small
for j = 1:mix.ncentres
if min(svd(mix.covars(:,:,j))) < MIN_COVAR
mix.covars(:,:,j) = init_covars(:,:,j);
end
end
end
case 'ppca'
for j = 1:mix.ncentres
diffs = x - (ones(ndata, 1) * mix.centres(j,:));
diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
[tempcovars, tempU, templambda] = ...
ppca((diffs'*diffs)/new_pr(j), mix.ppca_dim);
if length(templambda) ~= mix.ppca_dim
error('Unable to extract enough components');
else
mix.covars(j) = tempcovars;
mix.U(:, :, j) = tempU;
mix.lambda(j, :) = templambda;
end
end
if check_covars
if mix.covars(j) < MIN_COVAR
mix.covars(j) = init_covars(j);
end
end
otherwise
error(['Unknown covariance type ', mix.covar_type]);
end
end
options(8) = -sum(log(gmmprob(mix, x)));
if (display >= 0)
disp(maxitmess);
end
|
