2D Gaussian Mixture Model Expectation Maximum#

Synopsis#

2D Gaussian Mixture Model Expectation Maximization.

Results#

Output:

Cluster[0]
    Parameters:
         [101.40933830302448, 99.43004497807948, 1098.5993639665169, -107.16526601343287, -107.16526601343287, 913.9641556669595]
    Proportion:          0.495716
Cluster[1]
    Parameters:
         [196.3354813961237, 195.29542020949035, 991.7367739288584, 84.51759523418217, 84.51759523418217, 845.9604643808337]
    Proportion:          0.504284

Code#

C++#

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

int
main()
{
  unsigned int numberOfClasses = 2;
  using MeasurementVectorType = itk::Vector<double, 2>;
  using SampleType = itk::Statistics::ListSample<MeasurementVectorType>;
  auto sample = SampleType::New();


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

  // Create the first set of 2D Gaussian samples
  normalGenerator->Initialize(101);

  MeasurementVectorType mv;
  double                mean = 100;
  double                standardDeviation = 30;
  for (unsigned int i = 0; i < 100; ++i)
  {
    mv[0] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    mv[1] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    sample->PushBack(mv);
  }

  // Create the second set of 2D Gaussian samples
  normalGenerator->Initialize(3024);
  mean = 200;
  standardDeviation = 30;
  for (unsigned int i = 0; i < 100; ++i)
  {
    mv[0] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    mv[1] = (normalGenerator->GetVariate() * standardDeviation) + mean;
    sample->PushBack(mv);
  }

  using ParametersType = itk::Array<double>;
  ParametersType params(6);

  // Create the first set of initial parameters
  std::vector<ParametersType> initialParameters(numberOfClasses);
  params[0] = 110.0; // mean of dimension 1
  params[1] = 115.0; // mean of dimension 2
  params[2] = 800.0; // covariance(0,0)
  params[3] = 0;     // covariance(0,1)
  params[4] = 0;     // covariance(1,0)
  params[5] = 805.0; // covariance(1,1)
  initialParameters[0] = params;

  // Create the second set of initial parameters
  params[0] = 210.0; // mean of dimension 1
  params[1] = 215.0; // mean of dimension 2
  params[2] = 850.0; // covariance(0,0)
  params[3] = 0;     // covariance(0,1)
  params[4] = 0;     // covariance(1,0)
  params[5] = 855.0; // covariance(1,1)
  initialParameters[1] = params;

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

  // Create the components
  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>;
  auto estimator = EstimatorType::New();

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

  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();

  // Output the results
  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#

template<typename TSample> class ExpectationMaximizationMixtureModelEstimator : public itk::Object

This class generates the parameter estimates for a mixture model using expectation maximization strategy.

The first template argument is the type of the target sample data. This estimator expects one or more mixture model component objects of the classes derived from the MixtureModelComponentBase. The actual component (or module) parameters are updated by each component. Users can think this class as a strategy or a integration point for the EM procedure. The initial proportion (SetInitialProportions), the input sample (SetSample), the mixture model components (AddComponent), and the maximum iteration (SetMaximumIteration) are required. The EM procedure terminates when the current iteration reaches the maximum iteration or the model parameters converge.

Recent API changes: The static const macro to get the length of a measurement vector, MeasurementVectorSize has been removed to allow the length of a measurement vector to be specified at run time. It is now obtained at run time from the sample set as input. Please use the function GetMeasurementVectorSize() to get the length.

See

MixtureModelComponentBase, GaussianMixtureModelComponent

ITK Sphinx Examples:

All ITK Sphinx Examples

2D Gaussian Mixture Model Expectation Maximum

Distribution Of Pixels Using GMM EM

Distribute Sampling Using GMM EM

See itk::Statistics::ExpectationMaximizationMixtureModelEstimator for additional documentation.