34#include "discretepdf.h"
35#include "../wrappers/matrix/vector_wrapper.h"
36#include "../wrappers/matrix/matrix_wrapper.h"
37#include "../wrappers/rng/rng.h"
48 template <
typename T>
class Mixture :
public Pdf<T>
52 unsigned int _numComponents;
55 vector<Probability> *_componentWeights;
58 vector<Pdf<T>* > *_componentPdfs;
61 bool NormalizeWeights();
66 vector<double> _cumWeights;
69 bool CumWeightsUpdate();
73 void TestNotInit()
const;
79 Mixture(
const unsigned int dimension=0);
84 template <
typename U> Mixture(
const U &componentVector);
88 Mixture(
const Mixture & my_mixture);
94 virtual Mixture* Clone()
const;
97 unsigned int NumComponentsGet()
const;
100 Probability ProbabilityGet(
const T& state)
const;
102 bool SampleFrom (vector<Sample<T> >& list_samples,
103 const unsigned int num_samples,
104 const SampleMthd method = SampleMthd::DEFAULT,
105 void * args = NULL)
const;
107 bool SampleFrom (Sample<T>& one_sample,
const SampleMthd method = SampleMthd::DEFAULT,
void * args = NULL)
const;
109 T ExpectedValueGet()
const;
111 MatrixWrapper::SymmetricMatrix CovarianceGet ( )
const;
114 vector<Probability> WeightsGet()
const;
120 Probability WeightGet(
unsigned int componentNumber)
const;
128 bool WeightsSet(vector<Probability> & weights);
140 bool WeightSet(
unsigned int componentNumber, Probability w);
144 int MostProbableComponentGet()
const;
151 bool AddComponent(Pdf<T>& pdf );
157 bool AddComponent(Pdf<T>& pdf, Probability w);
163 bool DeleteComponent(
unsigned int componentNumber );
166 vector<Pdf<T>* > ComponentsGet()
const;
172 Pdf<T>* ComponentGet(
unsigned int componentNumber)
const;
182Mixture<T>::Mixture(
const unsigned int dimension):
185 , _cumWeights(_numComponents+1)
188 _componentWeights =
new vector<Probability>(this->NumComponentsGet());
191 _componentPdfs =
new vector< Pdf<T>* >(NumComponentsGet());
192#ifdef __CONSTRUCTOR__
193 cout <<
"Mixture constructor\n";
198template<
typename T>
template <
typename U>
199Mixture<T>::Mixture(
const U &componentVector): Pdf<T>(componentVector[0]->DimensionGet() )
200 , _numComponents(componentVector.size())
203 _componentWeights =
new vector<Probability>(NumComponentsGet());
204 for (
int i=0; i<NumComponentsGet();i++)
206 (*_componentWeights)[i] = (Probability)(1.0/NumComponentsGet());
208 _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
212 _componentPdfs =
new vector< Pdf<T>* >(NumComponentsGet());
214 for (
int i=0; i<NumComponentsGet();i++)
218 (*_componentPdfs)[i] = (componentVector[i])->Clone();
220#ifdef __CONSTRUCTOR__
221 cout <<
"Mixture constructor\n";
226Mixture<T>::Mixture(
const Mixture & my_mixture):Pdf<T>(my_mixture.DimensionGet() )
227 ,_numComponents(my_mixture.NumComponentsGet())
230 _componentWeights =
new vector<Probability>(this->NumComponentsGet());
231 (*_componentWeights) = my_mixture.WeightsGet();
232 _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
236 _componentPdfs =
new vector< Pdf<T>* >(NumComponentsGet());
237 for (
int i=0; i<NumComponentsGet();i++)
239 (*_componentPdfs)[i] = (my_mixture.ComponentGet(i))->Clone();
241#ifdef __CONSTRUCTOR__
242 cout <<
"Mixture copy constructor\n";
247 Mixture<T>::~Mixture()
249#ifdef __CONSTRUCTOR__
250 cout <<
"Mixture destructor\n";
253 delete _componentWeights;
254 for (
int i=0; i<NumComponentsGet();i++)
256 delete (*_componentPdfs)[i];
258 delete _componentPdfs;
262 Mixture<T>* Mixture<T>::Clone()
const
264 return new Mixture(*
this);
268 unsigned int Mixture<T>::NumComponentsGet()
const
270 return _numComponents;
274 Probability Mixture<T>::ProbabilityGet(
const T& state)
const
277 Probability prob(0.0);
278 for (
int i=0; i<NumComponentsGet();i++)
280 prob= prob + (*_componentWeights)[i] * (*_componentPdfs)[i]->ProbabilityGet(state);
286 bool Mixture<T>::SampleFrom (vector<Sample<T> >& list_samples,
287 const unsigned int num_samples,
288 const SampleMthd method,
294 case SampleMthd::DEFAULT:
295 return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
297 case SampleMthd::RIPLEY:
299 list_samples.resize(num_samples);
301 std::vector<double> unif_samples(num_samples);
302 for (
unsigned int i = 0; i < num_samples ; i++)
303 unif_samples[i] = runif();
306 unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
307 double (1.0/num_samples));
309 for (
int i = num_samples-2; i >= 0 ; i--)
310 unif_samples[i] = pow(unif_samples[i],
double (1.0/(i+1))) * unif_samples[i+1];
313 unsigned int index = 0;
314 unsigned int num_states = NumComponentsGet();
315 vector<double>::const_iterator CumPDFit = _cumWeights.begin();
316 typename vector<Sample<T> >::iterator sit = list_samples.begin();
318 for (
unsigned int i = 0; i < num_samples ; i++)
320 while ( unif_samples[i] > *CumPDFit )
323 assert(index <= num_states);
328 (*_componentPdfs)[index-1]->SampleFrom(*sit,method,args);
334 cerr <<
"Mixture::Samplefrom(T, void *): No such sampling method" << endl;
339 bool Mixture<T>::SampleFrom (Sample<T>& one_sample,
const SampleMthd method,
void * args)
const
344 case SampleMthd::DEFAULT:
347 double unif_sample; unif_sample = runif();
349 unsigned int index = 0;
350 while ( unif_sample > _cumWeights[index] )
352 assert(index <= NumComponentsGet());
357 (*_componentPdfs)[index-1]->SampleFrom(one_sample,method,args);
361 cerr <<
"Mixture::Samplefrom(T, void *): No such sampling method"
368 T Mixture<T>::ExpectedValueGet()
const
370 cerr <<
"Mixture ExpectedValueGet: not implemented for the template parameters you use."
371 << endl <<
"Use template specialization as shown in mixture.cpp " << endl;
374 return expectedValue;
378 MatrixWrapper::SymmetricMatrix Mixture<T>::CovarianceGet ( )
const
381 cerr <<
"Mixture CovarianceGet: not implemented since so far I don't believe its usefull"
382 << endl <<
"If you decide to implement is: Use template specialization as shown in mcpdf.cpp " << endl;
385 MatrixWrapper::SymmetricMatrix result;
390 vector<Probability> Mixture<T>::WeightsGet()
const
393 return *_componentWeights;
397 Probability Mixture<T>::WeightGet(
unsigned int componentNumber)
const
400 assert((
int)componentNumber >= 0 && componentNumber < NumComponentsGet());
401 return (*_componentWeights)[componentNumber];
405 bool Mixture<T>::WeightsSet(vector<Probability> & weights)
408 assert(weights.size() == NumComponentsGet());
409 *_componentWeights = weights;
411 return (NormalizeWeights() && CumWeightsUpdate());
415 bool Mixture<T>::WeightSet(
unsigned int componentNumber, Probability weight)
418 assert((
int)componentNumber >= 0 && componentNumber < NumComponentsGet());
419 assert((
double)weight<=1.0);
421 if (NumComponentsGet() == 1)
423 (*_componentWeights)[0] = weight;
430 Probability old_weight = WeightGet(componentNumber);
431 if ((
double)old_weight!=1.0) {
432 double normalization_factor = (1-weight)/(1-old_weight);
433 for (
int i=0; i<NumComponentsGet();i++)
435 (*_componentWeights)[i] = (Probability)( (
double)( (*_componentWeights)[i] )* normalization_factor);
439 for (
int i=0; i<NumComponentsGet();i++)
441 (*_componentWeights)[i] = (Probability)( (1-weight)/(NumComponentsGet()-1));
444 (*_componentWeights)[componentNumber] = weight;
446 return CumWeightsUpdate();
450 int Mixture<T>::MostProbableComponentGet()
const
453 int index_mostProbable= -1;
454 Probability prob_mostProbable= 0.0;
455 for (
int component = 0 ; component < NumComponentsGet() ; component++)
457 if ( (*_componentWeights)[component] > prob_mostProbable)
459 index_mostProbable= component;
460 prob_mostProbable= (*_componentWeights)[component];
463 return index_mostProbable;
467 bool Mixture<T>::AddComponent(Pdf<T>& pdf)
469 if (NumComponentsGet()==0)
470 return AddComponent(pdf, Probability(1.0));
474 (*_componentPdfs).push_back(pdf.Clone() );
476 (*_componentWeights).push_back(Probability(0.0));
477 _cumWeights.push_back(0.0);
479 assert(NumComponentsGet()==(*_componentPdfs).size());
480 assert(NumComponentsGet()==(*_componentWeights).size());
481 assert(NumComponentsGet()+1==_cumWeights.size());
482 return (NormalizeWeights() && CumWeightsUpdate());
487 bool Mixture<T>::AddComponent(Pdf<T>& pdf, Probability w)
489 if (NumComponentsGet()==0 && w!=1.0)
490 return AddComponent(pdf, Probability(1.0));
494 (*_componentPdfs).push_back(pdf.Clone() );
495 (*_componentWeights).push_back(Probability(0.0));
496 _cumWeights.resize(NumComponentsGet()+1);
498 assert(NumComponentsGet()==(*_componentPdfs).size());
499 assert(NumComponentsGet()==(*_componentWeights).size());
500 assert(NumComponentsGet()+1==_cumWeights.size());
501 WeightSet(_numComponents-1,w);
502 return (NormalizeWeights() && CumWeightsUpdate());
507 bool Mixture<T>::DeleteComponent(
unsigned int componentNumber)
510 assert(NumComponentsGet()==(*_componentPdfs).size());
511 assert(NumComponentsGet()==(*_componentWeights).size());
512 assert(NumComponentsGet()+1==_cumWeights.size());
515 assert((
int)componentNumber >= 0 && componentNumber < NumComponentsGet());
517 Pdf<T>* pointer = (*_componentPdfs)[componentNumber];
519 (*_componentPdfs).erase((*_componentPdfs).begin()+componentNumber);
520 (*_componentWeights).erase((*_componentWeights).begin()+componentNumber);
521 _cumWeights.resize(NumComponentsGet()+1);
523 assert(NumComponentsGet()==(*_componentPdfs).size());
524 assert(NumComponentsGet()==(*_componentWeights).size());
525 assert(NumComponentsGet()+1==_cumWeights.size());
526 if(_numComponents==0)
529 return (NormalizeWeights() && CumWeightsUpdate());
533 vector<Pdf<T>*> Mixture<T>::ComponentsGet()
const
536 return _componentPdfs;
540 Pdf<T>* Mixture<T>::ComponentGet(
unsigned int componentNumber)
const
543 return (*_componentPdfs)[componentNumber];
547 void Mixture<T>::TestNotInit()
const
549 if (NumComponentsGet() == 0)
551 cerr <<
"Mixture method called which requires that the number of components is not zero"
552 << endl <<
"Current number of components: " << NumComponentsGet() << endl;
558 bool Mixture<T>::NormalizeWeights()
560 double SumOfWeights = 0.0;
561 for (
unsigned int i = 0; i < NumComponentsGet() ; i++){
562 SumOfWeights += (*_componentWeights)[i];
564 if (SumOfWeights > 0){
565 for (
unsigned int i = 0; i < NumComponentsGet() ; i++){
566 (*_componentWeights)[i] = (Probability)( (
double) ( (*_componentWeights)[i]) /SumOfWeights);
571 cerr <<
"Mixture::NormalizeProbs(): SumOfWeights = " << SumOfWeights << endl;
577 bool Mixture<T>::CumWeightsUpdate()
581 static vector<double>::iterator CumWeightsit;
582 CumWeightsit = _cumWeights.begin();
586 for (
unsigned int i = 0; i < NumComponentsGet(); i++)
590 CumSum += ( (*_componentWeights)[i] );
591 *CumWeightsit = CumSum;
594 assert( (_cumWeights[NumComponentsGet()] >= 1.0 - NUMERIC_PRECISION) &&
595 (_cumWeights[NumComponentsGet()] <= 1.0 + NUMERIC_PRECISION) );
597 _cumWeights[NumComponentsGet()]=1;
603#include "mixture.cpp"