wolffd@0
|
1 function CPD = maximize_params(CPD, temp)
|
wolffd@0
|
2 % MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
|
wolffd@0
|
3 % CPD = maximize_params(CPD, temperature)
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % Temperature is currently ignored.
|
wolffd@0
|
6
|
wolffd@0
|
7 if ~adjustable_CPD(CPD), return; end
|
wolffd@0
|
8
|
wolffd@0
|
9 CPD1 = struct(new_maximize_params(CPD));
|
wolffd@0
|
10 CPD2 = struct(old_maximize_params(CPD));
|
wolffd@0
|
11 assert(approxeq(CPD1.mean, CPD2.mean))
|
wolffd@0
|
12 assert(approxeq(CPD1.cov, CPD2.cov))
|
wolffd@0
|
13 assert(approxeq(CPD1.weights, CPD2.weights))
|
wolffd@0
|
14
|
wolffd@0
|
15 CPD = new_maximize_params(CPD);
|
wolffd@0
|
16
|
wolffd@0
|
17 %%%%%%%
|
wolffd@0
|
18 function CPD = new_maximize_params(CPD)
|
wolffd@0
|
19
|
wolffd@0
|
20 if CPD.clamped_mean
|
wolffd@0
|
21 cl_mean = CPD.mean;
|
wolffd@0
|
22 else
|
wolffd@0
|
23 cl_mean = [];
|
wolffd@0
|
24 end
|
wolffd@0
|
25
|
wolffd@0
|
26 if CPD.clamped_cov
|
wolffd@0
|
27 cl_cov = CPD.cov;
|
wolffd@0
|
28 else
|
wolffd@0
|
29 cl_cov = [];
|
wolffd@0
|
30 end
|
wolffd@0
|
31
|
wolffd@0
|
32 if CPD.clamped_weights
|
wolffd@0
|
33 cl_weights = CPD.weights;
|
wolffd@0
|
34 else
|
wolffd@0
|
35 cl_weights = [];
|
wolffd@0
|
36 end
|
wolffd@0
|
37
|
wolffd@0
|
38 [ssz psz Q] = size(CPD.weights);
|
wolffd@0
|
39
|
wolffd@0
|
40 prior = repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]);
|
wolffd@0
|
41 [CPD.mean, CPD.cov, CPD.weights] = ...
|
wolffd@0
|
42 Mstep_clg('w', CPD.Wsum, 'YY', CPD.WYYsum, 'Y', CPD.WYsum, 'YTY', [], ...
|
wolffd@0
|
43 'XX', CPD.WXXsum, 'XY', CPD.WXYsum, 'X', CPD.WXsum, ...
|
wolffd@0
|
44 'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ...
|
wolffd@0
|
45 'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ...
|
wolffd@0
|
46 'tied_cov', CPD.tied_cov, ...
|
wolffd@0
|
47 'cov_prior', prior);
|
wolffd@0
|
48
|
wolffd@0
|
49
|
wolffd@0
|
50 %%%%%%%%%%%
|
wolffd@0
|
51
|
wolffd@0
|
52 function CPD = old_maximize_params(CPD)
|
wolffd@0
|
53
|
wolffd@0
|
54
|
wolffd@0
|
55 if ~adjustable_CPD(CPD), return; end
|
wolffd@0
|
56
|
wolffd@0
|
57 %assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));
|
wolffd@0
|
58 assert(~any(isnan(CPD.WXXsum)))
|
wolffd@0
|
59 assert(~any(isnan(CPD.WXYsum)))
|
wolffd@0
|
60 assert(~any(isnan(CPD.WYYsum)))
|
wolffd@0
|
61
|
wolffd@0
|
62 [self_size cpsize dpsize] = size(CPD.weights);
|
wolffd@0
|
63
|
wolffd@0
|
64 % Append 1s to the parents, and derive the corresponding cross products.
|
wolffd@0
|
65 % This is used when estimate the means and weights simultaneosuly,
|
wolffd@0
|
66 % and when estimatting Sigma.
|
wolffd@0
|
67 % Let x2 = [x 1]'
|
wolffd@0
|
68 XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)'
|
wolffd@0
|
69 XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)'
|
wolffd@0
|
70 YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)'
|
wolffd@0
|
71 for i=1:dpsize
|
wolffd@0
|
72 XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
|
wolffd@0
|
73 CPD.WYsum(:,i)']; % 1*Y
|
wolffd@0
|
74 % [x * [x' 1] = [xx' x
|
wolffd@0
|
75 % 1] x' 1]
|
wolffd@0
|
76 XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
|
wolffd@0
|
77 CPD.WXsum(:,i)' CPD.Wsum(i)];
|
wolffd@0
|
78 YY(:,:,i) = CPD.WYYsum(:,:,i);
|
wolffd@0
|
79 end
|
wolffd@0
|
80
|
wolffd@0
|
81 w = CPD.Wsum(:);
|
wolffd@0
|
82 % Set any zeros to one before dividing
|
wolffd@0
|
83 % This is valid because w(i)=0 => WYsum(:,i)=0, etc
|
wolffd@0
|
84 w = w + (w==0);
|
wolffd@0
|
85
|
wolffd@0
|
86 if CPD.clamped_mean
|
wolffd@0
|
87 % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent
|
wolffd@0
|
88 % to estimating B and then adding the clamped_mean to the last column.
|
wolffd@0
|
89 if ~CPD.clamped_weights
|
wolffd@0
|
90 B = zeros(self_size, cpsize, dpsize);
|
wolffd@0
|
91 for i=1:dpsize
|
wolffd@0
|
92 if det(CPD.WXXsum(:,:,i))==0
|
wolffd@0
|
93 B(:,:,i) = 0;
|
wolffd@0
|
94 else
|
wolffd@0
|
95 % Eqn 9 in table 2 of TR
|
wolffd@0
|
96 %B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
|
wolffd@0
|
97 B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
|
wolffd@0
|
98 end
|
wolffd@0
|
99 end
|
wolffd@0
|
100 %CPD.weights = reshape(B, [self_size cpsize dpsize]);
|
wolffd@0
|
101 CPD.weights = B;
|
wolffd@0
|
102 end
|
wolffd@0
|
103 elseif CPD.clamped_weights % KPM 1/25/02
|
wolffd@0
|
104 if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals
|
wolffd@0
|
105 for i=1:dpsize
|
wolffd@0
|
106 CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i);
|
wolffd@0
|
107 end
|
wolffd@0
|
108 end
|
wolffd@0
|
109 else % nothing is clamped, so estimate mean and weights simultaneously
|
wolffd@0
|
110 B2 = zeros(self_size, cpsize+1, dpsize);
|
wolffd@0
|
111 for i=1:dpsize
|
wolffd@0
|
112 if det(XX(:,:,i))==0 % fix by U. Sondhauss 6/27/99
|
wolffd@0
|
113 B2(:,:,i)=0;
|
wolffd@0
|
114 else
|
wolffd@0
|
115 % Eqn 9 in table 2 of TR
|
wolffd@0
|
116 %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
|
wolffd@0
|
117 B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
|
wolffd@0
|
118 end
|
wolffd@0
|
119 CPD.mean(:,i) = B2(:,cpsize+1,i);
|
wolffd@0
|
120 CPD.weights(:,:,i) = B2(:,1:cpsize,i);
|
wolffd@0
|
121 end
|
wolffd@0
|
122 end
|
wolffd@0
|
123
|
wolffd@0
|
124 % Let B2 = [W mu]
|
wolffd@0
|
125 if cpsize>0
|
wolffd@0
|
126 B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
|
wolffd@0
|
127 end
|
wolffd@0
|
128 B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
|
wolffd@0
|
129
|
wolffd@0
|
130 % To avoid singular covariance matrices,
|
wolffd@0
|
131 % we use the regularization method suggested in "A Quasi-Bayesian approach to estimating
|
wolffd@0
|
132 % parameters for mixtures of normal distributions", Hamilton 91.
|
wolffd@0
|
133 % If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma),
|
wolffd@0
|
134 % where gamma >=0 is a smoothing parameter (equivalent sample size of I prior)
|
wolffd@0
|
135
|
wolffd@0
|
136 gamma = CPD.cov_prior_weight;
|
wolffd@0
|
137
|
wolffd@0
|
138 if ~CPD.clamped_cov
|
wolffd@0
|
139 if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99
|
wolffd@0
|
140 Z = 1-temp;
|
wolffd@0
|
141 % When temp > 1, Z is negative, so we are dividing by a smaller
|
wolffd@0
|
142 % number, ie. increasing the variance.
|
wolffd@0
|
143 else
|
wolffd@0
|
144 Z = 0;
|
wolffd@0
|
145 end
|
wolffd@0
|
146 if CPD.tied_cov
|
wolffd@0
|
147 S = zeros(self_size, self_size);
|
wolffd@0
|
148 % Eqn 2 from table 2 in TR
|
wolffd@0
|
149 for i=1:dpsize
|
wolffd@0
|
150 S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
|
wolffd@0
|
151 end
|
wolffd@0
|
152 %denom = CPD.nsamples + gamma + Z;
|
wolffd@0
|
153 denom = CPD.nsamples + Z;
|
wolffd@0
|
154 S = (S + gamma*eye(self_size)) / denom;
|
wolffd@0
|
155 if strcmp(CPD.cov_type, 'diag')
|
wolffd@0
|
156 S = diag(diag(S));
|
wolffd@0
|
157 end
|
wolffd@0
|
158 CPD.cov = repmat(S, [1 1 dpsize]);
|
wolffd@0
|
159 else
|
wolffd@0
|
160 for i=1:dpsize
|
wolffd@0
|
161 % Eqn 1 from table 2 in TR
|
wolffd@0
|
162 S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i);
|
wolffd@0
|
163 %denom = w(i) + gamma + Z;
|
wolffd@0
|
164 denom = w(i) + Z;
|
wolffd@0
|
165 S = (S + gamma*eye(self_size)) / denom;
|
wolffd@0
|
166 CPD.cov(:,:,i) = S;
|
wolffd@0
|
167 end
|
wolffd@0
|
168 if strcmp(CPD.cov_type, 'diag')
|
wolffd@0
|
169 for i=1:dpsize
|
wolffd@0
|
170 CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
|
wolffd@0
|
171 end
|
wolffd@0
|
172 end
|
wolffd@0
|
173 end
|
wolffd@0
|
174 end
|
wolffd@0
|
175
|
wolffd@0
|
176
|
wolffd@0
|
177 check_covars = 0;
|
wolffd@0
|
178 min_covar = 1e-5;
|
wolffd@0
|
179 if check_covars % prevent collapsing to a point
|
wolffd@0
|
180 for i=1:dpsize
|
wolffd@0
|
181 if min(svd(CPD.cov(:,:,i))) < min_covar
|
wolffd@0
|
182 disp(['resetting singular covariance for node ' num2str(CPD.self)]);
|
wolffd@0
|
183 CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
|
wolffd@0
|
184 end
|
wolffd@0
|
185 end
|
wolffd@0
|
186 end
|
wolffd@0
|
187
|
wolffd@0
|
188
|
wolffd@0
|
189
|