Daniel@0: %DEMGMM4 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 with a mixture model using full
Daniel@0: %	covariance matrices.  The priors are 0.3, 0.5 and 0.2; the centres
Daniel@0: %	are (2, 3.5), (0, 0) and (0,2); the variances are (0.16, 0.64) axis
Daniel@0: %	aligned, (0.25, 1) rotated by 30 degrees and the identity matrix. The
Daniel@0: %	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), and centres (given
Daniel@0: %	as (x, y) pairs as the next two columns).  The covariance matrices
Daniel@0: %	are printed separately.
Daniel@0: %
Daniel@0: %	The second figure is a 3 dimensional view of the density function,
Daniel@0: %	while the third shows the axes of the 1-standard deviation ellipses
Daniel@0: %	for the three 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: 
Daniel@0: % Generate the data
Daniel@0: 
Daniel@0: ndata = 500;
Daniel@0: 
Daniel@0: % Fix the seeds for reproducible results
Daniel@0: randn('state', 42);
Daniel@0: rand('state', 42);
Daniel@0: data = randn(ndata, 2);
Daniel@0: prior = [0.3 0.5 0.2];
Daniel@0: % Mixture model swaps clusters 1 and 3
Daniel@0: datap = [0.2 0.5 0.3];
Daniel@0: datac = [0 2; 0 0; 2 3.5];
Daniel@0: datacov = repmat(eye(2), [1 1 3]);
Daniel@0: data1 = data(1:prior(1)*ndata,:);
Daniel@0: data2 = data(prior(1)*ndata+1:(prior(2)+prior(1))*ndata, :);
Daniel@0: data3 = data((prior(1)+prior(2))*ndata +1:ndata, :);
Daniel@0: 
Daniel@0: % First cluster has axis aligned variance and centre (2, 3.5)
Daniel@0: data1(:, 1) = data1(:, 1)*0.4 + 2.0;
Daniel@0: data1(:, 2) = data1(:, 2)*0.8 + 3.5;
Daniel@0: datacov(:, :, 3) = [0.4*0.4 0; 0 0.8*0.8];
Daniel@0: 
Daniel@0: % Second cluster has variance axes rotated by 30 degrees and centre (0, 0)
Daniel@0: rotn = [cos(pi/6) -sin(pi/6); sin(pi/6) cos(pi/6)];
Daniel@0: data2(:,1) = data2(:, 1)*0.5;
Daniel@0: data2 = data2*rotn;
Daniel@0: datacov(:, :, 2) = rotn' * [0.25 0; 0 1] * rotn;
Daniel@0: 
Daniel@0: % Third cluster is at (0,2)
Daniel@0: data3 = data3 + repmat([0 2], prior(3)*ndata, 1);
Daniel@0: 
Daniel@0: % Put the dataset together again
Daniel@0: data = [data1; data2; data3];
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('This demonstration illustrates the use of a Gaussian mixture model')
Daniel@0: disp('with full covariance matrices to approximate the unconditional ')
Daniel@0: disp('probability density of data in a two-dimensional space.')
Daniel@0: disp('We begin by generating the data from a mixture of three Gaussians and')
Daniel@0: disp('plotting it.')
Daniel@0: disp(' ')
Daniel@0: disp('The first cluster has axis aligned variance and centre (0, 2).')
Daniel@0: disp('The second cluster has variance axes rotated by 30 degrees')
Daniel@0: disp('and centre (0, 0).  The third cluster has unit variance and centre')
Daniel@0: disp('(2, 3.5).')
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: 
Daniel@0: % Set up mixture model
Daniel@0: ncentres = 3;
Daniel@0: input_dim = 2;
Daniel@0: mix = gmm(input_dim, ncentres, 'full');
Daniel@0: 
Daniel@0: % Initialise the model parameters from the data
Daniel@0: options = foptions;
Daniel@0: options(14) = 5;	% Just use 5 iterations of k-means in initialisation
Daniel@0: mix = gmminit(mix, data, options);
Daniel@0: 
Daniel@0: % Print out model
Daniel@0: clc
Daniel@0: disp('The mixture model has three components and full covariance')
Daniel@0: disp('matrices.  The model parameters after initialisation using the')
Daniel@0: disp('k-means algorithm are as follows')
Daniel@0: disp('    Priors        Centres')
Daniel@0: disp([mix.priors' mix.centres])
Daniel@0: disp('Covariance matrices are')
Daniel@0: disp(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) = 50;		% Number of iterations.
Daniel@0: 
Daniel@0: disp('We now train the model using the EM algorithm for 50 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 priors and centres:')
Daniel@0: disp('    Priors        Centres')
Daniel@0: disp([mix.priors' mix.centres])
Daniel@0: disp('The data generator has priors and centres')
Daniel@0: disp('    Priors        Centres')
Daniel@0: disp([datap' datac])
Daniel@0: disp('Model covariance matrices are')
Daniel@0: disp(mix.covars(:, :, 1))
Daniel@0: disp(mix.covars(:, :, 2))
Daniel@0: disp(mix.covars(:, :, 3))
Daniel@0: disp('Data generator covariance matrices are')
Daniel@0: disp(datacov(:, :, 1))
Daniel@0: disp(datacov(:, :, 2))
Daniel@0: disp(datacov(:, :, 3))
Daniel@0: disp('Note the close correspondence between these parameters and those')
Daniel@0: disp('of the distribution used to generate the data.  The match for')
Daniel@0: disp('covariance matrices is not that close, but would be improved with')
Daniel@0: disp('more iterations of the training algorithm.')
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: 
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: drawnow
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.  The axes')
Daniel@0: disp('of the ellipses of constant density are shown.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: 
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) - 30;
Daniel@0: fh2 = figure('Position', fig2_pos);
Daniel@0: 
Daniel@0: h3 = plot(data(:, 1), data(:, 2), 'bo');
Daniel@0: axis equal;
Daniel@0: hold on
Daniel@0: title('Plot of data and covariances')
Daniel@0: for i = 1:ncentres
Daniel@0:   [v,d] = eig(mix.covars(:,:,i));
Daniel@0:   for j = 1:2
Daniel@0:     % Ensure that eigenvector has unit length
Daniel@0:     v(:,j) = v(:,j)/norm(v(:,j));
Daniel@0:     start=mix.centres(i,:)-sqrt(d(j,j))*(v(:,j)');
Daniel@0:     endpt=mix.centres(i,:)+sqrt(d(j,j))*(v(:,j)');
Daniel@0:     linex = [start(1) endpt(1)];
Daniel@0:     liney = [start(2) endpt(2)];
Daniel@0:     line(linex, liney, 'Color', 'k', 'LineWidth', 3)
Daniel@0:   end
Daniel@0:   % Plot ellipses of one standard deviation
Daniel@0:   theta = 0:0.02:2*pi;
Daniel@0:   x = sqrt(d(1,1))*cos(theta);
Daniel@0:   y = sqrt(d(2,2))*sin(theta);
Daniel@0:   % Rotate ellipse axes
Daniel@0:   ellipse = (v*([x; y]))';
Daniel@0:   % Adjust centre
Daniel@0:   ellipse = ellipse + ones(length(theta), 1)*mix.centres(i,:);
Daniel@0:   plot(ellipse(:,1), ellipse(:,2), 'r-');
Daniel@0: end
Daniel@0: hold off
Daniel@0: 
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: