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")
142 init_k_means->
train(dotdata);
145 alpha=alpha_init(init_means);
161 while (iter<max_iter)
163 log_likelihood_prev=log_likelihood_cur;
164 log_likelihood_cur=0;
166 for (int32_t i=0; i<num_vectors; i++)
177 log_likelihood_cur+=logPx[i];
186 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
197 return log_likelihood_cur;
203 SG_ERROR(
"No features to train on.\n")
206 SG_ERROR(
"Can't run SMEM with less than 3 component mixture model.\n")
222 int32_t* split_ind=SG_MALLOC(int32_t,
m_components.size());
225 while (iter<max_iter)
230 for (int32_t i=0; i<num_vectors; i++)
252 for (int32_t k=j+1; k<int32_t(
m_components.size()); k++)
266 for (int32_t j=0; j<num_vectors; j++)
271 for (int32_t j=i+1; j<int32_t(
m_components.size()); j++)
281 bool better_found=
false;
282 int32_t candidates_checked=0;
287 if (merge_ind[j]/int32_t(
m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%
m_components.size()) != split_ind[i])
289 candidates_checked++;
292 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);
295 if (cand_likelihood>cur_likelihood)
297 cur_likelihood=cand_likelihood;
301 for (int32_t k=0; k<int32_t(candidate->
get_comp().size()); k++)
312 if (candidates_checked>=max_cand)
316 if (better_found || candidates_checked>=max_cand)
330 SG_FREE(logPostSum2);
331 SG_FREE(logPostSumSum);
335 return cur_likelihood;
338 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)
348 for (int32_t i=0; i<num_vectors; i++)
358 if (j!=comp1 && j!=comp2 && j!=comp3)
370 vector<CGaussian*> components(3);
378 float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
380 int32_t dim_n=components[0]->get_mean().vlen;
382 float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
383 float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
389 components[2]->get_mean().vector, dim_n);
391 for (int32_t i=0; i<dim_n; i++)
393 components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+
CMath::randn_double()*noise_mag;
394 components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+
CMath::randn_double()*noise_mag;
397 coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
398 coefficients.vector[2]=coefficients.vector[0]*0.5;
399 coefficients.vector[0]=coefficients.vector[0]*0.5;
401 if (components[0]->get_cov_type()==
FULL)
408 components[1]->set_u(c1);
411 for (int32_t i=0; i<dim_n; i++)
413 new_d+=
CMath::log(components[0]->get_d().vector[i]);
414 for (int32_t j=0; j<dim_n; j++)
418 components[0]->get_u().matrix[i*dim_n+j]=1;
419 components[2]->get_u().matrix[i*dim_n+j]=1;
423 components[0]->get_u().matrix[i*dim_n+j]=0;
424 components[2]->get_u().matrix[i*dim_n+j]=0;
429 for (int32_t i=0; i<dim_n; i++)
431 components[0]->get_d().vector[i]=new_d;
432 components[2]->get_d().vector[i]=new_d;
435 else if(components[0]->get_cov_type()==
DIAG)
438 alpha2, components[2]->get_d().vector, dim_n);
441 for (int32_t i=0; i<dim_n; i++)
443 new_d+=
CMath::log(components[0]->get_d().vector[i]);
446 for (int32_t i=0; i<dim_n; i++)
448 components[0]->get_d().vector[i]=new_d;
449 components[2]->get_d().vector[i]=new_d;
452 else if(components[0]->get_cov_type()==
SPHERICAL)
454 components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+
455 alpha2*components[2]->get_d().vector[0];
457 components[2]->get_d().vector[0]=components[0]->get_d().vector[0];
460 CGMM* partial_candidate=
new CGMM(components, coefficients);
471 while (iter<max_em_iter)
473 log_likelihood_prev=log_likelihood_cur;
474 log_likelihood_cur=0;
476 for (int32_t i=0; i<num_vectors; i++)
480 for (int32_t j=0; j<3; j++)
482 logPxy[i*3+j]=components[j]->compute_log_PDF(v)+
CMath::log(coefficients[j]);
486 logPx[i]=
CMath::log(logPx[i]+init_logPx_fix[i]);
487 log_likelihood_cur+=logPx[i];
489 for (int32_t j=0; j<3; j++)
492 alpha.matrix[i*3+j]=
CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
496 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
510 delete partial_candidate;
513 SG_FREE(init_logPxy);
515 SG_FREE(init_logPx_fix);
529 for (int32_t i=0; i<alpha.
num_cols; i++)
533 memset(mean_sum, 0, num_dim*
sizeof(
float64_t));
535 for (int32_t j=0; j<alpha.
num_rows; j++)
542 for (int32_t j=0; j<num_dim; j++)
543 mean_sum[j]/=alpha_sum;
551 cov_sum=SG_MALLOC(
float64_t, num_dim*num_dim);
552 memset(cov_sum, 0, num_dim*num_dim*
sizeof(
float64_t));
554 else if(cov_type==
DIAG)
557 memset(cov_sum, 0, num_dim*
sizeof(
float64_t));
565 for (int32_t j=0; j<alpha.
num_rows; j++)
574 1, (
double*) cov_sum, num_dim);
578 for (int32_t k=0; k<num_dim; k++)
585 for (int32_t k=0; k<num_dim; k++)
596 for (int32_t j=0; j<num_dim*num_dim; j++)
597 cov_sum[j]/=alpha_sum;
601 for (int32_t j=0; j<num_dim; j++)
609 for (int32_t j=0; j<num_dim; j++)
611 cov_sum[j]/=alpha_sum;
619 cov_sum[0]/=alpha_sum*num_dim;
628 alpha_sum_sum+=alpha_sum;
631 for (int32_t i=0; i<alpha.
num_cols; i++)
747 for (int32_t i=0; i<init_means.
num_cols; i++)
748 label_num.vector[i]=i;
757 for (int32_t i=0; i<num_vectors; i++)
774 if (cum_sum>=rand_num)
776 SG_DEBUG(
"Sampling from mixture component %d\n", i);
798 void CGMM::register_params()