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