48 for (
int k = 0;k <
K;k++) {
53 for (
int d = 0;d <
D;d++) {
54 double tmp = c_diag_cov[d];
55 c_diag_cov_inv_etc[d] = 1.0 / (2.0 * tmp);
70 for (
int k = 0;k <
K;k++) {
77 for (
int k = 0;k <
K;k++)
78 for (
int d = 0;d <
D;d++)
87 double acc_loglhood = 0.0;
89 for (
int k = 0;k <
K;k++) {
90 c_acc_loglhood_K[k] = 0.0;
92 double * c_acc_mean = c_acc_means[k];
93 double * c_acc_cov = c_acc_covs[k];
95 for (
int d = 0;d <
D;d++) { c_acc_mean[d] = 0.0; c_acc_cov[d] = 0.0; }
98 for (
int n = 0;n <
N;n++) {
99 double * c_x =
c_X[n];
102 for (
int k = 0;k <
K;k++) {
110 double log_sum = c_tmpvecK[0];
111 for (
int k = 1;k <
K;k++) log_sum =
log_add(log_sum, c_tmpvecK[k]);
112 acc_loglhood += log_sum;
114 for (
int k = 0;k <
K;k++) {
116 double * c_acc_mean = c_acc_means[k];
117 double * c_acc_cov = c_acc_covs[k];
119 double tmp_k =
trunc_exp(c_tmpvecK[k] - log_sum);
120 acc_loglhood_K[k] += tmp_k;
122 for (
int d = 0;d <
D;d++) {
123 double tmp_x = c_x[d];
124 c_acc_mean[d] += tmp_k * tmp_x;
125 c_acc_cov[d] += tmp_k * tmp_x * tmp_x;
132 for (
int k = 0;k <
K;k++) {
double tmp =
std::exp(c_tmpvecK[k]); c_tmpvecK[k] = tmp; sum += tmp; }
135 for (
int k = 0;k <
K;k++) {
137 double * c_acc_mean = c_acc_means[k];
138 double * c_acc_cov = c_acc_covs[k];
140 double tmp_k = c_tmpvecK[k] /
sum;
141 c_acc_loglhood_K[k] += tmp_k;
143 for (
int d = 0;d <
D;d++) {
144 double tmp_x = c_x[d];
145 c_acc_mean[d] += tmp_k * tmp_x;
146 c_acc_cov[d] += tmp_k * tmp_x * tmp_x;
152 for (
int k = 0;k <
K;k++) {
157 double * c_acc_mean = c_acc_means[k];
158 double * c_acc_cov = c_acc_covs[k];
160 double tmp_k = c_acc_loglhood_K[k];
164 for (
int d = 0;d <
D;d++) {
165 double tmp_mean = c_acc_mean[d] / tmp_k;
166 c_mean[d] = tmp_mean;
167 c_diag_cov[d] = c_acc_cov[d] / tmp_k - tmp_mean * tmp_mean;
171 return(acc_loglhood / N);
182 using std::noshowpos;
183 using std::scientific;
186 using std::setprecision;
193 cout <<
"MOG_diag_EM_sup::ml_iterate()" << endl;
194 cout << setw(14) <<
"iteration";
195 cout << setw(14) <<
"avg_loglhood";
196 cout << setw(14) <<
"delta";
197 cout << setw(10) <<
"toc";
201 for (
int i = 0; i <
max_iter; i++) {
209 double delta = avg_log_lhood_new - avg_log_lhood_old;
211 cout << noshowpos << fixed;
212 cout << setw(14) << i;
213 cout << showpos << scientific << setprecision(3);
214 cout << setw(14) << avg_log_lhood_new;
215 cout << setw(14) << delta;
216 cout << noshowpos << fixed;
217 cout << setw(10) << tt.
toc();
218 cout << endl <<
flush;
221 if (avg_log_lhood_new <= avg_log_lhood_old)
break;
223 avg_log_lhood_old = avg_log_lhood_new;
231 it_assert(model_in.
is_valid(),
"MOG_diag_EM_sup::ml(): initial model not valid");
233 it_assert((max_iter_in > 0),
"MOG_diag_EM_sup::ml(): 'max_iter' needs to be greater than zero");
243 init(means_in, diag_covs_in, weights_in);
247 weights_in.set_size(0);
250 it_warning(
"MOG_diag_EM_sup::ml(): WARNING: K > N");
254 it_warning(
"MOG_diag_EM_sup::ml(): WARNING: K > N/10");
270 acc_loglhood_K.set_size(
K);
273 for (
int k = 0;k <
K;k++) acc_means(k).
set_size(
D);
275 for (
int k = 0;k <
K;k++) acc_covs(k).
set_size(
D);
298 acc_loglhood_K.set_size(0);
307 double,
double,
bool)
309 it_error(
"MOG_diag_EM_sup::map(): not implemented yet");
319 EM.
ml(model_in, X_in, max_iter_in, var_floor_in, weight_floor_in, verbose_in);
325 it_error(
"MOG_diag_MAP(): not implemented yet");