| 482 |
482 |
|
| 483 |
483 |
// probability = (2*pi)^(-dims/2)*exp( -0.5 * cur )
|
| 484 |
484 |
cvConvertScale( &expo, &expo, -0.5 );
|
| 485 |
|
double factor = -double(dims)/2.0 * log(2.0*CV_PI);
|
| 486 |
|
cvAndS( &expo, cvScalar(factor), &expo );
|
|
485 |
const double factor = -0.5*double(dims)*CV_LOG2PI;
|
|
486 |
cvAddS( &expo, cvScalar(factor), &expo );
|
| 487 |
487 |
|
| 488 |
488 |
// Calculate the log-likelihood of the given sample -
|
| 489 |
489 |
// see Alex Smola's blog http://blog.smola.org/page/2 for
|
| 490 |
490 |
// details on the log-sum-exp trick
|
| 491 |
|
double mini,maxi,retval;
|
|
491 |
double mini, maxi;
|
| 492 |
492 |
cvMinMaxLoc( &expo, &mini, &maxi, 0, 0 );
|
| 493 |
|
CvMat *flp = cvCloneMat(&expo);
|
| 494 |
|
cvSubS( &expo, cvScalar(maxi), flp);
|
| 495 |
|
cvExp( flp, flp );
|
| 496 |
|
CvScalar ss = cvSum( flp );
|
| 497 |
|
retval = log(ss.val[0]) + maxi;
|
| 498 |
|
cvReleaseMat(&flp);
|
|
493 |
cvSubS( &expo, cvScalar(maxi), &expo);
|
|
494 |
cvExp( &expo, &expo );
|
|
495 |
const CvScalar ss = cvSum( &expo );
|
| 499 |
496 |
|
| 500 |
497 |
if( sample_data != _sample->data.fl )
|
| 501 |
498 |
cvFree( &sample_data );
|
| 502 |
499 |
|
| 503 |
|
return retval;
|
|
500 |
return log(ss.val[0]) + maxi;
|
| 504 |
501 |
}
|
| 505 |
502 |
|
| 506 |
503 |
/****************************************************************************************/
|
| ... | ... | |
| 561 |
558 |
}
|
| 562 |
559 |
}
|
| 563 |
560 |
|
| 564 |
|
// probability = (2*pi)^(-dims/2)*exp( -0.5 * cur )
|
| 565 |
|
cvConvertScale( &expo, &expo, -0.5 );
|
| 566 |
|
double factor = -double(dims)/2.0 * log(2.0*CV_PI);
|
| 567 |
|
cvAndS( &expo, cvScalar(factor), &expo );
|
| 568 |
|
|
| 569 |
561 |
// Calculate the posterior probability of each component
|
| 570 |
562 |
// given the sample data.
|
| 571 |
563 |
if( _probs )
|
| 572 |
564 |
{
|
| 573 |
|
cvExp( &expo, &expo );
|
| 574 |
|
if( _probs->cols == 1 )
|
| 575 |
|
cvReshape( &expo, &expo, 0, nclusters );
|
| 576 |
|
cvConvertScale( &expo, _probs, 1./cvSum( &expo ).val[0] );
|
|
565 |
cvScale( &expo, &expo, -0.5 );
|
|
566 |
|
|
567 |
// Use the same log-sum-exp trick like in calcLikelihood,
|
|
568 |
// but neither max_i(log(p_i)) must be added back (because of
|
|
569 |
// normalization) nor taking of log of sum is necessary
|
|
570 |
double mini, maxi;
|
|
571 |
cvMinMaxLoc( &expo, &mini, &maxi, 0, 0 );
|
|
572 |
cvSubS( &expo, cvScalar(maxi), &expo);
|
|
573 |
cvExp( &expo, &expo );
|
|
574 |
if( _probs->cols == 1 )
|
|
575 |
cvReshape( &expo, &expo, 0, nclusters );
|
|
576 |
cvConvertScale( &expo, _probs, 1./cvSum( &expo ).val[0] );
|
| 577 |
577 |
}
|
| 578 |
578 |
|
| 579 |
579 |
if( sample_data != _sample->data.fl )
|
| ... | ... | |
| 814 |
814 |
{
|
| 815 |
815 |
class_ranges->data.i[0] = 0;
|
| 816 |
816 |
class_ranges->data.i[1] = nsamples;
|
|
817 |
nclusters_found = 1;
|
| 817 |
818 |
}
|
| 818 |
819 |
|
| 819 |
820 |
for( i = 0; i < nclusters_found; i++ )
|
| ... | ... | |
| 1081 |
1082 |
|
| 1082 |
1083 |
int nsamples = train_data.count, dims = train_data.dims, nclusters = params.nclusters;
|
| 1083 |
1084 |
double min_variation = FLT_EPSILON;
|
| 1084 |
|
double min_det_value = MAX( DBL_MIN, pow( min_variation, dims ));
|
|
1085 |
//double min_det_value = MAX( DBL_MIN, pow( min_variation, dims ));
|
|
1086 |
double min_log_det_value = dims*std::log(min_variation);
|
| 1085 |
1087 |
double _log_likelihood = -DBL_MAX;
|
| 1086 |
1088 |
int start_step = params.start_step;
|
| 1087 |
1089 |
double sum_max_val;
|
| ... | ... | |
| 1095 |
1097 |
|
| 1096 |
1098 |
if( nclusters == 1 )
|
| 1097 |
1099 |
{
|
| 1098 |
|
double log_weight;
|
| 1099 |
1100 |
CV_CALL( cvSet( probs, cvScalar(1.)) );
|
| 1100 |
1101 |
|
| 1101 |
1102 |
if( params.cov_mat_type == COV_MAT_SPHERICAL )
|
| ... | ... | |
| 1103 |
1104 |
d = cvTrace(*covs).val[0]/dims;
|
| 1104 |
1105 |
d = MAX( d, FLT_EPSILON );
|
| 1105 |
1106 |
inv_eigen_values->data.db[0] = 1./d;
|
| 1106 |
|
log_weight = pow( d, dims*0.5 );
|
|
1107 |
det = dims * std::log(d);
|
| 1107 |
1108 |
}
|
| 1108 |
1109 |
else
|
| 1109 |
1110 |
{
|
| ... | ... | |
| 1115 |
1116 |
cvTranspose( cvGetDiag(*covs, &diag), inv_eigen_values );
|
| 1116 |
1117 |
|
| 1117 |
1118 |
cvMaxS( inv_eigen_values, FLT_EPSILON, inv_eigen_values );
|
| 1118 |
|
for( j = 0, det = 1.; j < dims; j++ )
|
| 1119 |
|
det *= w_data[j];
|
| 1120 |
|
log_weight = sqrt(det);
|
|
1119 |
for( j = 0, det = 0.; j < dims; j++ )
|
|
1120 |
det += std::log( w_data[j] );
|
| 1121 |
1121 |
cvDiv( 0, inv_eigen_values, inv_eigen_values );
|
| 1122 |
1122 |
}
|
| 1123 |
1123 |
|
| 1124 |
|
log_weight_div_det->data.db[0] = -2*log(weights->data.db[0]/log_weight);
|
|
1124 |
log_weight_div_det->data.db[0] = -2*std::log(weights->data.db[0]) + det;
|
| 1125 |
1125 |
log_likelihood = DBL_MAX/1000.;
|
| 1126 |
1126 |
EXIT;
|
| 1127 |
1127 |
}
|
| ... | ... | |
| 1170 |
1170 |
w_data = w->data.db;
|
| 1171 |
1171 |
for( j = 0, det = 0.; j < dims; j++ )
|
| 1172 |
1172 |
det += std::log(w_data[j]);
|
| 1173 |
|
if( det < std::log(min_det_value) )
|
| 1174 |
|
{
|
| 1175 |
|
if( start_step == START_AUTO_STEP )
|
| 1176 |
|
det = std::log(min_det_value);
|
| 1177 |
|
else
|
| 1178 |
|
EXIT;
|
| 1179 |
|
}
|
|
1173 |
if( det < min_log_det_value )
|
|
1174 |
// {
|
|
1175 |
// if( start_step == START_AUTO_STEP )
|
|
1176 |
det = min_log_det_value;
|
|
1177 |
// else
|
|
1178 |
// EXIT;
|
|
1179 |
// early exit at this point _silently_ leaves inv_eigen_values
|
|
1180 |
// uninitialized and log_weight_div_det with wrong
|
|
1181 |
// values; it seems to be a better solution to just let
|
|
1182 |
// it running, i.e. use pseudoinverce on the first step
|
|
1183 |
// and recompute covariances (which hopefully will be
|
|
1184 |
// nonsingular) and then proceed as usual
|
|
1185 |
// }
|
| 1180 |
1186 |
log_det->data.db[k] = det;
|
| 1181 |
1187 |
}
|
| 1182 |
1188 |
else // spherical
|
| 1183 |
1189 |
{
|
| 1184 |
1190 |
d = cvTrace(covs[k]).val[0]/(double)dims;
|
| 1185 |
1191 |
if( d < min_variation )
|
| 1186 |
|
{
|
| 1187 |
|
if( start_step == START_AUTO_STEP )
|
|
1192 |
// {
|
|
1193 |
// if( start_step == START_AUTO_STEP )
|
| 1188 |
1194 |
d = min_variation;
|
| 1189 |
|
else
|
| 1190 |
|
EXIT;
|
| 1191 |
|
}
|
|
1195 |
// else
|
|
1196 |
// EXIT;
|
|
1197 |
// }
|
| 1192 |
1198 |
cov_eigen_values->data.db[k] = d;
|
| 1193 |
|
log_det->data.db[k] = d;
|
|
1199 |
log_det->data.db[k] = dims * std::log(d); // = d;
|
| 1194 |
1200 |
}
|
| 1195 |
1201 |
}
|
| 1196 |
1202 |
|
| 1197 |
|
if( is_spherical )
|
| 1198 |
|
{
|
| 1199 |
|
cvLog( log_det, log_det );
|
| 1200 |
|
cvScale( log_det, log_det, dims );
|
| 1201 |
|
}
|
|
1203 |
// if( is_spherical )
|
|
1204 |
// {
|
|
1205 |
// cvLog( log_det, log_det );
|
|
1206 |
// cvScale( log_det, log_det, dims );
|
|
1207 |
// }
|
| 1202 |
1208 |
}
|
| 1203 |
1209 |
|
| 1204 |
1210 |
for( n = 0; n < params.term_crit.max_iter; n++ )
|
| ... | ... | |
| 1323 |
1329 |
d = w_data[0]/(double)dims;
|
| 1324 |
1330 |
d = MAX( d, min_variation );
|
| 1325 |
1331 |
w->data.db[0] = d;
|
| 1326 |
|
log_det->data.db[k] = d;
|
|
1332 |
log_det->data.db[k] = dims * std::log(d); // = d;
|
| 1327 |
1333 |
}
|
| 1328 |
1334 |
else
|
| 1329 |
1335 |
{
|
| ... | ... | |
| 1340 |
1346 |
cvConvertScale( weights, weights, 1./(double)nsamples, 0 );
|
| 1341 |
1347 |
cvMaxS( weights, DBL_MIN, weights );
|
| 1342 |
1348 |
|
| 1343 |
|
if( is_spherical )
|
| 1344 |
|
{
|
| 1345 |
|
cvLog( log_det, log_det );
|
| 1346 |
|
cvScale( log_det, log_det, dims );
|
| 1347 |
|
}
|
|
1349 |
// if( is_spherical )
|
|
1350 |
// {
|
|
1351 |
// cvLog( log_det, log_det );
|
|
1352 |
// cvScale( log_det, log_det, dims );
|
|
1353 |
// }
|
| 1348 |
1354 |
} // end of iteration process
|
| 1349 |
1355 |
|
| 1350 |
1356 |
//log_weight_div_det[k] = -2*log(weights_k/det(Sigma_k))^0.5) = -2*log(weights_k) + log(det(Sigma_k)))
|