Daniel@0: %DEMGTM2 Demonstrate GTM for visualisation.
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %	 This script demonstrates the use of a GTM with  a two-dimensional
Daniel@0: %	latent space to visualise data in a higher dimensional space. This is
Daniel@0: %	done through the use of the mean responsibility and magnification
Daniel@0: %	factors.
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	DEMGTM1, GTM, GTMEM, GTMPOST
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: 
Daniel@0: % Fix seeds for reproducible results
Daniel@0: rand('state', 420);
Daniel@0: randn('state', 420);
Daniel@0: 
Daniel@0: ndata = 300
Daniel@0: clc;
Daniel@0: disp('This demonstration shows how a Generative Topographic Mapping')
Daniel@0: disp('can be used to model and visualise high dimensional data.  The')
Daniel@0: disp('data is generated from a mixture of two spherical Gaussians in')
Daniel@0: dstring = ['four dimensional space. ', num2str(ndata), ...
Daniel@0:       ' data points are generated.'];
Daniel@0: disp(dstring);
Daniel@0: disp(' ');
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: % Create data
Daniel@0: data_dim = 4;
Daniel@0: latent_dim = 2;
Daniel@0: mix = gmm(data_dim, 2, 'spherical');
Daniel@0: mix.centres = [1 1 1 1; 0 0 0 0];
Daniel@0: mix.priors = [0.5 0.5];
Daniel@0: mix.covars = [0.1 0.1];
Daniel@0: 
Daniel@0: [data, labels] = gmmsamp(mix, ndata);
Daniel@0: 
Daniel@0: latent_shape = [15 15];  % Number of latent points in each dimension
Daniel@0: nlatent = prod(latent_shape);  % Number of latent points
Daniel@0: num_rbf_centres = 16;
Daniel@0: 
Daniel@0: clc;
Daniel@0: dstring = ['Next we generate and initialise the GTM.  There are ',...
Daniel@0:       num2str(nlatent), ' latent points'];
Daniel@0: disp(dstring);
Daniel@0: dstring = ['arranged in a square of ', num2str(latent_shape(1)), ...
Daniel@0:       ' points on a side.  There are ', num2str(num_rbf_centres), ...
Daniel@0:       ' centres in the'];
Daniel@0: disp(dstring);
Daniel@0: disp('RBF model, which has Gaussian activation functions.')
Daniel@0: disp(' ')
Daniel@0: disp('Once the model is created, the latent data sample')
Daniel@0: disp('and RBF centres are placed uniformly in the square [-1 1 -1 1].')
Daniel@0: disp('The output weights of the RBF are computed to map the latent');
Daniel@0: disp('space to the two dimensional PCA subspace of the data.');
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.');
Daniel@0: pause;
Daniel@0: 
Daniel@0: % Create and initialise GTM model
Daniel@0: net = gtm(latent_dim, nlatent, data_dim, num_rbf_centres, ...
Daniel@0:    'gaussian', 0.1);
Daniel@0: 
Daniel@0: options = foptions;
Daniel@0: options(1) = -1;
Daniel@0: options(7) = 1;    % Set width factor of RBF
Daniel@0: net = gtminit(net, options, data, 'regular', latent_shape, [4 4]);
Daniel@0: 
Daniel@0: options = foptions;
Daniel@0: options(14) = 30;
Daniel@0: options(1) = 1;
Daniel@0: 
Daniel@0: clc;
Daniel@0: dstring = ['We now train the model with ', num2str(options(14)), ...
Daniel@0:       ' iterations of'];
Daniel@0: disp(dstring)
Daniel@0: disp('the EM algorithm for the GTM.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause;
Daniel@0: 
Daniel@0: [net, options] = gtmem(net, data, options);
Daniel@0: 
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause;
Daniel@0: 
Daniel@0: clc;
Daniel@0: disp('We now visualise the data by plotting, for each data point,');
Daniel@0: disp('the posterior mean and mode (in latent space).  These give');
Daniel@0: disp('a summary of the entire posterior distribution in latent space.')
Daniel@0: disp('The corresponding values are joined by a line to aid the')
Daniel@0: disp('interpretation.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.');
Daniel@0: pause;
Daniel@0: % Plot posterior means
Daniel@0: means = gtmlmean(net, data);
Daniel@0: modes = gtmlmode(net, data);
Daniel@0: PointSize = 12;
Daniel@0: ClassSymbol1 = 'r.';
Daniel@0: ClassSymbol2 = 'b.';
Daniel@0: fh1 = figure;
Daniel@0: hold on;
Daniel@0: title('Visualisation in latent space')
Daniel@0: plot(means((labels==1),1), means(labels==1,2), ...
Daniel@0:   ClassSymbol1, 'MarkerSize', PointSize)
Daniel@0: plot(means((labels>1),1),means(labels>1,2),...
Daniel@0:    ClassSymbol2, 'MarkerSize', PointSize)
Daniel@0: 
Daniel@0: ClassSymbol1 = 'ro';
Daniel@0: ClassSymbol2 = 'bo';
Daniel@0: plot(modes(labels==1,1), modes(labels==1,2), ...
Daniel@0:   ClassSymbol1)
Daniel@0: plot(modes(labels>1,1),modes(labels>1,2),...
Daniel@0:    ClassSymbol2)
Daniel@0: 
Daniel@0: % Join up means and modes
Daniel@0: for n = 1:ndata
Daniel@0:    plot([means(n,1); modes(n,1)], [means(n,2); modes(n,2)], 'g-')
Daniel@0: end
Daniel@0: % Place legend outside data plot
Daniel@0: legend('Mean (class 1)', 'Mean (class 2)', 'Mode (class 1)',...
Daniel@0:    'Mode (class 2)', -1);
Daniel@0: 
Daniel@0: % Display posterior for a data point
Daniel@0: % Choose an interesting one with a large distance between mean and
Daniel@0: % mode
Daniel@0: [distance, point] = max(sum((means-modes).^2, 2));
Daniel@0: resp = gtmpost(net, data(point, :));
Daniel@0: 
Daniel@0: disp(' ')
Daniel@0: disp('For more detailed information, the full posterior distribution')
Daniel@0: disp('(or responsibility) can be plotted in latent space for a')
Daniel@0: disp('single data point.  This point has been chosen as the one')
Daniel@0: disp('with the largest distance between mean and mode.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.');
Daniel@0: pause;
Daniel@0: 
Daniel@0: R = reshape(resp, fliplr(latent_shape));
Daniel@0: XL = reshape(net.X(:,1), fliplr(latent_shape));
Daniel@0: YL = reshape(net.X(:,2), fliplr(latent_shape));
Daniel@0: 
Daniel@0: fh2 = figure;
Daniel@0: imagesc(net.X(:, 1), net.X(:,2), R);
Daniel@0: hold on;
Daniel@0: tstr = ['Responsibility for point ', num2str(point)];
Daniel@0: title(tstr);
Daniel@0: set(gca,'YDir','normal')
Daniel@0: colormap(hot);
Daniel@0: colorbar
Daniel@0: disp(' ');
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('Finally, we visualise the data with the posterior means in')
Daniel@0: disp('latent space as before, but superimpose the magnification')
Daniel@0: disp('factors to highlight the separation between clusters.')
Daniel@0: disp(' ')
Daniel@0: disp('Note the large magnitude factors down the centre of the')
Daniel@0: disp('graph, showing that the manifold is stretched more in')
Daniel@0: disp('this region than within each of the two clusters.')
Daniel@0: ClassSymbol1 = 'g.';
Daniel@0: ClassSymbol2 = 'b.';
Daniel@0: 
Daniel@0: fh3 = figure;
Daniel@0: mags = gtmmag(net, net.X);
Daniel@0: % Reshape into grid form
Daniel@0: Mags = reshape(mags, fliplr(latent_shape));
Daniel@0: imagesc(net.X(:, 1), net.X(:,2), Mags);
Daniel@0: hold on
Daniel@0: title('Dataset visualisation with magnification factors')
Daniel@0: set(gca,'YDir','normal')
Daniel@0: colormap(hot);
Daniel@0: colorbar
Daniel@0: hold on; % Else the magnification plot disappears
Daniel@0: plot(means(labels==1,1), means(labels==1,2), ...
Daniel@0:   ClassSymbol1, 'MarkerSize', PointSize)
Daniel@0: plot(means(labels>1,1), means(labels>1,2), ...
Daniel@0:   ClassSymbol2, 'MarkerSize', PointSize)
Daniel@0: 
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to exit.')
Daniel@0: pause
Daniel@0: 
Daniel@0: close(fh1);
Daniel@0: close(fh2);
Daniel@0: close(fh3);
Daniel@0: clear all;