wolffd@0: function samples=demmet1(plot_wait)
wolffd@0: %DEMMET1 Demonstrate Markov Chain Monte Carlo sampling on a Gaussian.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	The problem consists of generating data from a Gaussian in two
wolffd@0: %	dimensions using a Markov Chain Monte Carlo algorithm. The points are
wolffd@0: %	plotted one after another to show the path taken by the chain.
wolffd@0: %
wolffd@0: %	DEMMET1(PLOTWAIT) allows the user to set the time (in a whole number
wolffd@0: %	of seconds) between the plotting of points.  This is passed to PAUSE
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	DEMHMC1, METROP, GMM, DEMPOT
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: if nargin == 0 | plot_wait < 0
wolffd@0:   plot_wait = 0; % No wait if not specified or incorrect
wolffd@0: end
wolffd@0: dim = 2;            	% Data dimension
wolffd@0: ncentres = 1;		% 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 Markov chain Monte Carlo')
wolffd@0: disp('algorithm to sample from a Gaussian distribution.')
wolffd@0: disp('The mean is at [0 0].')
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: x = [0 4];  % Start vector
wolffd@0: 
wolffd@0: % Set up vector of options for hybrid Monte Carlo.
wolffd@0: 
wolffd@0: nsamples = 150;		% Number of retained samples.
wolffd@0: 
wolffd@0: options = foptions;     % Default options vector.
wolffd@0: options(1) = 0;		% Switch off diagnostics.
wolffd@0: options(14) = nsamples;	% Number of Monte Carlo samples returned. 
wolffd@0: options(18) = 0.1;
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('Next we take 150 samples from the distribution.')
wolffd@0: disp('Sampling starts at the point [0 4].')
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('Press any key to continue.')
wolffd@0: pause
wolffd@0: 
wolffd@0: [samples, energies] = metrop('dempot', x, options, '', mix);
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('The plot shows the samples generated by the MCMC function in order')
wolffd@0: disp('as an animation to show the path taken by the Markov chain.')
wolffd@0: disp('The different colours are used to show that the first few samples')
wolffd@0: disp('should be discarded as they lie too far from the mean.')
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: g1end = floor(nsamples/4);
wolffd@0: 
wolffd@0: for n = 1:nsamples
wolffd@0:   
wolffd@0:   if n < g1end
wolffd@0:     Marker = 'k.';
wolffd@0:     p1 = plot(samples(n,1), samples(n,2), Marker, ...
wolffd@0:       'EraseMode', 'none', 'MarkerSize', 12);
wolffd@0:     if n == 1
wolffd@0:       axis([-3 5 -2 5])
wolffd@0:     end
wolffd@0:   else
wolffd@0:     Marker = 'r.';
wolffd@0:     p2 = plot(samples(n,1), samples(n,2), Marker, ...
wolffd@0:       'EraseMode', 'none', 'MarkerSize', 12);
wolffd@0:   end
wolffd@0:   hold on
wolffd@0:   drawnow;  % Force drawing immediately
wolffd@0:   pause(plot_wait);
wolffd@0: end
wolffd@0: lstrings = char(['Samples 1-' int2str(g1end)], ...
wolffd@0:   ['Samples ' int2str(g1end+1) '-' int2str(nsamples)]);
wolffd@0: legend([p1 p2], lstrings, 1);
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to exit.')
wolffd@0: pause
wolffd@0: close(fh1);
wolffd@0: clear all;
wolffd@0: