wolffd@0: %DEMHMC1 Demonstrate Hybrid Monte Carlo sampling on mixture of two Gaussians.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	The problem consists of generating data from a mixture of two
wolffd@0: %	Gaussians in two dimensions using a hybrid Monte Carlo algorithm with
wolffd@0: %	persistence. A mixture model is then fitted to the sample to compare
wolffd@0: %	it with the  true underlying generator.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	DEMHMC3, HMC, DEMPOT, DEMGPOT
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: 
wolffd@0: dim = 2;            	% Data dimension
wolffd@0: ncentres = 2;		% Number of centres in mixture model
wolffd@0: 
wolffd@0: seed = 42;              % Seed for random weight initialization.
wolffd@0: randn('state', seed);
wolffd@0: rand('state', seed);
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('This demonstration illustrates the use of the hybrid Monte Carlo')
wolffd@0: disp('algorithm to sample from a mixture of two Gaussians.')
wolffd@0: disp('The means of the two components are [0 0] and [2 2].')
wolffd@0: disp(' ')
wolffd@0: disp('First we set up the parameters of the mixture model we are sampling')
wolffd@0: disp('from.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: 
wolffd@0: % Set up mixture model to sample from
wolffd@0: mix = gmm(dim, ncentres, 'spherical');
wolffd@0: mix.centres(1, :) = [0 0];
wolffd@0: mix.centres(2, :) = [2 2];
wolffd@0: x = [0 1];  % Start vector
wolffd@0: 
wolffd@0: % Set up vector of options for hybrid Monte Carlo.
wolffd@0: 
wolffd@0: nsamples = 160;		% Number of retained samples.
wolffd@0: 
wolffd@0: options = foptions;     % Default options vector.
wolffd@0: options(1) = 1;		% Switch on diagnostics.
wolffd@0: options(5) = 1;		% Use persistence
wolffd@0: options(7) = 50;	% Number of steps in trajectory.
wolffd@0: options(14) = nsamples;	% Number of Monte Carlo samples returned. 
wolffd@0: options(15) = 30;	% Number of samples omitted at start of chain.
wolffd@0: options(18) = 0.02;
wolffd@0: 
wolffd@0: clc
wolffd@0: disp(['Next we take ', num2str(nsamples),' samples from the distribution.'...
wolffd@0:     , 'The first ', num2str(options(15))])
wolffd@0: disp('samples at the start of the chain are omitted.  As persistence')
wolffd@0: disp('is used, the momentum has a small random component added at each step.')
wolffd@0: disp([num2str(options(7)), ...
wolffd@0:     ' iterations are used at each step and the step size is ',...
wolffd@0:     num2str(options(18))])
wolffd@0: disp('Sampling starts at the point [0 1].')
wolffd@0: disp('The new state is accepted if the threshold value is greater than')
wolffd@0: disp('a random number between 0 and 1.')
wolffd@0: disp(' ')
wolffd@0: disp('Negative step numbers indicate samples discarded from the start of the')
wolffd@0: disp('chain.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: 
wolffd@0: [samples, energies] = hmc('dempot', x, options, 'demgpot', mix);
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: clc
wolffd@0: disp('The plot shows the samples generated by the HMC function.')
wolffd@0: disp('The different colours are used to show how the samples move from')
wolffd@0: disp('one component to the other over time.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: probs = exp(-energies);
wolffd@0: fh1 = figure;
wolffd@0: % Plot data in 4 groups
wolffd@0: ngroups = 4;
wolffd@0: g1end = floor(nsamples/ngroups);
wolffd@0: g2end = floor(2*nsamples/ngroups);
wolffd@0: g3end = floor(3*nsamples/ngroups);
wolffd@0: p1 = plot(samples(1:g1end,1), samples(1:g1end,2), 'k.', 'MarkerSize', 12);
wolffd@0: hold on
wolffd@0: lstrings = char(['Samples 1-' int2str(g1end)], ...
wolffd@0:   ['Samples ' int2str(g1end+1) '-' int2str(g2end)], ...
wolffd@0:   ['Samples ' int2str(g2end+1) '-' int2str(g3end)], ...
wolffd@0:   ['Samples ' int2str(g3end+1) '-' int2str(nsamples)]);
wolffd@0: p2 = plot(samples(g1end+1:g2end,1), samples(g1end+1:g2end,2), ...
wolffd@0:   'r.', 'MarkerSize', 12);
wolffd@0: p3 = plot(samples(g2end+1:g3end,1), samples(g2end+1:g3end,2), ...
wolffd@0:   'g.', 'MarkerSize', 12);
wolffd@0: p4 = plot(samples(g3end+1:nsamples,1), samples(g3end+1:nsamples,2), ...
wolffd@0:   'b.', 'MarkerSize', 12);
wolffd@0: legend([p1 p2 p3 p4], lstrings, 2);
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('We now fit a Gaussian mixture model to the sampled data.')
wolffd@0: disp('The model has spherical covariance structure and the correct')
wolffd@0: disp('number of components.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: % Fit a mixture model to the sample
wolffd@0: newmix = gmm(dim, ncentres, 'spherical');
wolffd@0: options = foptions;
wolffd@0: options(1) = -1;	% Switch off all diagnostics
wolffd@0: options(14) = 5;	% Just use 5 iterations of k-means in initialisation
wolffd@0: % Initialise the model parameters from the samples
wolffd@0: newmix = gmminit(newmix, samples, options);
wolffd@0: 
wolffd@0: % Set up vector of options for EM trainer
wolffd@0: options = zeros(1, 18);
wolffd@0: options(1)  = 1;		% Prints out error values.
wolffd@0: options(14) = 15;		% Max. Number of iterations.
wolffd@0: 
wolffd@0: disp('We now train the model using the EM algorithm for 15 iterations')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue')
wolffd@0: pause
wolffd@0: [newmix, options, errlog] = gmmem(newmix, samples, options);
wolffd@0: 
wolffd@0: % Print out model
wolffd@0: disp(' ')
wolffd@0: disp('The trained model has parameters ')
wolffd@0: disp('    Priors        Centres         Variances')
wolffd@0: disp([newmix.priors' newmix.centres newmix.covars'])
wolffd@0: disp('Note the close correspondence between these parameters and those')
wolffd@0: disp('of the distribution used to generate the data')
wolffd@0: disp(' ')
wolffd@0: disp('    Priors        Centres         Variances')
wolffd@0: disp([mix.priors' mix.centres mix.covars'])
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to exit')
wolffd@0: pause
wolffd@0: 
wolffd@0: close(fh1);
wolffd@0: clear all;
wolffd@0: