image thumbnail

Free Split and Merge Expectation Maximization for MultiVariate Gaussian Mixture

by

 

18 Jan 2009 (Updated )

FSMEM can estimate MVGM parameters and number of conpounds via split/merge cluster moves

test_fsmem_mvgm

Contents

clear, clc, close all hidden

First example of FSMEM : Fit a MVGM with default parameters

Calling FSMEM with the data Z only, without the knowledge of the number of compounds

% Generate MVGM  %%

d                                       = 2;
N                                       = 2000;
K_true                                  = 5;
clust_spread                            = 0.4;

[M_true , S_true , P_true]              = gene_mvgm(d , K_true , clust_spread);
Z                                       = sample_mvgm(N , M_true , S_true , P_true);
logl_true                               = loglike_mvgm(Z , M_true , S_true , P_true);
mdl_true                                = logl_true - 0.5*log(N)*K_true*((d+1) + d*((d+1))/2);

% Display Orignal data %%

offsetmin                               = [1.1 ; 0.9];
offsetmax                               = [0.9 ; 1.1];

pdf_true                                = pdf_mvgm(Z , M_true , S_true , P_true);
minZ                                    = min(Z , [] , 2);
maxZ                                    = max(Z , [] , 2);
minZ                                    = minZ.*offsetmin((sign(minZ) > 0)+1);
maxZ                                    = maxZ.*offsetmax((sign(maxZ) > 0)+1);

vectx                                   = minZ(1):0.1:maxZ(1);
vecty                                   = minZ(2):0.1:maxZ(2);
[X , Y]                                 = meshgrid(vectx , vecty);
grid                                    = [X(:) , Y(:)]';
pdf_support                             = pdf_mvgm(grid , M_true , S_true , P_true);

figure(1)
set(gcf , 'renderer' , 'opengl');
g = surfc(X , Y , reshape(pdf_support , length(vecty) , length(vectx)));
alpha(g , 0.7);
hold on
h = stem3(Z(1 , :) , Z(2 , :) , pdf_true);
hold off
shading interp
lighting phong
light
legend(h , '\bf{Z}')
title(sprintf('Original data Z, K_{true} = %d, MDL(Z) = %8.3f' , K_true , mdl_true) , 'fontsize' , 12 , 'fontweight' , 'bold');
rotate3d on
drawnow

% FSMEM with unknown number of clusters  %

tic,[Mest , Sest , Pest , logl]          = fsmem_mvgm(Z);,toc;
K_est                                    = size(Mest , 3);


% Display Results %

[x_true , y_true]                        = ndellipse(M_true , S_true);
[x_est , y_est]                          = ndellipse(Mest , Sest);

figure(2)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_est , y_est , 'r' , 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true)] , 'Z' , sprintf('True, K=%d' , K_true),sprintf('FSMEM, K=%d' , K_est))
title(sprintf('MDL(FSMEM) = %8.3f' , logl(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');


disp(sprintf('\nK_true = %d , K_est = %d\n' , K_true , K_est))
disp(sprintf('pause\n'))
Elapsed time is 0.592168 seconds.

K_true = 5 , K_est = 5

pause

Second example of FSMEM : Fit a MVGM with default parameters & compare with the classical EM

Calling FSMEM with the data Z only, without the knowledge of the number of compounds
% Generate MVGM  %

d                                       = 2;
N                                       = 1500;
K_true                                  = 5;
clust_spread                            = 0.5;

[M_true , S_true , P_true]              = gene_mvgm(d , K_true , clust_spread);
Z                                       = sample_mvgm(N , M_true , S_true , P_true);
logl_true                               = loglike_mvgm(Z , M_true , S_true , P_true);
mdl_true                                = logl_true - 0.5*log(N)*K_true*((d+1) + d*((d+1))/2);


% Display Orignal data %

offsetmin                               = [1.1 ; 0.9];
offsetmax                               = [0.9 ; 1.1];

pdf_true                                = pdf_mvgm(Z , M_true , S_true , P_true);
minZ                                    = min(Z , [] , 2);
maxZ                                    = max(Z , [] , 2);
minZ                                    = minZ.*offsetmin((sign(minZ) > 0)+1);
maxZ                                    = maxZ.*offsetmax((sign(maxZ) > 0)+1);

vectx                                   = minZ(1):0.1:maxZ(1);
vecty                                   = minZ(2):0.1:maxZ(2);
[X , Y]                                 = meshgrid(vectx , vecty);
grid                                    = [X(:) , Y(:)]';
pdf_support                             = pdf_mvgm(grid , M_true , S_true , P_true);

figure(3)
set(gcf , 'renderer' , 'opengl');
g = surfc(X , Y , reshape(pdf_support , length(vecty) , length(vectx)));
alpha(g , 0.7);
hold on
h = stem3(Z(1 , :) , Z(2 , :) , pdf_true);
hold off
shading interp
lighting phong
light
legend(h , '\bf{Z}')
title(sprintf('Original data Z, K_{true} = %d, MDL(Z) = %8.3f' , K_true , mdl_true) , 'fontsize' , 12 , 'fontweight' , 'bold');
rotate3d on
drawnow

% FSMEM with unknown number of clusters %

tic,[Mest , Sest , Pest , logl ]         = fsmem_mvgm(Z);,toc;
K_est                                    = size(Mest , 3);


% Full EM with known and fixed number of clusters %
% Calling FSMEM with flag fail_exit = 0 in the options structure run a single Full EM

options.Kini                             = K_true;
options.fail_exit                        = 0.0;

tic,[Mestem , Sestem , Pestem , loglem]  = fsmem_mvgm(Z , options);,toc

% Display Results %

[x_true , y_true]                        = ndellipse(M_true , S_true);
[x_est , y_est]                          = ndellipse(Mest , Sest);
[x_estem , y_estem]                      = ndellipse(Mestem , Sestem);

figure(4)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_est , y_est , 'r' , 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true)] , 'Z' , sprintf('True, K=%d' , K_true),sprintf('FSMEM, K=%d' , K_est))
title(sprintf('MDL(FSMEM) = %8.3f' , logl(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');

figure(5)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_estem , y_estem , 'r' , 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true)] , 'Z' , sprintf('True, K=%d' , K_true), 'EM')
title(sprintf('MDL(EM) = %8.3f' , loglem(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');

figure(6)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_est , y_est , 'r' , x_estem , y_estem , 'g', 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true) , h(2+K_true + K_est)] , 'Z' , sprintf('True, K=%d' , K_true),sprintf('FSMEM, K=%d' , K_est) , 'EM')
title(sprintf('MDL(FSMEM) = %8.3f, MDL(EM) = %8.3f' , logl(1 , end) , loglem(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');
Elapsed time is 1.216768 seconds.
Elapsed time is 0.023764 seconds.

Third example of FSMEM : Fit a Elliptical MVGM with a given initial MVGM

Calling FSMEM with the data Z only, without knowing number of compounds

% Generate MVGM  %


d                                       = 2;
N                                       = 2000;
K_true                                  = 8;

P_true                                  = permute((1/K_true)*ones(1 , K_true) , [1 3 2]);
M_true                                  = permute([1.5 , 1 , 0 , -1 , -1.5 , -1 , 0 , 1 ; 0 1 , 1.5 , 1 , 0 , -1 , -1.5 , -1] , [1 3 2]);
S_true                                  = cat(3 , diag([0.01 , 0.1]) , diag([0.1 , 0.1]) , diag([0.1 , 0.01]) , diag([0.1 , 0.1]) , diag([0.01 , 0.1]) , diag([0.1 , 0.1]) , diag([0.1 , 0.01]) , diag([0.1 , 0.1]));
Z                                       = sample_mvgm(N , M_true , S_true , P_true);
logl_true                               = loglike_mvgm(Z , M_true , S_true , P_true);
mdl_true                                = logl_true - 0.5*log(N)*K_true*((d+1) + d*((d+1))/2);

% Display Orignal data %

offsetmin                               = [1.1 ; 0.9];
offsetmax                               = [0.9 ; 1.1];

pdf_true                                = pdf_mvgm(Z , M_true , S_true , P_true);
minZ                                    = min(Z , [] , 2);
maxZ                                    = max(Z , [] , 2);
minZ                                    = minZ.*offsetmin((sign(minZ) > 0)+1);
maxZ                                    = maxZ.*offsetmax((sign(maxZ) > 0)+1);

vectx                                   = minZ(1):0.1:maxZ(1);
vecty                                   = minZ(2):0.1:maxZ(2);
[X , Y]                                 = meshgrid(vectx , vecty);
grid                                    = [X(:) , Y(:)]';
pdf_support                             = pdf_mvgm(grid , M_true , S_true , P_true);

figure(7)
set(gcf , 'renderer' , 'opengl');
g = surfc(X , Y , reshape(pdf_support , length(vecty) , length(vectx)));
alpha(g , 0.7);
hold on
h = stem3(Z(1 , :) , Z(2 , :) , pdf_true);
hold off
shading interp
lighting phong
light
legend(h , '\bf{Z}')
title(sprintf('Original data Z, K_{true} = %d, MDL(Z) = %8.3f' , K_true , mdl_true) , 'fontsize' , 12 , 'fontweight' , 'bold');
rotate3d on

% FSMEM with unknown number of clusters and fixed initial parameters%
% [Mest , Sest , Pest , logl]    = fsmem_mvgm(Z , [] , M_ini , S_ini , P_ini)
% where K_ini clusters <> K_true are supposed for the initial parameter (M_ini , S_ini , P_ini)
% Covariances by default are supposed to be full


K_ini                                    = 3;
[M_ini , S_ini , P_ini]                  = init_mvgm(Z , K_ini);

tic,[Mest , Sest , Pest , logl]          = fsmem_mvgm(Z , [] , M_ini , S_ini , P_ini);,toc;
K_est                                    = size(Mest , 3);

% Display Results %

[x_true , y_true]                        = ndellipse(M_true , S_true);
[x_est , y_est]                          = ndellipse(Mest , Sest);

figure(8)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_est , y_est , 'r' , 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true)] , 'Z' , sprintf('True, K=%d' , K_true),sprintf('FSMEM, K=%d' , K_est))
title(sprintf('MDL(FSMEM) = %8.3f' , logl(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');

% FSMEM with unknown number of clusters and initial parameters %
% Here covariances are supposed to be elliptical <=> options.covtype   = 1 %
% [Mest , Sest , Pest , logl]   = fsmem_mvgm(Z , options , M_ini , S_ini , P_ini) %

clear options
options.covtype                          = 1;
tic,[Mest1 , Sest1 , Pest1 , logl1]      = fsmem_mvgm(Z , options , M_ini , S_ini , P_ini);,toc;
K_est1                                   = size(Mest1 , 3);


% Display Results %

[x_true , y_true]                        = ndellipse(M_true , S_true);
[x_est1 , y_est1]                        = ndellipse(Mest1 , Sest1);

figure(9)
h                                        = plot(Z(1 , :) , Z(2 , :) , '+' , x_true , y_true , 'k' , x_est1 , y_est1 , 'r' , 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_true)] , 'Z' , sprintf('True, K=%d' , K_true),sprintf('FSMEM, K=%d' , K_est1))
title(sprintf('MDL(FSMEM) = %8.3f' , logl1(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');
Elapsed time is 2.048787 seconds.
Elapsed time is 1.325096 seconds.

Fourth example of FSMEM : Fit spiral data %%

Calling FSMEM with special hyperparameters in options structure (a bit long even with the mex-file)

% Generate Spiral data  %


N                                       = 3000;

options.Kini                            = 20;
options.Kmax                            = 40;
options.maxite_fullem                   = 500;
options.maxite_partialem                = 500;
options.epsi_fullem                     = 1e-7;
options.epsi_partialem                  = 1e-7;
options.maxcands_split                  = 12;
options.maxcands_merge                  = 12;
options.covtype                         = 0;


Z                                       = spiral2d(N);

% Display Orignal data %%


figure(10)
plot(Z(1 , :) , Z(2 , :) , '+', 'linewidth' , 2)
title('Original data Z' , 'fontsize' , 12 , 'fontweight' , 'bold');
drawnow


% FSMEM with unknown number of clusters %
% [Mest , Sest , Pest , logl ]  = fsmem_mvgm(Z , options);

tic,[Mest , Sest , Pest , logl]        = fsmem_mvgm(Z , options);,toc
K_est                                  = size(Pest , 3);

% Display Results %

figure(11)

[xest , yest]                          = ndellipse(Mest , Sest);
[xaxes , yaxes]                        = ndellipse(Mest , Sest , [] , [] , 3);
xaxes(end , :)                         = [];
yaxes(end , :)                         = [];
h                                      = plot(Z(1 , :) , Z(2 , :) , '+' , xest , yest , 'r' , xaxes , yaxes , 'k', 'linewidth' , 2);
legend([h(1), h(2) , h(2+K_est)] , 'Z' , sprintf('Estimated, K=%d' , K_est) , 'PC');
title(sprintf('MDL(FSMEM) = %8.3f' , logl(1 , end)) , 'fontsize' , 12 , 'fontweight' , 'bold');
Elapsed time is 52.220716 seconds.

Contact us