Distribute Sampling Using GMM EM

Warning

Fix Problem Contains problems not fixed from original wiki.

Synopsis

Compute distributions of samples using GMM EM.

Results

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code

C++

#include "itkVector.h"
#include "itkListSample.h"
#include "itkGaussianMixtureModelComponent.h"
#include "itkExpectationMaximizationMixtureModelEstimator.h"
#include "itkNormalVariateGenerator.h"

int
main(int, char *[])
{
  unsigned int numberOfClasses = 2;
  using MeasurementVectorType = itk::Vector<double, 1>;
  using SampleType = itk::Statistics::ListSample<MeasurementVectorType>;
  SampleType::Pointer sample = SampleType::New();

  using NormalGeneratorType = itk::Statistics::NormalVariateGenerator;
  NormalGeneratorType::Pointer normalGenerator = NormalGeneratorType::New();

  normalGenerator->Initialize(101);

  MeasurementVectorType mv;
  double                mean = 100;
  double                standardDeviation = 30;
  for (unsigned int i = 0; i < 10; ++i)
  {
    mv[0] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    std::cout << "m[" << i << "] = " << mv[0] << std::endl;
    sample->PushBack(mv);
  }

  normalGenerator->Initialize(3024);
  mean = 200;
  standardDeviation = 30;
  for (unsigned int i = 0; i < 10; ++i)
  {
    mv[0] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    std::cout << "m[" << i << "] = " << mv[0] << std::endl;
    sample->PushBack(mv);
  }

  using ParametersType = itk::Array<double>;
  ParametersType params1(2);

  std::vector<ParametersType> initialParameters(numberOfClasses);
  params1[0] = 110.0;
  params1[1] = 50.0;
  initialParameters[0] = params1;

  ParametersType params2(2);
  params2[0] = 210.0;
  params2[1] = 50.0;
  initialParameters[1] = params2;

  using ComponentType = itk::Statistics::GaussianMixtureModelComponent<SampleType>;

  std::vector<ComponentType::Pointer> components;
  for (unsigned int i = 0; i < numberOfClasses; i++)
  {
    components.push_back(ComponentType::New());
    components[i]->SetSample(sample);
    components[i]->SetParameters(initialParameters[i]);
  }

  using EstimatorType = itk::Statistics::ExpectationMaximizationMixtureModelEstimator<SampleType>;
  EstimatorType::Pointer estimator = EstimatorType::New();

  estimator->SetSample(sample);
  estimator->SetMaximumIteration(500);

  itk::Array<double> initialProportions(numberOfClasses);
  initialProportions[0] = 0.5;
  initialProportions[1] = 0.5;

  estimator->SetInitialProportions(initialProportions);

  for (unsigned int i = 0; i < numberOfClasses; i++)
  {
    estimator->AddComponent((ComponentType::Superclass *)components[i].GetPointer());
  }

  estimator->Update();

  for (unsigned int i = 0; i < numberOfClasses; i++)
  {
    std::cout << "Cluster[" << i << "]" << std::endl;
    std::cout << "    Parameters:" << std::endl;
    std::cout << "         " << components[i]->GetFullParameters() << std::endl;
    std::cout << "    Proportion: ";
    std::cout << "         " << estimator->GetProportions()[i] << std::endl;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated