Daniel@0: %DEMGMM1 Demonstrate density modelling with a Gaussian mixture model.
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %	The problem consists of modelling data generated by a mixture of
Daniel@0: %	three Gaussians in 2 dimensions.  The priors are 0.3, 0.5 and 0.2;
Daniel@0: %	the centres are (2, 3.5), (0, 0) and (0,2); the variances are 0.2,
Daniel@0: %	0.5 and 1.0. The first figure contains a  scatter plot of the data.
Daniel@0: %
Daniel@0: %	A Gaussian mixture model with three components is trained using EM.
Daniel@0: %	The parameter vector is printed before training and after training.
Daniel@0: %	The user should press any key to continue at these points.  The
Daniel@0: %	parameter vector consists of priors (the column), centres (given as
Daniel@0: %	(x, y) pairs as the next two columns), and variances (the last
Daniel@0: %	column).
Daniel@0: %
Daniel@0: %	The second figure is a 3 dimensional view of the density function,
Daniel@0: %	while the third shows the 1-standard deviation circles for the three
Daniel@0: %	components of the mixture model.
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	GMM, GMMINIT, GMMEM, GMMPROB, GMMUNPAK
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: % Generate the data
Daniel@0: % Fix seeds for reproducible results
Daniel@0: randn('state', 42);
Daniel@0: rand('state', 42);
Daniel@0: 
Daniel@0: ndata = 500;
Daniel@0: [data, datac, datap, datasd] = dem2ddat(ndata);
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('This demonstration illustrates the use of a Gaussian mixture model')
Daniel@0: disp('to approximate the unconditional probability density of data in')
Daniel@0: disp('a two-dimensional space.  We begin by generating the data from')
Daniel@0: disp('a mixture of three Gaussians and plotting it.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue')
Daniel@0: pause
Daniel@0: 
Daniel@0: fh1 = figure;
Daniel@0: plot(data(:, 1), data(:, 2), 'o')
Daniel@0: set(gca, 'Box', 'on')
Daniel@0: % Set up mixture model
Daniel@0: ncentres = 3;
Daniel@0: input_dim = 2;
Daniel@0: mix = gmm(input_dim, ncentres, 'spherical');
Daniel@0: 
Daniel@0: options = foptions;
Daniel@0: options(14) = 5;	% Just use 5 iterations of k-means in initialisation
Daniel@0: % Initialise the model parameters from the data
Daniel@0: mix = gmminit(mix, data, options);
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('The data is drawn from a mixture with parameters')
Daniel@0: disp('    Priors        Centres         Variances')
Daniel@0: disp([datap' datac (datasd.^2)'])
Daniel@0: disp(' ')
Daniel@0: disp('The mixture model has three components and spherical covariance')
Daniel@0: disp('matrices.  The model parameters after initialisation using the')
Daniel@0: disp('k-means algorithm are as follows')
Daniel@0: % Print out model
Daniel@0: disp('    Priors        Centres         Variances')
Daniel@0: disp([mix.priors' mix.centres mix.covars'])
Daniel@0: disp('Press any key to continue')
Daniel@0: pause
Daniel@0: 
Daniel@0: % Set up vector of options for EM trainer
Daniel@0: options = zeros(1, 18);
Daniel@0: options(1)  = 1;		% Prints out error values.
Daniel@0: options(14) = 10;		% Max. Number of iterations.
Daniel@0: 
Daniel@0: disp('We now train the model using the EM algorithm for 10 iterations')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue')
Daniel@0: pause
Daniel@0: [mix, options, errlog] = gmmem(mix, data, options);
Daniel@0: 
Daniel@0: % Print out model
Daniel@0: disp(' ')
Daniel@0: disp('The trained model has parameters ')
Daniel@0: disp('    Priors        Centres         Variances')
Daniel@0: disp([mix.priors' mix.centres mix.covars'])
Daniel@0: disp('Note the close correspondence between these parameters and those')
Daniel@0: disp('of the distribution used to generate the data, which are repeated here.')
Daniel@0: disp('    Priors        Centres         Variances')
Daniel@0: disp([datap' datac (datasd.^2)'])
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue')
Daniel@0: pause
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('We now plot the density given by the mixture model as a surface plot')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue')
Daniel@0: pause
Daniel@0: % Plot the result
Daniel@0: x = -4.0:0.2:5.0;
Daniel@0: y = -4.0:0.2:5.0;
Daniel@0: [X, Y] = meshgrid(x,y);
Daniel@0: X = X(:);
Daniel@0: Y = Y(:);
Daniel@0: grid = [X Y];
Daniel@0: Z = gmmprob(mix, grid);
Daniel@0: Z = reshape(Z, length(x), length(y));
Daniel@0: c = mesh(x, y, Z);
Daniel@0: hold on
Daniel@0: title('Surface plot of probability density')
Daniel@0: hold off
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('The final plot shows the centres and widths, given by one standard')
Daniel@0: disp('deviation, of the three components of the mixture model.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: % Try to calculate a sensible position for the second figure, below the first
Daniel@0: fig1_pos = get(fh1, 'Position');
Daniel@0: fig2_pos = fig1_pos;
Daniel@0: fig2_pos(2) = fig2_pos(2) - fig1_pos(4);
Daniel@0: fh2 = figure;
Daniel@0: set(fh2, 'Position', fig2_pos)
Daniel@0: 
Daniel@0: hp1 = plot(data(:, 1), data(:, 2), 'bo');
Daniel@0: axis('equal');
Daniel@0: hold on
Daniel@0: hp2 = plot(mix.centres(:, 1), mix.centres(:,2), 'g+');
Daniel@0: set(hp2, 'MarkerSize', 10);
Daniel@0: set(hp2, 'LineWidth', 3);
Daniel@0: 
Daniel@0: title('Plot of data and mixture centres')
Daniel@0: angles = 0:pi/30:2*pi;
Daniel@0: for i = 1 : mix.ncentres
Daniel@0:   x_circle = mix.centres(i,1)*ones(1, length(angles)) + ...
Daniel@0:     sqrt(mix.covars(i))*cos(angles);
Daniel@0:   y_circle = mix.centres(i,2)*ones(1, length(angles)) + ...
Daniel@0:     sqrt(mix.covars(i))*sin(angles);
Daniel@0:   plot(x_circle, y_circle, 'r')
Daniel@0: end
Daniel@0: hold off
Daniel@0: disp('Note how the data cluster positions and widths are captured by')
Daniel@0: disp('the mixture model.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to end.')
Daniel@0: pause
Daniel@0: 
Daniel@0: close(fh1);
Daniel@0: close(fh2);
Daniel@0: clear all; 
Daniel@0: