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