wolffd@0: %DEMGMM5 Demonstrate density modelling with a PPCA mixture model.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	 The problem consists of modelling data generated by a mixture of
wolffd@0: %	three Gaussians in 2 dimensions with a mixture model using full
wolffd@0: %	covariance matrices.  The priors are 0.3, 0.5 and 0.2; the centres
wolffd@0: %	are (2, 3.5), (0, 0) and (0,2); the variances are (0.16, 0.64) axis
wolffd@0: %	aligned, (0.25, 1) rotated by 30 degrees and the identity matrix. The
wolffd@0: %	first figure contains a scatter plot of the data.
wolffd@0: %
wolffd@0: %	A mixture model with three one-dimensional PPCA components is trained
wolffd@0: %	using EM.  The parameter vector is printed before training and after
wolffd@0: %	training.  The parameter vector consists of priors (the column), and
wolffd@0: %	centres (given as (x, y) pairs as the next two columns).
wolffd@0: %
wolffd@0: %	The second figure is a 3 dimensional view of the density function,
wolffd@0: %	while the third shows the axes of the 1-standard deviation ellipses
wolffd@0: %	for the three components of the mixture model together with the one
wolffd@0: %	standard deviation along the principal component of each mixture
wolffd@0: %	model component.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GMM, GMMINIT, GMMEM, GMMPROB, PPCA
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: 
wolffd@0: ndata = 500;
wolffd@0: data = randn(ndata, 2);
wolffd@0: prior = [0.3 0.5 0.2];
wolffd@0: % Mixture model swaps clusters 1 and 3
wolffd@0: datap = [0.2 0.5 0.3];
wolffd@0: datac = [0 2; 0 0; 2 3.5];
wolffd@0: datacov = repmat(eye(2), [1 1 3]);
wolffd@0: data1 = data(1:prior(1)*ndata,:);
wolffd@0: data2 = data(prior(1)*ndata+1:(prior(2)+prior(1))*ndata, :);
wolffd@0: data3 = data((prior(1)+prior(2))*ndata +1:ndata, :);
wolffd@0: 
wolffd@0: % First cluster has axis aligned variance and centre (2, 3.5)
wolffd@0: data1(:, 1) = data1(:, 1)*0.1 + 2.0;
wolffd@0: data1(:, 2) = data1(:, 2)*0.8 + 3.5;
wolffd@0: datacov(:, :, 3) = [0.1*0.1 0; 0 0.8*0.8];
wolffd@0: 
wolffd@0: % Second cluster has variance axes rotated by 30 degrees and centre (0, 0)
wolffd@0: rotn = [cos(pi/6) -sin(pi/6); sin(pi/6) cos(pi/6)];
wolffd@0: data2(:,1) = data2(:, 1)*0.2;
wolffd@0: data2 = data2*rotn;
wolffd@0: datacov(:, :, 2) = rotn' * [0.04 0; 0 1] * rotn;
wolffd@0: 
wolffd@0: % Third cluster is at (0,2)
wolffd@0: data3(:, 2) = data3(:, 2)*0.1;
wolffd@0: data3 = data3 + repmat([0 2], prior(3)*ndata, 1);
wolffd@0: 
wolffd@0: % Put the dataset together again
wolffd@0: data = [data1; data2; data3];
wolffd@0: 
wolffd@0: ndata = 100;			% Number of data points.
wolffd@0: noise = 0.2;			% Standard deviation of noise distribution.
wolffd@0: x = [0:1/(2*(ndata - 1)):0.5]';
wolffd@0: randn('state', 1);
wolffd@0: rand('state', 1);
wolffd@0: t = sin(2*pi*x) + noise*randn(ndata, 1);
wolffd@0: 
wolffd@0: % Fit three one-dimensional PPCA models
wolffd@0: ncentres = 3;
wolffd@0: ppca_dim = 1;
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('This demonstration illustrates the use of a Gaussian mixture model')
wolffd@0: disp('with a probabilistic PCA covariance structure to approximate the')
wolffd@0: disp('unconditional probability density of data in a two-dimensional space.')
wolffd@0: disp('We begin by generating the data from a mixture of three Gaussians and')
wolffd@0: disp('plotting it.')
wolffd@0: disp(' ')
wolffd@0: disp('The first cluster has axis aligned variance and centre (0, 2).')
wolffd@0: disp('The variance parallel to the x-axis is significantly greater')
wolffd@0: disp('than that parallel to the y-axis.')
wolffd@0: disp('The second cluster has variance axes rotated by 30 degrees')
wolffd@0: disp('and centre (0, 0).  The third cluster has significant variance')
wolffd@0: disp('parallel to the y-axis and centre (2, 3.5).')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: 
wolffd@0: fh1 = figure;
wolffd@0: plot(data(:, 1), data(:, 2), 'o')
wolffd@0: set(gca, 'Box', 'on')
wolffd@0: axis equal
wolffd@0: hold on
wolffd@0: 
wolffd@0: mix = gmm(2, ncentres, 'ppca', ppca_dim);
wolffd@0: options = foptions;
wolffd@0: options(14) = 10;
wolffd@0: options(1) = -1;  % Switch off all warnings
wolffd@0: 
wolffd@0: % Just use 10 iterations of k-means in initialisation
wolffd@0: % Initialise the model parameters from the data
wolffd@0: mix = gmminit(mix, data, options);
wolffd@0: disp('The mixture model has three components with 1-dimensional')
wolffd@0: disp('PPCA subspaces.  The model parameters after initialisation using')
wolffd@0: disp('the k-means algorithm are as follows')
wolffd@0: disp('    Priors        Centres')
wolffd@0: disp([mix.priors' mix.centres])
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue')
wolffd@0: pause
wolffd@0: 
wolffd@0: options(1)  = 1;		% Prints out error values.
wolffd@0: options(14) = 30;		% Number of iterations.
wolffd@0: 
wolffd@0: disp('We now train the model using the EM algorithm for up to 30 iterations.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: 
wolffd@0: [mix, options, errlog] = gmmem(mix, data, options);
wolffd@0: disp('The trained model has priors and centres:')
wolffd@0: disp('    Priors        Centres')
wolffd@0: disp([mix.priors' mix.centres])
wolffd@0: 
wolffd@0: % Now plot the result
wolffd@0: for i = 1:ncentres
wolffd@0:   % Plot the PC vectors
wolffd@0:   v = mix.U(:,:,i);
wolffd@0:   start=mix.centres(i,:)-sqrt(mix.lambda(i))*(v');
wolffd@0:   endpt=mix.centres(i,:)+sqrt(mix.lambda(i))*(v');
wolffd@0:   linex = [start(1) endpt(1)];
wolffd@0:   liney = [start(2) endpt(2)];
wolffd@0:   line(linex, liney, 'Color', 'k', 'LineWidth', 3)
wolffd@0:   % Plot ellipses of one standard deviation
wolffd@0:   theta = 0:0.02:2*pi;
wolffd@0:   x = sqrt(mix.lambda(i))*cos(theta);
wolffd@0:   y = sqrt(mix.covars(i))*sin(theta);
wolffd@0:   % Rotate ellipse axes
wolffd@0:   rot_matrix = [v(1) -v(2); v(2) v(1)];
wolffd@0:   ellipse = (rot_matrix*([x; y]))';
wolffd@0:   % Adjust centre
wolffd@0:   ellipse = ellipse + ones(length(theta), 1)*mix.centres(i,:);
wolffd@0:   plot(ellipse(:,1), ellipse(:,2), 'r-')
wolffd@0: end
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to exit')
wolffd@0: pause
wolffd@0: close (fh1);
wolffd@0: clear all;