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
182
183
184
185
186
187
188
189
|
function CPD = maximize_params(CPD, temp)
% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
% CPD = maximize_params(CPD, temperature)
%
% Temperature is currently ignored.
if ~adjustable_CPD(CPD), return; end
CPD1 = struct(new_maximize_params(CPD));
CPD2 = struct(old_maximize_params(CPD));
assert(approxeq(CPD1.mean, CPD2.mean))
assert(approxeq(CPD1.cov, CPD2.cov))
assert(approxeq(CPD1.weights, CPD2.weights))
CPD = new_maximize_params(CPD);
%%%%%%%
function CPD = new_maximize_params(CPD)
if CPD.clamped_mean
cl_mean = CPD.mean;
else
cl_mean = [];
end
if CPD.clamped_cov
cl_cov = CPD.cov;
else
cl_cov = [];
end
if CPD.clamped_weights
cl_weights = CPD.weights;
else
cl_weights = [];
end
[ssz psz Q] = size(CPD.weights);
prior = repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]);
[CPD.mean, CPD.cov, CPD.weights] = ...
Mstep_clg('w', CPD.Wsum, 'YY', CPD.WYYsum, 'Y', CPD.WYsum, 'YTY', [], ...
'XX', CPD.WXXsum, 'XY', CPD.WXYsum, 'X', CPD.WXsum, ...
'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ...
'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ...
'tied_cov', CPD.tied_cov, ...
'cov_prior', prior);
%%%%%%%%%%%
function CPD = old_maximize_params(CPD)
if ~adjustable_CPD(CPD), return; end
%assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));
assert(~any(isnan(CPD.WXXsum)))
assert(~any(isnan(CPD.WXYsum)))
assert(~any(isnan(CPD.WYYsum)))
[self_size cpsize dpsize] = size(CPD.weights);
% Append 1s to the parents, and derive the corresponding cross products.
% This is used when estimate the means and weights simultaneosuly,
% and when estimatting Sigma.
% Let x2 = [x 1]'
XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)'
XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)'
YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)'
for i=1:dpsize
XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
CPD.WYsum(:,i)']; % 1*Y
% [x * [x' 1] = [xx' x
% 1] x' 1]
XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
CPD.WXsum(:,i)' CPD.Wsum(i)];
YY(:,:,i) = CPD.WYYsum(:,:,i);
end
w = CPD.Wsum(:);
% Set any zeros to one before dividing
% This is valid because w(i)=0 => WYsum(:,i)=0, etc
w = w + (w==0);
if CPD.clamped_mean
% Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent
% to estimating B and then adding the clamped_mean to the last column.
if ~CPD.clamped_weights
B = zeros(self_size, cpsize, dpsize);
for i=1:dpsize
if det(CPD.WXXsum(:,:,i))==0
B(:,:,i) = 0;
else
% Eqn 9 in table 2 of TR
%B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
end
end
%CPD.weights = reshape(B, [self_size cpsize dpsize]);
CPD.weights = B;
end
elseif CPD.clamped_weights % KPM 1/25/02
if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals
for i=1:dpsize
CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i);
end
end
else % nothing is clamped, so estimate mean and weights simultaneously
B2 = zeros(self_size, cpsize+1, dpsize);
for i=1:dpsize
if det(XX(:,:,i))==0 % fix by U. Sondhauss 6/27/99
B2(:,:,i)=0;
else
% Eqn 9 in table 2 of TR
%B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
end
CPD.mean(:,i) = B2(:,cpsize+1,i);
CPD.weights(:,:,i) = B2(:,1:cpsize,i);
end
end
% Let B2 = [W mu]
if cpsize>0
B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
end
B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
% To avoid singular covariance matrices,
% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating
% parameters for mixtures of normal distributions", Hamilton 91.
% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma),
% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior)
gamma = CPD.cov_prior_weight;
if ~CPD.clamped_cov
if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99
Z = 1-temp;
% When temp > 1, Z is negative, so we are dividing by a smaller
% number, ie. increasing the variance.
else
Z = 0;
end
if CPD.tied_cov
S = zeros(self_size, self_size);
% Eqn 2 from table 2 in TR
for i=1:dpsize
S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
end
%denom = CPD.nsamples + gamma + Z;
denom = CPD.nsamples + Z;
S = (S + gamma*eye(self_size)) / denom;
if strcmp(CPD.cov_type, 'diag')
S = diag(diag(S));
end
CPD.cov = repmat(S, [1 1 dpsize]);
else
for i=1:dpsize
% Eqn 1 from table 2 in TR
S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i);
%denom = w(i) + gamma + Z;
denom = w(i) + Z;
S = (S + gamma*eye(self_size)) / denom;
CPD.cov(:,:,i) = S;
end
if strcmp(CPD.cov_type, 'diag')
for i=1:dpsize
CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
end
end
end
end
check_covars = 0;
min_covar = 1e-5;
if check_covars % prevent collapsing to a point
for i=1:dpsize
if min(svd(CPD.cov(:,:,i))) < min_covar
disp(['resetting singular covariance for node ' num2str(CPD.self)]);
CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
end
end
end
|
