25 using namespace shogun;
39 for (int32_t i=0; i<n; i++)
51 ASSERT(int32_t(components.size())==coefficients.
vlen)
57 for (int32_t i=0; i<int32_t(components.size()); i++)
67 for (int32_t i=0; i<int32_t(components.size()); i++)
71 m_components[i]->set_cov_type(components[i]->get_cov_type());
83 if (components[i]->get_cov_type()==
FULL)
106 for (int32_t i = 0; i < int32_t(
m_components.size()); i++)
121 SG_ERROR(
"Specified features are not of type CDotFeatures\n")
131 SG_ERROR(
"No features to train on.\n")
141 init_k_means->
train(dotdata);
144 alpha=alpha_init(init_means);
160 while (iter<max_iter)
162 log_likelihood_prev=log_likelihood_cur;
163 log_likelihood_cur=0;
165 for (int32_t i=0; i<num_vectors; i++)
176 log_likelihood_cur+=logPx[i];
185 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
196 return log_likelihood_cur;
202 SG_ERROR(
"No features to train on.\n")
205 SG_ERROR(
"Can't run SMEM with less than 3 component mixture model.\n")
221 int32_t* split_ind=SG_MALLOC(int32_t,
m_components.size());
224 while (iter<max_iter)
229 for (int32_t i=0; i<num_vectors; i++)
251 for (int32_t k=j+1; k<int32_t(
m_components.size()); k++)
265 for (int32_t j=0; j<num_vectors; j++)
270 for (int32_t j=i+1; j<int32_t(
m_components.size()); j++)
280 bool better_found=
false;
281 int32_t candidates_checked=0;
286 if (merge_ind[j]/int32_t(
m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%
m_components.size()) != split_ind[i])
288 candidates_checked++;
291 candidate->partial_em(split_ind[i], merge_ind[j]/int32_t(
m_components.size()), merge_ind[j]%int32_t(
m_components.size()), min_cov, max_em_iter, min_change);
294 if (cand_likelihood>cur_likelihood)
296 cur_likelihood=cand_likelihood;
300 for (int32_t k=0; k<int32_t(candidate->
get_comp().size()); k++)
311 if (candidates_checked>=max_cand)
315 if (better_found || candidates_checked>=max_cand)
329 SG_FREE(logPostSum2);
330 SG_FREE(logPostSumSum);
334 return cur_likelihood;
337 void CGMM::partial_em(int32_t comp1, int32_t comp2, int32_t comp3,
float64_t min_cov, int32_t max_em_iter,
float64_t min_change)
347 for (int32_t i=0; i<num_vectors; i++)
357 if (j!=comp1 && j!=comp2 && j!=comp3)
369 vector<CGaussian*> components(3);
377 float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
379 int32_t dim_n=components[0]->get_mean().vlen;
381 float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
382 float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
388 components[2]->get_mean().vector, dim_n);
390 for (int32_t i=0; i<dim_n; i++)
392 components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+
CMath::randn_double()*noise_mag;
393 components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+
CMath::randn_double()*noise_mag;
396 coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
397 coefficients.vector[2]=coefficients.vector[0]*0.5;
398 coefficients.vector[0]=coefficients.vector[0]*0.5;
400 if (components[0]->get_cov_type()==
FULL)
407 components[1]->set_u(c1);
410 for (int32_t i=0; i<dim_n; i++)
412 new_d+=
CMath::log(components[0]->get_d().vector[i]);
413 for (int32_t j=0; j<dim_n; j++)
417 components[0]->get_u().matrix[i*dim_n+j]=1;
418 components[2]->get_u().matrix[i*dim_n+j]=1;
422 components[0]->get_u().matrix[i*dim_n+j]=0;
423 components[2]->get_u().matrix[i*dim_n+j]=0;
428 for (int32_t i=0; i<dim_n; i++)
430 components[0]->get_d().vector[i]=new_d;
431 components[2]->get_d().vector[i]=new_d;
434 else if(components[0]->get_cov_type()==
DIAG)
437 alpha2, components[2]->get_d().vector, dim_n);
440 for (int32_t i=0; i<dim_n; i++)
442 new_d+=
CMath::log(components[0]->get_d().vector[i]);
445 for (int32_t i=0; i<dim_n; i++)
447 components[0]->get_d().vector[i]=new_d;
448 components[2]->get_d().vector[i]=new_d;
451 else if(components[0]->get_cov_type()==
SPHERICAL)
453 components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+
454 alpha2*components[2]->get_d().vector[0];
456 components[2]->get_d().vector[0]=components[0]->get_d().vector[0];
459 CGMM* partial_candidate=
new CGMM(components, coefficients);
470 while (iter<max_em_iter)
472 log_likelihood_prev=log_likelihood_cur;
473 log_likelihood_cur=0;
475 for (int32_t i=0; i<num_vectors; i++)
479 for (int32_t j=0; j<3; j++)
481 logPxy[i*3+j]=components[j]->compute_log_PDF(v)+
CMath::log(coefficients[j]);
485 logPx[i]=
CMath::log(logPx[i]+init_logPx_fix[i]);
486 log_likelihood_cur+=logPx[i];
488 for (int32_t j=0; j<3; j++)
491 alpha.matrix[i*3+j]=
CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
495 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
509 delete partial_candidate;
512 SG_FREE(init_logPxy);
514 SG_FREE(init_logPx_fix);
528 for (int32_t i=0; i<alpha.
num_cols; i++)
532 memset(mean_sum, 0, num_dim*
sizeof(
float64_t));
534 for (int32_t j=0; j<alpha.
num_rows; j++)
541 for (int32_t j=0; j<num_dim; j++)
542 mean_sum[j]/=alpha_sum;
550 cov_sum=SG_MALLOC(
float64_t, num_dim*num_dim);
551 memset(cov_sum, 0, num_dim*num_dim*
sizeof(
float64_t));
553 else if(cov_type==
DIAG)
556 memset(cov_sum, 0, num_dim*
sizeof(
float64_t));
564 for (int32_t j=0; j<alpha.
num_rows; j++)
573 1, (
double*) cov_sum, num_dim);
577 for (int32_t k=0; k<num_dim; k++)
584 for (int32_t k=0; k<num_dim; k++)
595 for (int32_t j=0; j<num_dim*num_dim; j++)
596 cov_sum[j]/=alpha_sum;
600 for (int32_t j=0; j<num_dim; j++)
608 for (int32_t j=0; j<num_dim; j++)
610 cov_sum[j]/=alpha_sum;
618 cov_sum[0]/=alpha_sum*num_dim;
627 alpha_sum_sum+=alpha_sum;
630 for (int32_t i=0; i<alpha.
num_cols; i++)
724 for (int32_t i=0; i<init_means.
num_cols; i++)
725 label_num.vector[i]=i;
734 for (int32_t i=0; i<num_vectors; i++)
750 if (cum_sum>=rand_num)
771 void CGMM::register_params()