SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMM.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Alesis Novik
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 #include <shogun/lib/config.h>
11 
12 #ifdef HAVE_LAPACK
13 
14 #include <shogun/clustering/GMM.h>
17 #include <shogun/base/Parameter.h>
20 #include <shogun/features/Labels.h>
21 #include <shogun/classifier/KNN.h>
22 
23 using namespace shogun;
24 
25 CGMM::CGMM() : CDistribution(), m_components(), m_coefficients()
26 {
27  register_params();
28 }
29 
30 CGMM::CGMM(int32_t n, ECovType cov_type) : CDistribution(), m_components(), m_coefficients()
31 {
36 
37  for (int32_t i=0; i<n; i++)
38  {
39  m_components.vector[i]=new CGaussian();
41  m_components.vector[i]->set_cov_type(cov_type);
42  }
43 
44  register_params();
45 }
46 
47 CGMM::CGMM(SGVector<CGaussian*> components, SGVector<float64_t> coefficients, bool copy) : CDistribution()
48 {
49  ASSERT(components.vlen==coefficients.vlen);
50 
51  if (!copy)
52  {
53  m_components=components;
54  m_coefficients=coefficients;
55  for (int32_t i=0; i<components.vlen; i++)
56  {
58  }
59  }
60  else
61  {
63  m_coefficients.vlen=components.vlen;
65  m_components.vlen=components.vlen;
66 
67  for (int32_t i=0; i<components.vlen; i++)
68  {
69  m_components.vector[i]=new CGaussian();
71  m_components.vector[i]->set_cov_type(components.vector[i]->get_cov_type());
72 
73  SGVector<float64_t> old_mean=components.vector[i]->get_mean();
74  SGVector<float64_t> new_mean(old_mean.vlen);
75  memcpy(new_mean.vector, old_mean.vector, old_mean.vlen*sizeof(float64_t));
76  m_components.vector[i]->set_mean(new_mean);
77 
78  SGVector<float64_t> old_d=components.vector[i]->get_d();
79  SGVector<float64_t> new_d(old_d.vlen);
80  memcpy(new_d.vector, old_d.vector, old_d.vlen*sizeof(float64_t));
81  m_components.vector[i]->set_d(new_d);
82 
83  if (components.vector[i]->get_cov_type()==FULL)
84  {
85  SGMatrix<float64_t> old_u=components.vector[i]->get_u();
86  SGMatrix<float64_t> new_u(old_u.num_rows, old_u.num_cols);
87  memcpy(new_u.matrix, old_u.matrix, old_u.num_rows*old_u.num_cols*sizeof(float64_t));
88  m_components.vector[i]->set_u(new_u);
89  }
90 
91  m_coefficients.vector[i]=coefficients.vector[i];
92  }
93  }
94 
95  register_params();
96 }
97 
99 {
100  if (m_components.vector)
101  cleanup();
102 }
103 
105 {
106  for (int32_t i = 0; i < m_components.vlen; i++)
108 
111 }
112 
114 {
115  ASSERT(m_components.vlen != 0);
116 
118  if (data)
119  {
120  if (!data->has_property(FP_DOT))
121  SG_ERROR("Specified features are not of type CDotFeatures\n");
122  set_features(data);
123  }
124 
125  return true;
126 }
127 
128 float64_t CGMM::train_em(float64_t min_cov, int32_t max_iter, float64_t min_change)
129 {
130  if (!features)
131  SG_ERROR("No features to train on.\n");
132 
133  CDotFeatures* dotdata=(CDotFeatures *) features;
134  int32_t num_vectors=dotdata->get_num_vectors();
135 
136  SGMatrix<float64_t> alpha;
137 
138  if (m_components.vector[0]->get_mean().vector==NULL)
139  {
140  CKMeans* init_k_means=new CKMeans(m_components.vlen, new CEuclidianDistance());
141  init_k_means->train(dotdata);
142  SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers();
143 
144  alpha=alpha_init(init_means);
145 
146  SG_UNREF(init_k_means);
147 
148  max_likelihood(alpha, min_cov);
149  }
150  else
151  {
152  alpha.matrix=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
153  alpha.num_rows=num_vectors;
154  alpha.num_cols=m_components.vlen;
155  }
156 
157  int32_t iter=0;
158  float64_t log_likelihood_prev=0;
159  float64_t log_likelihood_cur=0;
160  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
161  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
162  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
163 
164  while (iter<max_iter)
165  {
166  log_likelihood_prev=log_likelihood_cur;
167  log_likelihood_cur=0;
168 
169  for (int32_t i=0; i<num_vectors; i++)
170  {
171  logPx[i]=0;
173  for (int32_t j=0; j<m_components.vlen; j++)
174  {
176  logPx[i]+=CMath::exp(logPxy[i*m_components.vlen+j]);
177  }
178 
179  logPx[i]=CMath::log(logPx[i]);
180  log_likelihood_cur+=logPx[i];
181  v.free_vector();
182 
183  for (int32_t j=0; j<m_components.vlen; j++)
184  {
185  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
186  alpha.matrix[i*m_components.vlen+j]=CMath::exp(logPxy[i*m_components.vlen+j]-logPx[i]);
187  }
188  }
189 
190  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
191  break;
192 
193  max_likelihood(alpha, min_cov);
194 
195  iter++;
196  }
197 
198  SG_FREE(logPxy);
199  SG_FREE(logPx);
200  //SG_FREE(logPost);
201  alpha.free_matrix();
202 
203  return log_likelihood_cur;
204 }
205 
206 float64_t CGMM::train_smem(int32_t max_iter, int32_t max_cand, float64_t min_cov, int32_t max_em_iter, float64_t min_change)
207 {
208  if (!features)
209  SG_ERROR("No features to train on.\n");
210 
211  if (m_components.vlen<3)
212  SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n");
213 
214  CDotFeatures* dotdata=(CDotFeatures *) features;
215  int32_t num_vectors=dotdata->get_num_vectors();
216 
217  float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change);
218 
219  int32_t iter=0;
220  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
221  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
222  float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
228  int32_t* split_ind=SG_MALLOC(int32_t, m_components.vlen);
229  int32_t* merge_ind=SG_MALLOC(int32_t, m_components.vlen*(m_components.vlen-1)/2);
230 
231  while (iter<max_iter)
232  {
233  memset(logPostSum, 0, m_components.vlen*sizeof(float64_t));
234  memset(logPostSum2, 0, m_components.vlen*sizeof(float64_t));
235  memset(logPostSumSum, 0, (m_components.vlen*(m_components.vlen-1)/2)*sizeof(float64_t));
236  for (int32_t i=0; i<num_vectors; i++)
237  {
238  logPx[i]=0;
240  for (int32_t j=0; j<m_components.vlen; j++)
241  {
243  logPx[i]+=CMath::exp(logPxy[i*m_components.vlen+j]);
244  }
245 
246  logPx[i]=CMath::log(logPx[i]);
247  v.free_vector();
248 
249  for (int32_t j=0; j<m_components.vlen; j++)
250  {
251  logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
252  logPostSum[j]+=CMath::exp(logPost[i*m_components.vlen+j]);
253  logPostSum2[j]+=CMath::exp(2*logPost[i*m_components.vlen+j]);
254  }
255 
256  int32_t counter=0;
257  for (int32_t j=0; j<m_components.vlen; j++)
258  {
259  for (int32_t k=j+1; k<m_components.vlen; k++)
260  {
261  logPostSumSum[counter]+=CMath::exp(logPost[i*m_components.vlen+j]+logPost[i*m_components.vlen+k]);
262  counter++;
263  }
264  }
265  }
266 
267  int32_t counter=0;
268  for (int32_t i=0; i<m_components.vlen; i++)
269  {
270  logPostSum[i]=CMath::log(logPostSum[i]);
271  split_crit[i]=0;
272  split_ind[i]=i;
273  for (int32_t j=0; j<num_vectors; j++)
274  {
275  split_crit[i]+=(logPost[j*m_components.vlen+i]-logPostSum[i]-logPxy[j*m_components.vlen+i]+CMath::log(m_coefficients.vector[i]))*
276  (CMath::exp(logPost[j*m_components.vlen+i])/CMath::exp(logPostSum[i]));
277  }
278  for (int32_t j=i+1; j<m_components.vlen; j++)
279  {
280  merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j]));
281  merge_ind[counter]=i*m_components.vlen+j;
282  counter++;
283  }
284  }
285  CMath::qsort_backward_index(split_crit, split_ind, m_components.vlen);
286  CMath::qsort_backward_index(merge_crit, merge_ind, m_components.vlen*(m_components.vlen-1)/2);
287 
288  bool better_found=false;
289  int32_t candidates_checked=0;
290  for (int32_t i=0; i<m_components.vlen; i++)
291  {
292  for (int32_t j=0; j<m_components.vlen*(m_components.vlen-1)/2; j++)
293  {
294  if (merge_ind[j]/m_components.vlen != split_ind[i] && merge_ind[j]%m_components.vlen != split_ind[i])
295  {
296  candidates_checked++;
297  CGMM* candidate=new CGMM(m_components, m_coefficients, true);
298  candidate->train(features);
299  candidate->partial_em(split_ind[i], merge_ind[j]/m_components.vlen, merge_ind[j]%m_components.vlen, min_cov, max_em_iter, min_change);
300  float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change);
301 
302  if (cand_likelihood>cur_likelihood)
303  {
304  cur_likelihood=cand_likelihood;
305  set_comp(candidate->get_comp());
306  set_coef(candidate->get_coef());
307 
308  for (int32_t k=0; k<candidate->get_comp().vlen; k++)
309  {
310  SG_UNREF(candidate->get_comp().vector[k]);
311  }
312 
313  better_found=true;
314  break;
315  }
316  else
317  delete candidate;
318 
319  if (candidates_checked>=max_cand)
320  break;
321  }
322  }
323  if (better_found || candidates_checked>=max_cand)
324  break;
325  }
326  if (!better_found)
327  break;
328  iter++;
329  }
330 
331  SG_FREE(logPxy);
332  SG_FREE(logPx);
333  SG_FREE(logPost);
334  SG_FREE(split_crit);
335  SG_FREE(merge_crit);
336  SG_FREE(logPostSum);
337  SG_FREE(logPostSum2);
338  SG_FREE(logPostSumSum);
339  SG_FREE(split_ind);
340  SG_FREE(merge_ind);
341 
342  return cur_likelihood;
343 }
344 
345 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)
346 {
347  CDotFeatures* dotdata=(CDotFeatures *) features;
348  int32_t num_vectors=dotdata->get_num_vectors();
349 
350  float64_t* init_logPxy=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
351  float64_t* init_logPx=SG_MALLOC(float64_t, num_vectors);
352  float64_t* init_logPx_fix=SG_MALLOC(float64_t, num_vectors);
353  float64_t* post_add=SG_MALLOC(float64_t, num_vectors);
354 
355  for (int32_t i=0; i<num_vectors; i++)
356  {
357  init_logPx[i]=0;
358  init_logPx_fix[i]=0;
359 
361  for (int32_t j=0; j<m_components.vlen; j++)
362  {
364  init_logPx[i]+=CMath::exp(init_logPxy[i*m_components.vlen+j]);
365  if (j!=comp1 && j!=comp2 && j!=comp3)
366  {
367  init_logPx_fix[i]+=CMath::exp(init_logPxy[i*m_components.vlen+j]);
368  }
369  }
370 
371  init_logPx[i]=CMath::log(init_logPx[i]);
372  post_add[i]=CMath::log(CMath::exp(init_logPxy[i*m_components.vlen+comp1]-init_logPx[i])+
373  CMath::exp(init_logPxy[i*m_components.vlen+comp2]-init_logPx[i])+
374  CMath::exp(init_logPxy[i*m_components.vlen+comp3]-init_logPx[i]));
375  v.free_vector();
376  }
377 
378  SGVector<CGaussian*> components(3);
379  SGVector<float64_t> coefficients(3);
380  components.vector[0]=m_components.vector[comp1];
381  components.vector[1]=m_components.vector[comp2];
382  components.vector[2]=m_components.vector[comp3];
383  coefficients.vector[0]=m_coefficients.vector[comp1];
384  coefficients.vector[1]=m_coefficients.vector[comp2];
385  coefficients.vector[2]=m_coefficients.vector[comp3];
386  float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
387 
388  int32_t dim_n=components.vector[0]->get_mean().vlen;
389 
390  float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
391  float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
392 
393  float64_t noise_mag=CMath::twonorm(components.vector[0]->get_mean().vector, dim_n)*0.1/
394  CMath::sqrt((float64_t)dim_n);
395 
396  CMath::add(components.vector[1]->get_mean().vector, alpha1, components.vector[1]->get_mean().vector, alpha2,
397  components.vector[2]->get_mean().vector, dim_n);
398 
399  for (int32_t i=0; i<dim_n; i++)
400  {
401  components.vector[2]->get_mean().vector[i]=components.vector[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
402  components.vector[0]->get_mean().vector[i]=components.vector[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
403  }
404 
405  coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
406  coefficients.vector[2]=coefficients.vector[0]*0.5;
407  coefficients.vector[0]=coefficients.vector[0]*0.5;
408 
409  if (components.vector[0]->get_cov_type()==FULL)
410  {
411  SGMatrix<float64_t> c1=components.vector[1]->get_cov();
412  SGMatrix<float64_t> c2=components.vector[2]->get_cov();
413  CMath::add(c1.matrix, alpha1, c1.matrix, alpha2, c2.matrix, dim_n*dim_n);
414 
415  components.vector[1]->set_d(SGVector<float64_t>(CMath::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n));
416  components.vector[1]->set_u(c1);
417 
418  c2.destroy_matrix();
419 
420  float64_t new_d=0;
421  for (int32_t i=0; i<dim_n; i++)
422  {
423  new_d+=CMath::log(components.vector[0]->get_d().vector[i]);
424  for (int32_t j=0; j<dim_n; j++)
425  {
426  if (i==j)
427  {
428  components.vector[0]->get_u().matrix[i*dim_n+j]=1;
429  components.vector[2]->get_u().matrix[i*dim_n+j]=1;
430  }
431  else
432  {
433  components.vector[0]->get_u().matrix[i*dim_n+j]=0;
434  components.vector[2]->get_u().matrix[i*dim_n+j]=0;
435  }
436  }
437  }
438  new_d=CMath::exp(new_d*(1./dim_n));
439  for (int32_t i=0; i<dim_n; i++)
440  {
441  components.vector[0]->get_d().vector[i]=new_d;
442  components.vector[2]->get_d().vector[i]=new_d;
443  }
444  }
445  else if(components.vector[0]->get_cov_type()==DIAG)
446  {
447  CMath::add(components.vector[1]->get_d().vector, alpha1, components.vector[1]->get_d().vector,
448  alpha2, components.vector[2]->get_d().vector, dim_n);
449 
450  float64_t new_d=0;
451  for (int32_t i=0; i<dim_n; i++)
452  {
453  new_d+=CMath::log(components.vector[0]->get_d().vector[i]);
454  }
455  new_d=CMath::exp(new_d*(1./dim_n));
456  for (int32_t i=0; i<dim_n; i++)
457  {
458  components.vector[0]->get_d().vector[i]=new_d;
459  components.vector[2]->get_d().vector[i]=new_d;
460  }
461  }
462  else if(components.vector[0]->get_cov_type()==SPHERICAL)
463  {
464  components.vector[1]->get_d().vector[0]=alpha1*components.vector[1]->get_d().vector[0]+
465  alpha2*components.vector[2]->get_d().vector[0];
466 
467  components.vector[2]->get_d().vector[0]=components.vector[0]->get_d().vector[0];
468  }
469 
470  CGMM* partial_candidate=new CGMM(components, coefficients);
471  partial_candidate->train(features);
472 
473  float64_t log_likelihood_prev=0;
474  float64_t log_likelihood_cur=0;
475  int32_t iter=0;
476  SGMatrix<float64_t> alpha(num_vectors, 3);
477  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*3);
478  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
479  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
480 
481  while (iter<max_em_iter)
482  {
483  log_likelihood_prev=log_likelihood_cur;
484  log_likelihood_cur=0;
485 
486  for (int32_t i=0; i<num_vectors; i++)
487  {
488  logPx[i]=0;
490  for (int32_t j=0; j<3; j++)
491  {
492  logPxy[i*3+j]=components.vector[j]->compute_log_PDF(v)+CMath::log(coefficients.vector[j]);
493  logPx[i]+=CMath::exp(logPxy[i*3+j]);
494  }
495 
496  logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]);
497  log_likelihood_cur+=logPx[i];
498  v.free_vector();
499 
500  for (int32_t j=0; j<3; j++)
501  {
502  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
503  alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
504  }
505  }
506 
507  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
508  break;
509 
510  partial_candidate->max_likelihood(alpha, min_cov);
511  partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum;
512  partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum;
513  partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum;
514  iter++;
515  }
516 
517  m_coefficients.vector[comp1]=coefficients.vector[0];
518  m_coefficients.vector[comp2]=coefficients.vector[1];
519  m_coefficients.vector[comp3]=coefficients.vector[2];
520 
521  delete partial_candidate;
522  alpha.destroy_matrix();
523  SG_FREE(logPxy);
524  SG_FREE(logPx);
525  SG_FREE(init_logPxy);
526  SG_FREE(init_logPx);
527  SG_FREE(init_logPx_fix);
528  SG_FREE(post_add);
529 }
530 
532 {
533  CDotFeatures* dotdata=(CDotFeatures *) features;
534  int32_t num_dim=dotdata->get_dim_feature_space();
535 
536  float64_t alpha_sum;
537  float64_t alpha_sum_sum=0;
538  float64_t* mean_sum;
539  float64_t* cov_sum=NULL;
540 
541  for (int32_t i=0; i<alpha.num_cols; i++)
542  {
543  alpha_sum=0;
544  mean_sum=SG_MALLOC(float64_t, num_dim);
545  memset(mean_sum, 0, num_dim*sizeof(float64_t));
546 
547  for (int32_t j=0; j<alpha.num_rows; j++)
548  {
549  alpha_sum+=alpha.matrix[j*alpha.num_cols+i];
551  CMath::add<float64_t>(mean_sum, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, mean_sum, v.vlen);
552  v.free_vector();
553  }
554 
555  for (int32_t j=0; j<num_dim; j++)
556  mean_sum[j]/=alpha_sum;
557 
558  m_components.vector[i]->set_mean(SGVector<float64_t>(mean_sum, num_dim));
559 
560  ECovType cov_type=m_components.vector[i]->get_cov_type();
561 
562  if (cov_type==FULL)
563  {
564  cov_sum=SG_MALLOC(float64_t, num_dim*num_dim);
565  memset(cov_sum, 0, num_dim*num_dim*sizeof(float64_t));
566  }
567  else if(cov_type==DIAG)
568  {
569  cov_sum=SG_MALLOC(float64_t, num_dim);
570  memset(cov_sum, 0, num_dim*sizeof(float64_t));
571  }
572  else if(cov_type==SPHERICAL)
573  {
574  cov_sum=SG_MALLOC(float64_t, 1);
575  cov_sum[0]=0;
576  }
577 
578  for (int32_t j=0; j<alpha.num_rows; j++)
579  {
581  CMath::add<float64_t>(v.vector, 1, v.vector, -1, mean_sum, v.vlen);
582 
583  switch (cov_type)
584  {
585  case FULL:
586  cblas_dger(CblasRowMajor, num_dim, num_dim, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, v.vector,
587  1, (double*) cov_sum, num_dim);
588 
589  break;
590  case DIAG:
591  for (int32_t k=0; k<num_dim; k++)
592  cov_sum[k]+=v.vector[k]*v.vector[k]*alpha.matrix[j*alpha.num_cols+i];
593 
594  break;
595  case SPHERICAL:
596  float64_t temp=0;
597 
598  for (int32_t k=0; k<num_dim; k++)
599  temp+=v.vector[k]*v.vector[k];
600 
601  cov_sum[0]+=temp*alpha.matrix[j*alpha.num_cols+i];
602  break;
603  }
604 
605  v.free_vector();
606  }
607 
608  switch (cov_type)
609  {
610  case FULL:
611  for (int32_t j=0; j<num_dim*num_dim; j++)
612  cov_sum[j]/=alpha_sum;
613 
614  float64_t* d0;
615  d0=CMath::compute_eigenvectors(cov_sum, num_dim, num_dim);
616  for (int32_t j=0; j<num_dim; j++)
617  d0[j]=CMath::max(min_cov, d0[j]);
618 
619  m_components.vector[i]->set_d(SGVector<float64_t>(d0, num_dim));
620  m_components.vector[i]->set_u(SGMatrix<float64_t>(cov_sum, num_dim, num_dim));
621 
622  break;
623  case DIAG:
624  for (int32_t j=0; j<num_dim; j++)
625  {
626  cov_sum[j]/=alpha_sum;
627  cov_sum[j]=CMath::max(min_cov, cov_sum[j]);
628  }
629 
630  m_components.vector[i]->set_d(SGVector<float64_t>(cov_sum, num_dim));
631 
632  break;
633  case SPHERICAL:
634  cov_sum[0]/=alpha_sum*num_dim;
635  cov_sum[0]=CMath::max(min_cov, cov_sum[0]);
636 
637  m_components.vector[i]->set_d(SGVector<float64_t>(cov_sum, 1));
638 
639  break;
640  }
641 
642  m_coefficients.vector[i]=alpha_sum;
643  alpha_sum_sum+=alpha_sum;
644  }
645 
646  for (int32_t i=0; i<alpha.num_cols; i++)
647  m_coefficients.vector[i]/=alpha_sum_sum;
648 }
649 
651 {
652  return 1;
653 }
654 
656 {
657  ASSERT(num_param==1);
658 
659  return CMath::log(m_components.vlen);
660 }
661 
662 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example)
663 {
665  return 0;
666 }
667 
669 {
671  return 1;
672 }
673 
675 {
676  return CMath::exp(get_log_likelihood_example(num_example));
677 }
678 
680 {
681  ASSERT(num<m_components.vlen);
682  return m_components.vector[num]->get_mean();
683 }
684 
686 {
687  ASSERT(num<m_components.vlen);
688  m_components.vector[num]->set_mean(mean);
689 }
690 
692 {
693  ASSERT(num<m_components.vlen);
694  return m_components.vector[num]->get_cov();
695 }
696 
698 {
699  ASSERT(num<m_components.vlen);
700  m_components.vector[num]->set_cov(cov);
701 }
702 
704 {
705  return m_coefficients;
706 }
707 
709 {
711  m_coefficients=coefficients;
712 }
713 
715 {
716  return m_components;
717 }
718 
720 {
721  for (int32_t i=0; i<m_components.vlen; i++)
722  {
724  }
725 
727  m_components=components;
728 
729  for (int32_t i=0; i<m_components.vlen; i++)
730  {
732  }
733 }
734 
735 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means)
736 {
737  CDotFeatures* dotdata=(CDotFeatures *) features;
738  int32_t num_vectors=dotdata->get_num_vectors();
739 
740  SGVector<float64_t> label_num(init_means.num_cols);
741 
742  for (int32_t i=0; i<init_means.num_cols; i++)
743  label_num.vector[i]=i;
744 
745  CKNN* knn=new CKNN(1, new CEuclidianDistance(), new CLabels(label_num));
746  knn->train(new CSimpleFeatures<float64_t>(init_means));
747  CLabels* init_labels=knn->apply(features);
748 
749  SGMatrix<float64_t> alpha(num_vectors, m_components.vlen);
750  memset(alpha.matrix, 0, num_vectors*m_components.vlen*sizeof(float64_t));
751 
752  for (int32_t i=0; i<num_vectors; i++)
753  alpha.matrix[i*m_components.vlen+init_labels->get_int_label(i)]=1;
754 
755  SG_UNREF(init_labels);
756 
757  return alpha;
758 }
759 
761 {
763  float64_t rand_num=CMath::random(float64_t(0), float64_t(1));
764  float64_t cum_sum=0;
765  for (int32_t i=0; i<m_coefficients.vlen; i++)
766  {
767  cum_sum+=m_coefficients.vector[i];
768  if (cum_sum>=rand_num)
769  return m_components.vector[i]->sample();
770  }
772 }
773 
775 {
777  answer.vector[m_components.vlen]=0;
778 
779  for (int32_t i=0; i<m_components.vlen; i++)
780  {
782  answer.vector[m_components.vlen]+=CMath::exp(answer.vector[i]);
783  }
785 
786  return answer;
787 }
788 
789 void CGMM::register_params()
790 {
791  m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components");
792  m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients.");
793 }
794 
795 #endif

SHOGUN Machine Learning Toolbox - Documentation