diff toolboxes/FullBNT-1.0.7/netlab3.3/demgmm1.m @ 0:e9a9cd732c1e tip

first hg version after svn

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolboxes/FullBNT-1.0.7/netlab3.3/demgmm1.m	Tue Feb 10 15:05:51 2015 +0000
@@ -0,0 +1,173 @@
+%DEMGMM1 Demonstrate EM for Gaussian mixtures.
+%
+%	Description
+%	This script demonstrates the use of the EM algorithm to fit a mixture
+%	of Gaussians to a set of data using maximum likelihood. A colour
+%	coding scheme is used to illustrate the evaluation of the posterior
+%	probabilities in the E-step of the EM algorithm.
+%
+%	See also
+%	DEMGMM2, DEMGMM3, DEMGMM4, GMM, GMMEM, GMMPOST
+%
+
+%	Copyright (c) Ian T Nabney (1996-2001)
+
+clc;
+disp('This demonstration illustrates the use of the EM (expectation-')
+disp('maximization) algorithm for fitting of a mixture of Gaussians to a')
+disp('data set by maximum likelihood.')
+disp(' ')
+disp('The data set consists of 40 data points in a 2-dimensional')
+disp('space, generated by sampling from a mixture of 2 Gaussian')
+disp('distributions.')
+disp(' ')
+disp('Press any key to see a plot of the data.')
+pause;
+
+% Generate the data
+randn('state', 0); rand('state', 0);
+gmix = gmm(2, 2, 'spherical');
+ndat1 = 20; ndat2 = 20; ndata = ndat1+ndat2;
+gmix.centres =  [0.3 0.3; 0.7 0.7]; 
+gmix.covars = [0.01 0.01];
+x = gmmsamp(gmix, ndata);
+
+h = figure;
+hd = plot(x(:, 1), x(:, 2), '.g', 'markersize', 30);
+hold on; axis([0 1 0 1]); axis square; set(gca, 'box', 'on');
+ht = text(0.5, 1.05, 'Data', 'horizontalalignment', 'center');
+disp(' ');
+disp('Press any key to continue.')
+pause; clc;
+
+disp('We next create and initialize a mixture model consisting of a mixture')
+disp('of 2 Gaussians having ''spherical'' covariance matrices, using the')
+disp('function GMM. The Gaussian components can be displayed on the same')
+disp('plot as the data by drawing a contour of constant probability density')
+disp('for each component having radius equal to the corresponding standard')
+disp('deviation. Component 1 is coloured red and component 2 is coloured')
+disp('blue.')
+disp(' ')
+disp('Note that a particulary poor choice of initial parameters has been')
+disp('made in order to illustrate more effectively the operation of the')
+disp('EM algorithm.')
+disp(' ')
+disp('Press any key to see the initial configuration of the mixture model.')
+pause;
+
+% Set up mixture model
+ncentres = 2; input_dim = 2;
+mix = gmm(input_dim, ncentres, 'spherical');
+
+% Initialise the mixture model
+mix.centres = [0.2 0.8; 0.8, 0.2];
+mix.covars = [0.01 0.01];
+
+% Plot the initial model
+ncirc = 30; theta = linspace(0, 2*pi, ncirc);
+xs = cos(theta); ys = sin(theta);
+xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+set(ht, 'string', 'Initial Configuration');
+figure(h);
+disp(' ')
+disp('Press any key to continue'); 
+pause; clc;
+
+disp('Now we adapt the parameters of the mixture model iteratively using the')
+disp('EM algorithm. Each cycle of the EM algorithm consists of an E-step')
+disp('followed by an M-step.  We start with the E-step, which involves the')
+disp('evaluation of the posterior probabilities (responsibilities) which the')
+disp('two components have for each of the data points.')
+disp(' ')
+disp('Since we have labelled the two components using the colours red and')
+disp('blue, a convenient way to indicate the value of a posterior')
+disp('probability for a given data point is to colour the point using a')
+disp('scale ranging from pure red (corresponding to a posterior probability')
+disp('of 1.0 for the red component and 0.0 for the blue component) through')
+disp('to pure blue.')
+disp(' ')
+disp('Press any key to see the result of applying the first E-step.')
+pause;
+
+% Initial E-step.
+set(ht, 'string', 'E-step');
+post = gmmpost(mix, x);
+dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
+delete(hd); 
+for i = 1 : ndata
+  hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
+          'marker', '.', 'markersize', 30);
+end
+figure(h);
+
+disp(' ');
+disp('Press any key to continue')
+pause; clc;
+
+disp('Next we perform the corresponding M-step. This involves replacing the')
+disp('centres of the component Gaussians by the corresponding weighted means')
+disp('of the data. Thus the centre of the red component is replaced by the')
+disp('mean of the data set, in which each data point is weighted according to')
+disp('the amount of red ink (corresponding to the responsibility of')
+disp('component 1 for explaining that data point). The variances and mixing')
+disp('proportions of the two components are similarly re-estimated.')
+disp(' ')
+disp('Press any key to see the result of applying the first M-step.')
+pause;
+
+% M-step.
+set(ht, 'string', 'M-step');
+options = foptions; 
+options(14) = 1; % A single iteration
+options(1) = -1; % Switch off all messages, including warning
+mix = gmmem(mix, x, options);
+delete(hc);
+xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+figure(h);
+disp(' ')
+disp('Press any key to continue')
+pause; clc;
+
+disp('We can continue making alternate E and M steps until the changes in')
+disp('the log likelihood at each cycle become sufficiently small.')
+disp(' ')
+disp('Press any key to see an animation of a further 9 EM cycles.')
+pause;
+figure(h);
+
+% Loop over EM iterations.
+numiters = 9;
+for n = 1 : numiters
+
+  set(ht, 'string', 'E-step');
+  post = gmmpost(mix, x);
+  dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
+  delete(hd); 
+  for i = 1 : ndata
+    hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
+                 'marker', '.', 'markersize', 30);
+  end
+  pause(1)
+
+  set(ht, 'string', 'M-step');
+  [mix, options] = gmmem(mix, x, options);
+  fprintf(1, 'Cycle %4d  Error %11.6f\n', n, options(8));
+  delete(hc);
+  xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+  yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+  hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+  hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+  pause(1)
+
+end
+
+disp(' ');
+disp('Press any key to end.')
+pause; clc; close(h); clear all
+